請務必注意 R
程式碼英文字大小寫有區別!
## Page-01
u <- c(0,4)
v <- c(3,0)
w <- c(0,-2)
a <- c(-2,-4)
b <- c(3,0)
c <- c(-6,4)
(tot <- u+v+w+a+b+c)
## [1] -2 2
# 計算小羊從頭到尾移動的效果並印出
O <- c(0,0)
G <- c(-2,2)
(O-G)
## [1] 2 -2
# 計算小羊如何由G點回到原點並印出
D <- c(3,2)
(D+tot)
## [1] 1 4
# 計算小羊若由D點開始出發，最後將位居何處，並印出
## Page-02
u <- c(6,8)
v <- c(-4,14)
(u+v)
## [1] 2 22
(u+1/2*v)
## [1] 4 15
sqrt(u%*%u)
## [,1]
## [1,] 10
# 向量 u 的長度(內積開根號)
sqrt(v%*%v)
## [,1]
## [1,] 14.56022
# 向量 v 的長度(內積開根號)
# ru+sv=c(0,8r+14s) => 6r=4s => r:s = 2:3 = 1:1.5
(1*u+1.5*v)
## [1] 0 29
(w <- 1.5*v)
## [1] -6 21
sqrt(w%*%w)
## [,1]
## [1,] 21.84033
# 調整為大約 10:21.84033 (10: sqrt(477))
A <- c(88,95,92)
B <- c(96,95,91)
C <- c(91,84,88)
D <- c(85,90,89)
E <- c(95,88,91)
(A+B+C+D+E)
## [1] 455 452 451
# 總分(甲錄取)
(3*A+2*B+C+D+E)
## [1] 727 737 726
# 總裁加權後之總分(乙錄取)
## Page-03
A <- c(0,0,0)
B <- c(2,6,5)
C <- c(4,5,3)
D <- c(5,5,5)
u <- B-A # 或直接 u <- c(2,6,5)
v <- C-B # 或直接 v <- c(2,-1,-2)
w <- D-C # 或直接 w <- c(1,0,2)
(tot <- sqrt(u%*%u)+ sqrt(v%*%v)+ sqrt(w%*%w)) # 總距離
## [,1]
## [1,] 13.29833
s <- u+v+w
(sa<-sqrt(s%*%s)) # 直接飛到D點的飛行距離
## [,1]
## [1,] 8.660254
tot-sa
## [,1]
## [1,] 4.638072
# 節省的距離
uw <- c(0.88,0.75,0.94,0.84,0.93)
uc <- c(0.92,0.80,0.79,0.82,0.95)
(sqrt(uw%*%uw))
## [,1]
## [1,] 1.947049
# 老王的幸福向量長度
(sqrt(uc%*%uc))
## [,1]
## [1,] 1.91974
# 老張的幸福向量長度
up <- c(1,1,1,1,1) # 幸福最高境界
dw <- up-uw # 老王與幸福最高境界的差距
dc <- up-uc # 老張與幸福最高境界的差距
(sqrt(dw%*%dw)) # 老王與幸福最高境界的差距之長度
## [,1]
## [1,] 0.3331666
(sqrt(dc%*%dc)) # 老張與幸福最高境界的差距之長度
## [,1]
## [1,] 0.3541186
## Page-04~05
A <- matrix(c(20,60,40,12,30,20), nrow=3)
# 亦可使用 A <- matrix(c(20,12,60,30,40,20), nrow=3, byrow=T)
B <- matrix(c(25,55,30,18,40,30), nrow=3)
C <- matrix(c(25,50,45,20,35,30), nrow=3)
(M <- A+B+C)
## [,1] [,2]
## [1,] 70 50
## [2,] 165 105
## [3,] 115 80
# 計算矩陣相加結果並印出
D <- matrix(c(2200,1800,1000,850,1500,900),nrow=2)
D
## [,1] [,2] [,3]
## [1,] 2200 1000 1500
## [2,] 1800 850 900
D[1,]%*%M[,1]
## [,1]
## [1,] 491500
# 將矩陣D的第一列與矩陣M的第一行做內積(視為兩向量)
D[2,]%*%M[,2]
## [,1]
## [1,] 251250
# 將矩陣D的第二列與矩陣M的第二行做內積(視為兩向量)
# 以下是將矩陣中的向量取出，再轉換為矩陣的做法
D1 <- as.matrix(D[1,])
# 將矩陣D的第一列視為矩陣D1
M1 <- as.matrix(M[,1])
# 將矩陣M的第一行視為矩陣M1
dim(D1)
## [1] 3 1
# 注意D1為3x1矩陣
class(D1)
## [1] "matrix" "array"
# 確認D1之類型
class(M1)
## [1] "matrix" "array"
# 確認M1之類型
t(D1)%*%M1
## [,1]
## [1,] 491500
# 將D1轉置為1x3矩陣後再與矩陣M1(3x1)相乘
## Page-06~07
A1 <- matrix(c(1,3,2,4),nrow=2)
A2 <- matrix(c(3,1,4,2),nrow=2)
(A1%*%A2)
## [,1] [,2]
## [1,] 5 8
## [2,] 13 20
B1 <- matrix(c(2,5,-1,3,4,0,-1,0,2),nrow=3)
B2 <- matrix(c(0,2,-1,1,0,1),nrow=3)
(B1%*%B2)
## [,1] [,2]
## [1,] 7 1
## [2,] 8 5
## [3,] -2 1
C1 <- matrix(c(2,1,0,-1),nrow=2)
C2 <- matrix(c(3,-2,0,2),nrow=2)
C3 <- matrix(c(1,0,2,3),nrow=2)
(C1%*%C2%*%C3)
## [,1] [,2]
## [1,] 6 12
## [2,] 5 4
D1 <- matrix(c(0,2,1,2),nrow=2)
D2 <- matrix(c(0,1,3,2,-1,0),nrow=2)
D3 <- matrix(c(0,2,1,1,-1,0),nrow=3)
(D1%*%D2%*%D3)
## [,1] [,2]
## [1,] 4 -1
## [2,] 18 -8
MA <- matrix(c(2,1,0,-1),nrow=2)
MB <- matrix(c(3,2,0,-2),nrow=2)
MC <- matrix(c(1,0,2,3),nrow=2)
(MA+2*MB-3*MC)
## [,1] [,2]
## [1,] 5 -6
## [2,] 5 -14
K1 <- matrix(c(-1,0,0,1),nrow=2)
K2 <- matrix(c(0,1,1,0),nrow=2)
K3 <- matrix(c(0.5,0,0,0.5),nrow=2)
MM <- matrix(c(1,3,2,5,2.5,4.5),nrow=2)
MM
## [,1] [,2] [,3]
## [1,] 1 2 2.5
## [2,] 3 5 4.5
K1
## [,1] [,2]
## [1,] -1 0
## [2,] 0 1
K1%*%MM
## [,1] [,2] [,3]
## [1,] -1 -2 -2.5
## [2,] 3 5 4.5
K2
## [,1] [,2]
## [1,] 0 1
## [2,] 1 0
K2%*%MM
## [,1] [,2] [,3]
## [1,] 3 5 4.5
## [2,] 1 2 2.5
K3
## [,1] [,2]
## [1,] 0.5 0.0
## [2,] 0.0 0.5
K3%*%MM
## [,1] [,2] [,3]
## [1,] 0.5 1.0 1.25
## [2,] 1.5 2.5 2.25
(K21 <- K2%*%K1)
## [,1] [,2]
## [1,] 0 1
## [2,] -1 0
(K21%*%MM)
## [,1] [,2] [,3]
## [1,] 3 5 4.5
## [2,] -1 -2 -2.5
# 合併以下兩步驟
#(1)以Y軸做對稱軸做轉換 #(2)以Y=X做對稱軸做轉換
(K12 <- K1%*%K2)
## [,1] [,2]
## [1,] 0 -1
## [2,] 1 0
(K12%*%MM)
## [,1] [,2] [,3]
## [1,] -3 -5 -4.5
## [2,] 1 2 2.5
# 合併以下兩步驟
# (1)以Y=X做對稱軸做轉換
# (2)以Y軸做對稱軸做轉換
(K1%*%K1)
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
# 以Y軸做對稱軸做兩次轉換
(K2%*%K2)
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
# 以Y=X做對稱軸做兩次轉換
## Page-08~09
A <- matrix(c(8,4,4,7),nrow=2)
solve(A)
## [,1] [,2]
## [1,] 0.175 -0.1
## [2,] -0.100 0.2
B <- matrix(c(0,3,2,1),nrow=2)
solve(B)
## [,1] [,2]
## [1,] -0.1666667 0.3333333
## [2,] 0.5000000 0.0000000
polyroot(c(-8,8,1))
## [1] 0.8989795+0i -8.8989795+0i
# 解 k^2+8k-8 = 0，注意多項式係數升冪順序
C <- matrix(c(0,-1,1,0),nrow=2)
solve(C)
## [,1] [,2]
## [1,] 0 -1
## [2,] 1 0
(D <- matrix(c(-2,0,0,-2),nrow=2))
## [,1] [,2]
## [1,] -2 0
## [2,] 0 -2
(E <- solve(D))
## [,1] [,2]
## [1,] -0.5 0.0
## [2,] 0.0 -0.5
M <- matrix(c(2,8,4,7,4,5),nrow=2)
# 三個點座標構成的矩陣M
(MN <- D%*%M)
## [,1] [,2] [,3]
## [1,] -4 -8 -8
## [2,] -16 -14 -10
# 三個點透過矩陣D轉換得新座標矩陣MN
(E%*%MN)
## [,1] [,2] [,3]
## [1,] 2 4 4
## [2,] 8 7 5
# 新座標再透過D的反矩陣E轉換回原來的座標
## Page-11
A <- matrix(c(1,2,1,4),nrow=2)
AX <- matrix(c(10,28,1,4),nrow=2)
AY <- matrix(c(1,2,10,28),nrow=2)
cat("X=",det(AX)/det(A), "\n")
## X= 6
# 計算 x 之值並印出，換行
cat("Y=",det(AY)/det(A), "\n")
## Y= 4
# 計算 y 之值並印出，換行
B <- matrix(c(4, 6, 5, -3),nrow=2)
BX <- matrix(c(7, 21, 5, -3),nrow=2)
BY <- matrix(c(4,6, 7, 21),nrow=2)
cat("X=",det(BX)/det(B), "\n")
## X= 3
# 計算 x 之值並印出，換行
cat("Y=",det(BY)/det(B), "\n")
## Y= -1
# 計算 y 之值並印出，換行
C <- matrix(c(0.8,0.3,0.6,0.4),nrow=2)
CX <- matrix(c(4,3,0.6,0.4),nrow=2)
CY <- matrix(c(0.8,0.3,4,3),nrow=2)
cat("X=",det(CX)/det(C), "\n")
## X= -1.428571
# 計算 x 之值並印出，換行
cat("Y=",det(CY)/det(C), "\n")
## Y= 8.571429
# 計算 y 之值並印出，換行
D <- matrix(c(3/5,1/4,1/2,4/3),nrow=2)
DX <- matrix(c(2,4,1/2,4/3),nrow=2)
DY <- matrix(c(3/5,1/4,2,4),nrow=2)
cat("X=",det(DX)/det(D), "\n")
## X= 0.9876543
# 計算 x 之值並印出，換行
cat("Y=",det(DY)/det(D), "\n")
## Y= 2.814815
# 計算 y 之值並印出，換行
D <- matrix(c(1,2,2,1,4,0,1,6,4),nrow=3)
DX <- matrix(c(10,38,14,1,4,0,1,6,4),nrow=3)
DY <- matrix(c(1,2,2,10,38,14,1,6,4),nrow=3)
DZ <- matrix(c(1,2,2,1,4,0,10,38,14),nrow=3)
cat("X=",det(DX)/det(D), "\n")
## X= 3
# 計算 x 之值並印出，換行
cat("Y=",det(DY)/det(D), "\n")
## Y= 5
# 計算 y 之值並印出，換行
cat("Y=",det(DZ)/det(D), "\n")
## Y= 2
# 計算 z 之值並印出，換行
## Page-13
M <- matrix(c(3,4,2,8,1,0,5,3,2),nrow=3)
MX <- matrix(c(60,46,24,8,1,0,5,3,2),nrow=3)
MY <- matrix(c(3,4,2,60,46,24,5,3,2),nrow=3)
MZ <- matrix(c(3,4,2,8,1,0,60,46,24),nrow=3)
cat("X=",det(MX)/det(M), "\n")
## X= 8
# 計算 x 之值並印出，換行
cat("Y=",det(MY)/det(M), "\n")
## Y= 2
# 計算 y 之值並印出，換行
cat("Y=",det(MZ)/det(M), "\n")
## Y= 4
# 計算 z 之值並印出，換行
b <- matrix(c(60,46,24),nrow=3)
MS <- solve(M)
# 計算 M 的反矩陣
cat("x y z=", MS%*%b, "\n")
## x y z= 8 2 4
# 解出(x,y,z)向量，並印出
## Page-14~15
A <- matrix(c(1,2,2,1,4,0,1,6,4),nrow=3)
b <- matrix(c(10,38,14),nrow=3)
(Ab <-cbind(A,b))
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 10
## [2,] 2 4 6 38
## [3,] 2 0 4 14
# 利用cbind指令將矩陣A與b合併(並列)
Ab[2,] <- Ab[1,]*(-2)+Ab[2,]
# 將合併矩陣之第一列乘以-2加到第二列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 10
## [2,] 0 2 4 18
## [3,] 2 0 4 14
Ab[3,] <- Ab[1,]*(-2)+Ab[3,]
# 將合併矩陣之第一列乘以-2加到第三列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 10
## [2,] 0 2 4 18
## [3,] 0 -2 2 -6
Ab[2,] <- Ab[2,]*(1/2)
# 將合併矩陣之第二列乘以1/2
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 10
## [2,] 0 1 2 9
## [3,] 0 -2 2 -6
Ab[1,] <- Ab[2,]*(-1)+Ab[1,]
# 將合併矩陣之第二列乘以-1加到第一列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 -1 1
## [2,] 0 1 2 9
## [3,] 0 -2 2 -6
Ab[3,] <- Ab[2,]*2+Ab[3,]
# 將合併矩陣之第二列乘以2加到第三列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 -1 1
## [2,] 0 1 2 9
## [3,] 0 0 6 12
Ab[3,] <- Ab[3,]*(1/6)
# 將合併矩陣之第三列乘以1/6
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 -1 1
## [2,] 0 1 2 9
## [3,] 0 0 1 2
Ab[1,] <- Ab[3,]*1+Ab[1,]
# 將合併矩陣之第三列乘以1加到第一列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 3
## [2,] 0 1 2 9
## [3,] 0 0 1 2
Ab[2,] <- Ab[3,]*(-2)+Ab[2,]
# 將合併矩陣之第三列乘以-2 加到第二列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 3
## [2,] 0 1 0 5
## [3,] 0 0 1 2
A <- matrix(c(1,2,2,1,4,0,1,6,4),nrow=3)
(AI <-cbind(A,diag(3)))
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 1 1 1 0 0
## [2,] 2 4 6 0 1 0
## [3,] 2 0 4 0 0 1
# 利用cbind指令將矩陣A與I3合併(並列)
AI[2,] <- AI[1,]*(-2)+AI[2,]
# 將合併矩陣之第一列乘以-2加到第二列
AI
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 1 1 1 0 0
## [2,] 0 2 4 -2 1 0
## [3,] 2 0 4 0 0 1
AI[3,] <- AI[1,]*(-2)+AI[3,]
# 將合併矩陣之第一列乘以-2加到第三列
AI
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 1 1 1 0 0
## [2,] 0 2 4 -2 1 0
## [3,] 0 -2 2 -2 0 1
AI[2,] <- AI[2,]*(1/2)
# 將合併矩陣之第二列乘以1/2
AI
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 1 1 1 0.0 0
## [2,] 0 1 2 -1 0.5 0
## [3,] 0 -2 2 -2 0.0 1
AI[1,] <- AI[2,]*(-1)+AI[1,]
# 將合併矩陣之第二列乘以-1加到第一列
AI
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 -1 2 -0.5 0
## [2,] 0 1 2 -1 0.5 0
## [3,] 0 -2 2 -2 0.0 1
AI[3,] <- AI[2,]*2+AI[3,]
# 將合併矩陣之第二列乘以2加到第三列
AI
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 -1 2 -0.5 0
## [2,] 0 1 2 -1 0.5 0
## [3,] 0 0 6 -4 1.0 1
AI[3,] <- AI[3,]*(1/6)
# 將合併矩陣之第三列乘以1/6
AI
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 -1 2.0000000 -0.5000000 0.0000000
## [2,] 0 1 2 -1.0000000 0.5000000 0.0000000
## [3,] 0 0 1 -0.6666667 0.1666667 0.1666667
AI[1,] <- AI[3,]*1+AI[1,]
# 將合併矩陣之第三列乘以1加到第一列
AI
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 1.3333333 -0.3333333 0.1666667
## [2,] 0 1 2 -1.0000000 0.5000000 0.0000000
## [3,] 0 0 1 -0.6666667 0.1666667 0.1666667
AI[2,] <- AI[3,]*(-2)+AI[2,]
# 將合併矩陣之第三列乘以-2 加到第二列
AI
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 1.3333333 -0.3333333 0.1666667
## [2,] 0 1 0 0.3333333 0.1666667 -0.3333333
## [3,] 0 0 1 -0.6666667 0.1666667 0.1666667
# 亦可直接利用 A 的逆矩陣 A^{-1} 求解
(AIN <- AI[,4:6])
## [,1] [,2] [,3]
## [1,] 1.3333333 -0.3333333 0.1666667
## [2,] 0.3333333 0.1666667 -0.3333333
## [3,] -0.6666667 0.1666667 0.1666667
# 合併矩陣之後三欄即為欲求之反矩陣
cat("x y z=", AIN%*%b, "\n")
## x y z= 3 5 2
# 印出(x,y,z)之解
## Page-16
A <- matrix(c(1,2,4,5,4,-2,-2,3,1),nrow=3)
b <- matrix(c(21,17,9),nrow=3)
(Ab <-cbind(A,b))
## [,1] [,2] [,3] [,4]
## [1,] 1 5 -2 21
## [2,] 2 4 3 17
## [3,] 4 -2 1 9
# 利用cbind指令將矩陣A與b合併(並列)
Ab[2,] <- Ab[1,]*(-2)+Ab[2,]
# 將合併矩陣之第一列乘以-2加到第二列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 5 -2 21
## [2,] 0 -6 7 -25
## [3,] 4 -2 1 9
Ab[3,] <- Ab[1,]*(-4)+Ab[3,]
# 將合併矩陣之第一列乘以-4加到第三列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 5 -2 21
## [2,] 0 -6 7 -25
## [3,] 0 -22 9 -75
Ab[2,] <- Ab[2,]*(-1/6)
# 將合併矩陣之第二列乘以 -1/6
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 5 -2.000000 21.000000
## [2,] 0 1 -1.166667 4.166667
## [3,] 0 -22 9.000000 -75.000000
Ab[1,] <- Ab[2,]*(-5)+Ab[1,]
# 將合併矩陣之第二列乘以-5加到第一列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 3.833333 0.1666667
## [2,] 0 1 -1.166667 4.1666667
## [3,] 0 -22 9.000000 -75.0000000
Ab[3,] <- Ab[2,]*(22)+Ab[3,]
# 將合併矩陣之第二列乘以22加到第三列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 3.833333 0.1666667
## [2,] 0 1 -1.166667 4.1666667
## [3,] 0 0 -16.666667 16.6666667
Ab[3,] <- Ab[3,]*(-1)/(16+2/3)
# 將合併矩陣之第三列乘以-1/(16+2/3)
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 3.833333 0.1666667
## [2,] 0 1 -1.166667 4.1666667
## [3,] 0 0 1.000000 -1.0000000
Ab[1,] <- Ab[3,]*(-(3+5/6))+Ab[1,]
# 將合併矩陣之第三列乘以(-(3+5/6))加到第一列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 -4.440892e-16 4.000000
## [2,] 0 1 -1.166667e+00 4.166667
## [3,] 0 0 1.000000e+00 -1.000000
Ab[2,] <- Ab[3,]*(1+1/6)+Ab[2,]
# 將合併矩陣之第三列乘以(1+1/6)加到第二列
Ab
## [,1] [,2] [,3] [,4]
## [1,] 1 0 -4.440892e-16 4
## [2,] 0 1 0.000000e+00 3
## [3,] 0 0 1.000000e+00 -1
A <- matrix(c(1,2,1,-1,3,4,2,1,-1),nrow=3)
b1 <- matrix(c(4,6,2),nrow=3)
b2 <- matrix(c(4,6,5),nrow=3)
(Ab12 <-cbind(A,b1,b2))
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -1 2 4 4
## [2,] 2 3 1 6 6
## [3,] 1 4 -1 2 5
# 利用cbind指令將矩陣A與b1、b2合併(並列)
Ab12[2,] <- Ab12[1,]*(-2)+Ab12[2,]
# 將合併矩陣之第一列乘以-2加到第二列
Ab12
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -1 2 4 4
## [2,] 0 5 -3 -2 -2
## [3,] 1 4 -1 2 5
Ab12[3,] <- Ab12[1,]*(-1)+Ab12[3,]
# 將合併矩陣之第一列乘以-1加到第三列
Ab12
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -1 2 4 4
## [2,] 0 5 -3 -2 -2
## [3,] 0 5 -3 -2 1
Ab12[2,] <- Ab12[2,]*(1/5)
# 將合併矩陣之第二列乘以1/5
Ab12
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -1 2.0 4.0 4.0
## [2,] 0 1 -0.6 -0.4 -0.4
## [3,] 0 5 -3.0 -2.0 1.0
Ab12[1,] <- Ab12[2,]*(1)+Ab12[1,]
# 將合併矩陣之第二列乘以1加到第一列
Ab12
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 1.4 3.6 3.6
## [2,] 0 1 -0.6 -0.4 -0.4
## [3,] 0 5 -3.0 -2.0 1.0
Ab12[3,] <- Ab12[2,]*(-5)+Ab12[3,]
# 將合併矩陣之第二列乘以-5加到第三列
Ab12
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 1.400000e+00 3.6 3.6
## [2,] 0 1 -6.000000e-01 -0.4 -0.4
## [3,] 0 0 4.440892e-16 0.0 3.0
# 對於b1，得 x1+1.4*x3=3.6 , x2-0.6*x2=-0.4 , 0*x1+0*x2+0*x3=0 => 無窮多解
# 對於b2，得 x1+1.4*x3=3.6 , x2-0.6*x2=-0.4 , 0*x1+0*x2+0*x3=3.0 => 無解
A <- matrix(c(1,2,4,0,1,0,-1,1,0,1,3,2,-1,3,1,2),nrow=4)
(AI <-cbind(A,diag(4)))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 1 0 -1 1 0 0 0
## [2,] 2 0 1 3 0 1 0 0
## [3,] 4 -1 3 1 0 0 1 0
## [4,] 0 1 2 2 0 0 0 1
# 利用cbind指令將矩陣A與I4合併(並列)
AI[2,] <- AI[1,]*(-2)+AI[2,]
# 將合併矩陣之第一列乘以-2加到第二列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 1 0 -1 1 0 0 0
## [2,] 0 -2 1 5 -2 1 0 0
## [3,] 4 -1 3 1 0 0 1 0
## [4,] 0 1 2 2 0 0 0 1
AI[3,] <- AI[1,]*(-4)+AI[3,]
# 將合併矩陣之第一列乘以-4加到第三列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 1 0 -1 1 0 0 0
## [2,] 0 -2 1 5 -2 1 0 0
## [3,] 0 -5 3 5 -4 0 1 0
## [4,] 0 1 2 2 0 0 0 1
AI[2,] <- AI[2,]*(-1/2)
# 將合併矩陣之第二列乘以-1/2
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 1 0.0 -1.0 1 0.0 0 0
## [2,] 0 1 -0.5 -2.5 1 -0.5 0 0
## [3,] 0 -5 3.0 5.0 -4 0.0 1 0
## [4,] 0 1 2.0 2.0 0 0.0 0 1
AI[1,] <- AI[2,]*(-1)+AI[1,]
# 將合併矩陣之第二列乘以-1加到第一列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0.5 1.5 0 0.5 0 0
## [2,] 0 1 -0.5 -2.5 1 -0.5 0 0
## [3,] 0 -5 3.0 5.0 -4 0.0 1 0
## [4,] 0 1 2.0 2.0 0 0.0 0 1
AI[3,] <- AI[2,]*(5)+AI[3,]
# 將合併矩陣之第二列乘以5加到第三列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0.5 1.5 0 0.5 0 0
## [2,] 0 1 -0.5 -2.5 1 -0.5 0 0
## [3,] 0 0 0.5 -7.5 1 -2.5 1 0
## [4,] 0 1 2.0 2.0 0 0.0 0 1
AI[4,] <- AI[2,]*(-1)+AI[4,]
# 將合併矩陣之第二列乘以-1加到第四列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0.5 1.5 0 0.5 0 0
## [2,] 0 1 -0.5 -2.5 1 -0.5 0 0
## [3,] 0 0 0.5 -7.5 1 -2.5 1 0
## [4,] 0 0 2.5 4.5 -1 0.5 0 1
AI[3,] <- AI[3,]*(2)
# 將合併矩陣之第三列乘以2
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0.5 1.5 0 0.5 0 0
## [2,] 0 1 -0.5 -2.5 1 -0.5 0 0
## [3,] 0 0 1.0 -15.0 2 -5.0 2 0
## [4,] 0 0 2.5 4.5 -1 0.5 0 1
AI[1,] <- AI[3,]*(-0.5)+AI[1,]
# 將合併矩陣之第三列乘以-0.5加到第一列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0.0 9.0 -1 3.0 -1 0
## [2,] 0 1 -0.5 -2.5 1 -0.5 0 0
## [3,] 0 0 1.0 -15.0 2 -5.0 2 0
## [4,] 0 0 2.5 4.5 -1 0.5 0 1
AI[2,] <- AI[3,]*(0.5)+AI[2,]
# 將合併矩陣之第三列乘以0.5加到第二列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0.0 9.0 -1 3.0 -1 0
## [2,] 0 1 0.0 -10.0 2 -3.0 1 0
## [3,] 0 0 1.0 -15.0 2 -5.0 2 0
## [4,] 0 0 2.5 4.5 -1 0.5 0 1
AI[4,] <- AI[3,]*(-2.5)+AI[4,]
# 將合併矩陣之第三列乘以-2.5加到第四列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0 9 -1 3 -1 0
## [2,] 0 1 0 -10 2 -3 1 0
## [3,] 0 0 1 -15 2 -5 2 0
## [4,] 0 0 0 42 -6 13 -5 1
AI[4,] <- AI[4,]*(1/42)
# 將合併矩陣之第四列乘以1/42
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0 9 -1.0000000 3.0000000 -1.0000000 0.00000000
## [2,] 0 1 0 -10 2.0000000 -3.0000000 1.0000000 0.00000000
## [3,] 0 0 1 -15 2.0000000 -5.0000000 2.0000000 0.00000000
## [4,] 0 0 0 1 -0.1428571 0.3095238 -0.1190476 0.02380952
AI[1,] <- AI[4,]*(-9)+AI[1,]
# 將合併矩陣之第四列乘以-9加到第一列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0 0 0.2857143 0.2142857 0.07142857 -0.21428571
## [2,] 0 1 0 -10 2.0000000 -3.0000000 1.00000000 0.00000000
## [3,] 0 0 1 -15 2.0000000 -5.0000000 2.00000000 0.00000000
## [4,] 0 0 0 1 -0.1428571 0.3095238 -0.11904762 0.02380952
AI[2,] <- AI[4,]*(10)+AI[2,]
# 將合併矩陣之第四列乘以10加到第二列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0 0 0.2857143 0.2142857 0.07142857 -0.21428571
## [2,] 0 1 0 0 0.5714286 0.0952381 -0.19047619 0.23809524
## [3,] 0 0 1 -15 2.0000000 -5.0000000 2.00000000 0.00000000
## [4,] 0 0 0 1 -0.1428571 0.3095238 -0.11904762 0.02380952
AI[3,] <- AI[4,]*(15)+AI[3,]
# 將合併矩陣之第四列乘以15加到第三列
AI
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0 0 0.2857143 0.2142857 0.07142857 -0.21428571
## [2,] 0 1 0 0 0.5714286 0.0952381 -0.19047619 0.23809524
## [3,] 0 0 1 0 -0.1428571 -0.3571429 0.21428571 0.35714286
## [4,] 0 0 0 1 -0.1428571 0.3095238 -0.11904762 0.02380952
## Page-18
A <- matrix(c(-1, 0, 0, 1),nrow=2)
eigen(A)
## eigen() decomposition
## $values
## [1] 1 -1
##
## $vectors
## [,1] [,2]
## [1,] 0 -1
## [2,] 1 0
A <- matrix(c(0.5, 0, 0, 0.5),nrow=2)
eigen(A)
## eigen() decomposition
## $values
## [1] 0.5 0.5
##
## $vectors
## [,1] [,2]
## [1,] 0 -1
## [2,] 1 0
A <- matrix(c(-4, 3, -6, 5),nrow=2)
eigen(A)
## eigen() decomposition
## $values
## [1] 2 -1
##
## $vectors
## [,1] [,2]
## [1,] 0.7071068 -0.8944272
## [2,] -0.7071068 0.4472136
A <- matrix(c(5, 4, 2, 4, 5, 2, 2, 2, 2),nrow=3)
eigen(A)
## eigen() decomposition
## $values
## [1] 10 1 1
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.6666667 0.0000000 0.7453560
## [2,] -0.6666667 -0.4472136 -0.5962848
## [3,] -0.3333333 0.8944272 -0.2981424
## Page-23~25
(A <- matrix(c(1,0,0,0,0,0,0,4,0,3,0,0,0,0,0,0,2,0,0,0),nrow=4))
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 2
## [2,] 0 0 3 0 0
## [3,] 0 0 0 0 0
## [4,] 0 4 0 0 0
svd(A,nu=4,nv=5) #求矩陣A之SVD分解，設定U為4x4，V為5x5
## $d
## [1] 4.000000 3.000000 2.236068 0.000000
##
## $u
## [,1] [,2] [,3] [,4]
## [1,] 0 0 1 0
## [2,] 0 1 0 0
## [3,] 0 0 0 1
## [4,] -1 0 0 0
##
## $v
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0.4472136 0 -0.8944272
## [2,] -1 0 0.0000000 0 0.0000000
## [3,] 0 1 0.0000000 0 0.0000000
## [4,] 0 0 0.0000000 1 0.0000000
## [5,] 0 0 0.8944272 0 0.4472136
(A <- matrix(c(1,0,0,3,0,0,4,0,0,1,0,2),nrow=4))
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 0 1
## [3,] 0 4 0
## [4,] 3 0 2
svd(A,nu=4,nv=3) #求矩陣A之SVD分解，設定U為4x4，V為3x3
## $d
## [1] 4.000000 3.741657 1.000000
##
## $u
## [,1] [,2] [,3] [,4]
## [1,] 0 -0.2223748 5.547002e-01 -0.8017837
## [2,] 0 -0.1482499 -8.320503e-01 -0.5345225
## [3,] 1 0.0000000 0.000000e+00 0.0000000
## [4,] 0 -0.9636241 -1.110223e-16 0.2672612
##
## $v
## [,1] [,2] [,3]
## [1,] 0 -0.8320503 0.5547002
## [2,] 1 0.0000000 0.0000000
## [3,] 0 -0.5547002 -0.8320503
## 自行構建一個10x10矩陣(影片例)
(A <- matrix(c(rep(0,20),rep(c(0,1,1,1,1,1,1,1,1,0),2),rep(c(0,1,1,0,0,0,0,1,1,0),2),rep(c(0,1,1,1,1,1,1,1,1,0),2),rep(0,20)),nrow=10))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 1 1 1 1 1 1 0 0
## [3,] 0 0 1 1 1 1 1 1 0 0
## [4,] 0 0 1 1 0 0 1 1 0 0
## [5,] 0 0 1 1 0 0 1 1 0 0
## [6,] 0 0 1 1 0 0 1 1 0 0
## [7,] 0 0 1 1 0 0 1 1 0 0
## [8,] 0 0 1 1 1 1 1 1 0 0
## [9,] 0 0 1 1 1 1 1 1 0 0
## [10,] 0 0 0 0 0 0 0 0 0 0
image(A,col=topo.colors(256))
# 用image指令畫圖
s <- svd(A,nu=10,nv=10)
# 將svd分解的結果存入s
(D <- diag(s$d))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 6.040896 0.000000 0.000000e+00 0.000000e+00 0.0000e+00 0.000000e+00 0
## [2,] 0.000000 1.872853 0.000000e+00 0.000000e+00 0.0000e+00 0.000000e+00 0
## [3,] 0.000000 0.000000 7.918465e-16 0.000000e+00 0.0000e+00 0.000000e+00 0
## [4,] 0.000000 0.000000 0.000000e+00 1.275686e-16 0.0000e+00 0.000000e+00 0
## [5,] 0.000000 0.000000 0.000000e+00 0.000000e+00 3.9807e-33 0.000000e+00 0
## [6,] 0.000000 0.000000 0.000000e+00 0.000000e+00 0.0000e+00 1.670685e-49 0
## [7,] 0.000000 0.000000 0.000000e+00 0.000000e+00 0.0000e+00 0.000000e+00 0
## [8,] 0.000000 0.000000 0.000000e+00 0.000000e+00 0.0000e+00 0.000000e+00 0
## [9,] 0.000000 0.000000 0.000000e+00 0.000000e+00 0.0000e+00 0.000000e+00 0
## [10,] 0.000000 0.000000 0.000000e+00 0.000000e+00 0.0000e+00 0.000000e+00 0
## [,8] [,9] [,10]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 0
## [5,] 0 0 0
## [6,] 0 0 0
## [7,] 0 0 0
## [8,] 0 0 0
## [9,] 0 0 0
## [10,] 0 0 0
# 取出奇異值做成矩陣
image(D,col=topo.colors(256))
# 用image指令畫圖
AA <- t(A)%*%A
# 計算A的轉置矩陣乘以A
g <- eigen(AA) # 將特徵值與特徵向量的計算結果存入g
(V <- diag(g$values)) # 取出特徵值做成矩陣
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 36.49242 0.000000 0.000000e+00 0.00000e+00 0.000000e+00 0 0 0
## [2,] 0.00000 3.507577 0.000000e+00 0.00000e+00 0.000000e+00 0 0 0
## [3,] 0.00000 0.000000 3.647876e-16 0.00000e+00 0.000000e+00 0 0 0
## [4,] 0.00000 0.000000 0.000000e+00 1.97373e-16 0.000000e+00 0 0 0
## [5,] 0.00000 0.000000 0.000000e+00 0.00000e+00 4.930381e-32 0 0 0
## [6,] 0.00000 0.000000 0.000000e+00 0.00000e+00 0.000000e+00 0 0 0
## [7,] 0.00000 0.000000 0.000000e+00 0.00000e+00 0.000000e+00 0 0 0
## [8,] 0.00000 0.000000 0.000000e+00 0.00000e+00 0.000000e+00 0 0 0
## [9,] 0.00000 0.000000 0.000000e+00 0.00000e+00 0.000000e+00 0 0 0
## [10,] 0.00000 0.000000 0.000000e+00 0.00000e+00 0.000000e+00 0 0 0
## [,9] [,10]
## [1,] 0 0.000000e+00
## [2,] 0 0.000000e+00
## [3,] 0 0.000000e+00
## [4,] 0 0.000000e+00
## [5,] 0 0.000000e+00
## [6,] 0 0.000000e+00
## [7,] 0 0.000000e+00
## [8,] 0 0.000000e+00
## [9,] 0 0.000000e+00
## [10,] 0 -4.749287e-16
image(V,col=topo.colors(256)) # 用image指令畫圖
## 擷取矩陣的一部分
(RA <- A[3:8,3:8])
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 1 1 1 1 1
## [2,] 1 1 0 0 1 1
## [3,] 1 1 0 0 1 1
## [4,] 1 1 0 0 1 1
## [5,] 1 1 0 0 1 1
## [6,] 1 1 1 1 1 1
image(RA,col=topo.colors(256))
# 用image指令畫圖
s <- svd(RA,nu=6,nv=6)
# 將svd分解的結果存入s
(D <- diag(s$d))
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 5.048676 0.000000 0.000000e+00 0.000000e+00 0.000000e+00 0
## [2,] 0.000000 1.584574 0.000000e+00 0.000000e+00 0.000000e+00 0
## [3,] 0.000000 0.000000 6.613243e-16 0.000000e+00 0.000000e+00 0
## [4,] 0.000000 0.000000 0.000000e+00 1.695558e-17 0.000000e+00 0
## [5,] 0.000000 0.000000 0.000000e+00 0.000000e+00 5.908398e-34 0
## [6,] 0.000000 0.000000 0.000000e+00 0.000000e+00 0.000000e+00 0
# 取出奇異值做成矩陣
image(D,col=topo.colors(256))
# 用image指令畫圖
## 透過函數繪製影像圖(修改自 R 範例)
x <- y <- seq(-4*pi, 4*pi, len = 500)
r <- sqrt(outer(x^2, y^2, "+"))
image(z = z <- cos(r^2)*exp(-r/6), col = gray((0:256)/256))
image(z, axes = FALSE, main = "Math can be beautiful ...",
xlab = expression(cos(r^2) * e^{-r/6}))
image(z, axes = FALSE, main = "Math can be beautiful ...",
xlab = expression(cos(r^2) * e^{-r/6}),col=cm.colors(16))
## Page37
## 模擬醉漢走路，路寬左右各十米處有懸崖
x<-0 ; count<-0
repeat{x<- x + sample(c(-1,1),1) ; print(x); count<-count+1; if(x>=10 || x<=-10) {print(count);break}}
## [1] 1
## [1] 0
## [1] 1
## [1] 0
## [1] -1
## [1] 0
## [1] -1
## [1] -2
## [1] -1
## [1] -2
## [1] -3
## [1] -4
## [1] -5
## [1] -6
## [1] -5
## [1] -6
## [1] -5
## [1] -6
## [1] -7
## [1] -6
## [1] -5
## [1] -4
## [1] -5
## [1] -6
## [1] -7
## [1] -8
## [1] -9
## [1] -8
## [1] -7
## [1] -6
## [1] -5
## [1] -6
## [1] -7
## [1] -8
## [1] -9
## [1] -10
## [1] 36
## 另一種寫法，紀錄過程
x<-c() ; count<-0
repeat{x<- rbind(x, sample(c(-1,1),1)) ; count<-count+1; if(sum(x)>=10 || sum(x)<=-10) {print(x); print(count); break}}
## [,1]
## [1,] 1
## [2,] 1
## [3,] -1
## [4,] 1
## [5,] -1
## [6,] -1
## [7,] -1
## [8,] 1
## [9,] 1
## [10,] 1
## [11,] -1
## [12,] 1
## [13,] -1
## [14,] -1
## [15,] -1
## [16,] -1
## [17,] -1
## [18,] -1
## [19,] -1
## [20,] 1
## [21,] 1
## [22,] 1
## [23,] -1
## [24,] -1
## [25,] 1
## [26,] 1
## [27,] -1
## [28,] -1
## [29,] -1
## [30,] 1
## [31,] 1
## [32,] -1
## [33,] 1
## [34,] 1
## [35,] 1
## [36,] 1
## [37,] -1
## [38,] -1
## [39,] 1
## [40,] -1
## [41,] -1
## [42,] -1
## [43,] 1
## [44,] 1
## [45,] 1
## [46,] 1
## [47,] 1
## [48,] 1
## [49,] 1
## [50,] 1
## [51,] -1
## [52,] 1
## [53,] -1
## [54,] 1
## [55,] -1
## [56,] -1
## [57,] -1
## [58,] -1
## [59,] 1
## [60,] 1
## [61,] -1
## [62,] -1
## [63,] 1
## [64,] -1
## [65,] 1
## [66,] -1
## [67,] -1
## [68,] 1
## [69,] -1
## [70,] 1
## [71,] -1
## [72,] 1
## [73,] 1
## [74,] 1
## [75,] 1
## [76,] -1
## [77,] -1
## [78,] 1
## [79,] 1
## [80,] -1
## [81,] -1
## [82,] -1
## [83,] -1
## [84,] -1
## [85,] 1
## [86,] -1
## [87,] 1
## [88,] -1
## [89,] -1
## [90,] -1
## [91,] 1
## [92,] 1
## [93,] 1
## [94,] 1
## [95,] 1
## [96,] 1
## [97,] 1
## [98,] 1
## [99,] -1
## [100,] 1
## [101,] 1
## [102,] 1
## [103,] -1
## [104,] -1
## [105,] -1
## [106,] 1
## [107,] -1
## [108,] 1
## [109,] -1
## [110,] -1
## [111,] -1
## [112,] 1
## [113,] 1
## [114,] -1
## [115,] 1
## [116,] -1
## [117,] -1
## [118,] 1
## [119,] -1
## [120,] 1
## [121,] -1
## [122,] -1
## [123,] -1
## [124,] 1
## [125,] -1
## [126,] 1
## [127,] -1
## [128,] 1
## [129,] 1
## [130,] 1
## [131,] 1
## [132,] -1
## [133,] 1
## [134,] -1
## [135,] 1
## [136,] 1
## [137,] 1
## [138,] 1
## [139,] -1
## [140,] -1
## [141,] 1
## [142,] 1
## [143,] -1
## [144,] 1
## [145,] 1
## [146,] 1
## [147,] -1
## [148,] -1
## [149,] -1
## [150,] 1
## [151,] 1
## [152,] 1
## [153,] -1
## [154,] -1
## [155,] 1
## [156,] 1
## [157,] -1
## [158,] 1
## [159,] -1
## [160,] -1
## [161,] -1
## [162,] 1
## [163,] -1
## [164,] -1
## [165,] -1
## [166,] -1
## [167,] -1
## [168,] -1
## [169,] -1
## [170,] -1
## [171,] 1
## [172,] -1
## [173,] 1
## [174,] -1
## [175,] 1
## [176,] -1
## [177,] -1
## [178,] 1
## [179,] 1
## [180,] -1
## [181,] 1
## [182,] -1
## [183,] 1
## [184,] -1
## [185,] 1
## [186,] 1
## [187,] -1
## [188,] -1
## [189,] 1
## [190,] -1
## [191,] 1
## [192,] 1
## [193,] -1
## [194,] -1
## [195,] 1
## [196,] -1
## [197,] 1
## [198,] -1
## [199,] 1
## [200,] -1
## [201,] -1
## [202,] -1
## [203,] -1
## [204,] -1
## [205,] -1
## [206,] -1
## [207,] 1
## [208,] -1
## [209,] 1
## [210,] 1
## [211,] -1
## [212,] -1
## [213,] -1
## [214,] 1
## [215,] 1
## [216,] -1
## [217,] -1
## [218,] -1
## [1] 218
## 為什麼你通常會輸光? (模擬100次，賭徒與莊家財力 10:100)
for ( i in 1:100) {
x<-0 ; count<-0
repeat{x<- x + sample(c(-1,1),1) ; count<-count+1; if(x>=100 || x<=-10) {cat("x=",x,"count=",count, "\n"); break}} }
## x= -10 count= 34
## x= -10 count= 36
## x= -10 count= 424
## x= -10 count= 40
## x= -10 count= 214
## x= 100 count= 2760
## x= 100 count= 2396
## x= -10 count= 380
## x= -10 count= 16
## x= -10 count= 130
## x= 100 count= 4600
## x= -10 count= 4440
## x= -10 count= 24
## x= -10 count= 162
## x= 100 count= 3920
## x= 100 count= 980
## x= -10 count= 108
## x= 100 count= 3442
## x= 100 count= 3832
## x= -10 count= 72
## x= -10 count= 136
## x= -10 count= 10
## x= -10 count= 270
## x= -10 count= 160
## x= -10 count= 54
## x= -10 count= 24
## x= -10 count= 146
## x= -10 count= 504
## x= -10 count= 970
## x= 100 count= 9184
## x= -10 count= 244
## x= -10 count= 30
## x= -10 count= 14442
## x= -10 count= 6340
## x= -10 count= 68
## x= -10 count= 44
## x= -10 count= 972
## x= -10 count= 114
## x= -10 count= 406
## x= -10 count= 60
## x= -10 count= 210
## x= -10 count= 28
## x= -10 count= 170
## x= 100 count= 774
## x= -10 count= 836
## x= -10 count= 70
## x= -10 count= 1238
## x= -10 count= 932
## x= -10 count= 278
## x= -10 count= 270
## x= -10 count= 4260
## x= -10 count= 156
## x= -10 count= 248
## x= -10 count= 974
## x= 100 count= 5068
## x= -10 count= 554
## x= -10 count= 2272
## x= -10 count= 104
## x= -10 count= 54
## x= -10 count= 3184
## x= -10 count= 120
## x= 100 count= 3644
## x= -10 count= 56
## x= -10 count= 176
## x= -10 count= 1560
## x= -10 count= 496
## x= -10 count= 26
## x= -10 count= 270
## x= -10 count= 212
## x= -10 count= 30
## x= 100 count= 1396
## x= -10 count= 204
## x= -10 count= 92
## x= -10 count= 22
## x= -10 count= 64
## x= -10 count= 14
## x= -10 count= 732
## x= -10 count= 2944
## x= 100 count= 2818
## x= -10 count= 2968
## x= -10 count= 156
## x= 100 count= 4048
## x= -10 count= 618
## x= -10 count= 294
## x= -10 count= 22
## x= -10 count= 14
## x= -10 count= 1736
## x= -10 count= 100
## x= -10 count= 776
## x= -10 count= 1074
## x= -10 count= 96
## x= -10 count= 4494
## x= -10 count= 794
## x= -10 count= 54
## x= -10 count= 26
## x= -10 count= 94
## x= -10 count= 62
## x= -10 count= 68
## x= -10 count= 50
## x= -10 count= 1328
## 為什麼你通常會輸光? (模擬100次，調高莊家的勝率，賭徒與莊家財力 10:100)
for ( i in 1:100) {
x<-0 ; count<-0
repeat{x<- x + sample(c(-1,1,-1,1,-1,1,-1,1,-1,1,-1),1) ; count<-count+1; if(x>=100 || x<=-10) {cat("x=",x,"count=",count, "\n"); break}} }
## x= -10 count= 414
## x= -10 count= 26
## x= -10 count= 148
## x= -10 count= 80
## x= -10 count= 140
## x= -10 count= 60
## x= -10 count= 196
## x= -10 count= 42
## x= -10 count= 74
## x= -10 count= 330
## x= -10 count= 68
## x= -10 count= 26
## x= -10 count= 72
## x= -10 count= 80
## x= -10 count= 98
## x= -10 count= 52
## x= -10 count= 126
## x= -10 count= 134
## x= -10 count= 14
## x= -10 count= 60
## x= -10 count= 152
## x= -10 count= 30
## x= -10 count= 76
## x= -10 count= 50
## x= -10 count= 136
## x= -10 count= 84
## x= -10 count= 32
## x= -10 count= 46
## x= -10 count= 70
## x= -10 count= 50
## x= -10 count= 188
## x= -10 count= 12
## x= -10 count= 24
## x= -10 count= 34
## x= -10 count= 84
## x= -10 count= 12
## x= -10 count= 64
## x= -10 count= 140
## x= -10 count= 176
## x= -10 count= 68
## x= -10 count= 100
## x= -10 count= 38
## x= -10 count= 44
## x= -10 count= 156
## x= -10 count= 14
## x= -10 count= 78
## x= -10 count= 130
## x= -10 count= 52
## x= -10 count= 132
## x= -10 count= 184
## x= -10 count= 120
## x= -10 count= 286
## x= -10 count= 34
## x= -10 count= 236
## x= -10 count= 56
## x= -10 count= 158
## x= -10 count= 256
## x= -10 count= 134
## x= -10 count= 148
## x= -10 count= 174
## x= -10 count= 16
## x= -10 count= 46
## x= -10 count= 28
## x= -10 count= 46
## x= -10 count= 324
## x= -10 count= 18
## x= -10 count= 48
## x= -10 count= 194
## x= -10 count= 192
## x= -10 count= 66
## x= -10 count= 12
## x= -10 count= 72
## x= -10 count= 248
## x= -10 count= 124
## x= -10 count= 26
## x= -10 count= 32
## x= -10 count= 14
## x= -10 count= 204
## x= -10 count= 80
## x= -10 count= 82
## x= -10 count= 34
## x= -10 count= 86
## x= -10 count= 568
## x= -10 count= 114
## x= -10 count= 68
## x= -10 count= 26
## x= -10 count= 56
## x= -10 count= 80
## x= -10 count= 60
## x= -10 count= 40
## x= -10 count= 24
## x= -10 count= 66
## x= -10 count= 52
## x= -10 count= 40
## x= -10 count= 294
## x= -10 count= 52
## x= -10 count= 48
## x= -10 count= 166
## x= -10 count= 152
## x= -10 count= 272
## 為什麼你通常會輸光? (模擬100次，調整莊家的勝率至501/1000，賭徒與莊家財力 10:100)
for ( i in 1:100) {
x<-0 ; count<-0
repeat{x<- x + sample(c(rep(1,499),rep(-1,501)),1) ; count<-count+1; if(x>=100 || x<=-10) {cat("x=",x,"count=",count, "\n"); break}} }
## x= -10 count= 488
## x= -10 count= 48
## x= -10 count= 64
## x= -10 count= 94
## x= -10 count= 142
## x= -10 count= 284
## x= -10 count= 42
## x= -10 count= 22
## x= -10 count= 140
## x= -10 count= 28
## x= 100 count= 1216
## x= 100 count= 3838
## x= -10 count= 210
## x= -10 count= 50
## x= -10 count= 122
## x= -10 count= 532
## x= -10 count= 32
## x= -10 count= 36
## x= 100 count= 10318
## x= -10 count= 40
## x= -10 count= 2550
## x= -10 count= 170
## x= -10 count= 42
## x= -10 count= 66
## x= -10 count= 146
## x= -10 count= 182
## x= -10 count= 322
## x= -10 count= 1118
## x= -10 count= 52
## x= -10 count= 322
## x= -10 count= 340
## x= -10 count= 56
## x= -10 count= 742
## x= -10 count= 54
## x= 100 count= 1372
## x= -10 count= 2270
## x= -10 count= 140
## x= -10 count= 788
## x= -10 count= 210
## x= -10 count= 254
## x= 100 count= 2030
## x= -10 count= 724
## x= -10 count= 466
## x= -10 count= 108
## x= -10 count= 62
## x= -10 count= 9652
## x= 100 count= 2108
## x= -10 count= 46
## x= 100 count= 3500
## x= -10 count= 318
## x= -10 count= 154
## x= 100 count= 2954
## x= -10 count= 2086
## x= -10 count= 674
## x= -10 count= 1154
## x= -10 count= 246
## x= -10 count= 176
## x= -10 count= 3492
## x= -10 count= 86
## x= -10 count= 16
## x= -10 count= 86
## x= -10 count= 862
## x= -10 count= 26
## x= -10 count= 168
## x= -10 count= 98
## x= -10 count= 394
## x= -10 count= 80
## x= -10 count= 36
## x= -10 count= 116
## x= -10 count= 80
## x= -10 count= 3678
## x= -10 count= 564
## x= -10 count= 418
## x= -10 count= 38
## x= -10 count= 86
## x= -10 count= 906
## x= -10 count= 76
## x= -10 count= 16
## x= -10 count= 20
## x= -10 count= 1354
## x= -10 count= 5756
## x= -10 count= 68
## x= -10 count= 308
## x= -10 count= 1808
## x= -10 count= 1998
## x= -10 count= 360
## x= -10 count= 5186
## x= -10 count= 92
## x= -10 count= 800
## x= -10 count= 2552
## x= -10 count= 108
## x= -10 count= 66
## x= -10 count= 676
## x= -10 count= 108
## x= 100 count= 1228
## x= -10 count= 56
## x= -10 count= 396
## x= -10 count= 408
## x= -10 count= 104
## x= -10 count= 16
## Page38
A<-matrix(c(.5,.4,.2,.3,.3,.4,.2,.3,.4),nrow=3, dimnames=list(c("晴","陰","雨"),c("晴","陰","雨")))
A
## 晴 陰 雨
## 晴 0.5 0.3 0.2
## 陰 0.4 0.3 0.3
## 雨 0.2 0.4 0.4
# Ex: 0.41 = 0.5*0.5+0.3*0.4+0.2*0.2
A%*%A
## 晴 陰 雨
## 晴 0.41 0.32 0.27
## 陰 0.38 0.33 0.29
## 雨 0.34 0.34 0.32
A%*%A%*%A
## 晴 陰 雨
## 晴 0.387 0.327 0.286
## 陰 0.380 0.329 0.291
## 雨 0.370 0.332 0.298
P <- matrix(c(.7,0,.5,0,.3,0,.5,0,0,.4,0,.2,0,.6,0,.8),nrow=4)
row.names(P) <- c("RR","RN","NR","NN")
colnames(P) <- c("RR","RN","NR","NN")
P
## RR RN NR NN
## RR 0.7 0.3 0.0 0.0
## RN 0.0 0.0 0.4 0.6
## NR 0.5 0.5 0.0 0.0
## NN 0.0 0.0 0.2 0.8
# "RRNRNNR" => RR -> RN -> NR -> RN -> NN -> NR
(P1 <- 0.3*0.4*0.5*0.6*0.2)
## [1] 0.0072
# "RRRNNRN" => RR -> RR -> RN -> NN -> NR -> RN
(P2 <- 0.7*0.3*0.6*0.2*0.5)
## [1] 0.0126
## Page45
# Entropy的計算可用定義函數的方式簡化，例如：
# ent <- function(x) {sum(-x*log2(x))}
(H <- -(8*(1/16)*log2(1/16)+24*(1/48)*log2(1/46)))
## [1] 4.761781
# 若32隊奪冠機會完全相同，則由哪對奪冠不確定性高，此時 H=32*-(1/32)*log2(1/32)=5 達最大值)
## Page47
(HX <- -((1/2)*log2(1/2)+(1/4)*log2(1/4)+2*(1/8)*log2(1/8)))
## [1] 1.75
(HY <- -(4*(1/4)*log2(1/4)))
## [1] 2
(HYIX <- -((1/8)*log2(1/4)+2*(1/16)*log2(1/8)+(1/4)*log2(1/2)+(1/8)*log2(1/2)+2*(1/16)*log2(1/4)+2*(1/32)*log2(1/4)+(1/16)*log2(1/2)+2*(1/32)*log2(1/4)+(1/16)*log2(1/2)))
## [1] 1.625
(HXIY <- -((1/8)*log2(1/2)+(1/16)*log2(1/4)+2*(1/32)*log2(1/8)+(1/16)*log2(1/4)+(1/8)*log2(1/2)+2*(1/32)*log2(1/8)+4*(1/16)*log2(1/4)+(1/4)*log2(1)))
## [1] 1.375
(IXY <- -((1/8)*log2(1)+(1/16)*log2(1)+2*(1/32)*log2(1)+(1/16)*log2(1/2)+(1/8)*log2(1/2)+2*(1/32)*log2(1)+(1/16)*log2(1/2)+(1/16)*log2(1)+2*(1/16)*log2(2)+(1/4)*log2(1/2)))
## [1] 0.375
(IXY <- HX-HXIY)
## [1] 0.375
# 另一作法
(IYX <- HY-HYIX)
## [1] 0.375
(IXY==IYX)
## [1] TRUE
# 可看出 I(X;Y) = I(Y;X)
(HXY <- -(2*(1/8)*log2(1/8)+(1/4)*log2(1/4)+6*(1/16)*log2(1/16)+4*(1/32)*log2(1/32)))
## [1] 3.375
(IXY==HX+HY-HXY)
## [1] TRUE
# 可看出 I(X;Y) = H(X)+H(Y)-H(X,Y)
(IXX= ((1/2)*log2((1/2)/(1/4))+(1/4)*log2((1/4)/(1/16))+2*(1/8)*log2((1/8)/(1/64))))
## [1] 1.75
(IXX==HX)
## [1] TRUE
# 可看出 I(X;X) = H(X)
## Page48 (Entropy maximized when p=1/2)
(fx <- deriv(~ -(x*log2(x)+(1-x)*log2(1-x)),"x"))
## expression({
## .expr1 <- log2(x)
## .expr3 <- 1 - x
## .expr4 <- log2(.expr3)
## .expr8 <- log(2)
## .value <- -(x * .expr1 + .expr3 * .expr4)
## .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
## .grad[, "x"] <- -(.expr1 + x * (1/(x * .expr8)) - (.expr3 *
## (1/(.expr3 * .expr8)) + .expr4))
## attr(.value, "gradient") <- .grad
## .value
## })
x <- 0.5
eval(fx)
## [1] 1
## attr(,"gradient")
## x
## [1,] 0
## Page50
A <- matrix(c(1,3,1,2,3,4,0,5,2,1,2,3,2,5,-1,3),nrow=4)
B <- matrix(c(2,1,2,3,1,3,1,4,0,1,-2,2,4,0,0,-1),nrow=4)
(AB <- A%*%B)
## [,1] [,2] [,3] [,4]
## [1,] 15 20 3 2
## [2,] 27 36 12 7
## [3,] 3 -1 -6 5
## [4,] 24 32 5 5
# 將 A、B矩陣相乘並印出結果
# 將A拆成4個小矩陣，分別乘以B
A1B <- A[1,]%*%B
A2B <- A[2,]%*%B
A3B <- A[3,]%*%B
A4B <- A[4,]%*%B
(JA <- rbind(A1B,A2B,A3B,A4B))
## [,1] [,2] [,3] [,4]
## [1,] 15 20 3 2
## [2,] 27 36 12 7
## [3,] 3 -1 -6 5
## [4,] 24 32 5 5
# 將個四乘出來的部分合併
(AB==A%*%B)
## [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE
# 與A、B矩陣相乘結果比較，完全一致
(A23B <- A[2,1]*B[1,3]+A[2,2]*B[2,3]+A[2,3]*B[3,3]+A[2,4]*B[4,3])
## [1] 12
# MapReduce分治
(A23B==AB[2,3])
## [1] TRUE
# 與A、B矩陣相乘結果之第二列第三行元素，完全一致
K <- matrix(c(1,2,2,5),nrow=2)
MATH <- matrix(c(13,1,20,8),nrow=2)
ENC <- K%*%MATH
# 用 K 做 Encryption
KI <- solve(K)
# As a key to Decryption
DEC <- KI%*%ENC
(DEC==MATH)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
# 判斷解碼後結果與原始文字所代表矩陣是否一致
###### 測試 MapReduce分治算法
A <- matrix(rnorm(1000000),nrow=1000)
B <- matrix(rnorm(1000000),nrow=1000)
run1<-function(){
AB <- A%*%B ; print(AB[3,2])
}
run2<-function(){
s<-0
{for (i in 1:nrow(A)) {s <- s + A[3,i]*B[i,2] }
}
print(s)
}
system.time(run1())
## [1] -39.95986
## 使用者 系統 流逝
## 1.36 0.01 1.37
system.time(run2())
## [1] -39.95986
## 使用者 系統 流逝
## 0.02 0.00 0.02
######