迴圈控制 : for, while, repeat Textbook reading: Chapter 7. 固定次數迴圈 : for 指令. 假設有 k 個指令 C(1),..., C(k), 而我們要依序完成其中的 C(m),..., C(n), 語法為 for ( i in m:n) { C(i) 通常在 for 迴圈開始前, 會定義一個存結果的物件, 而在迴圈執行時更新物件. Example 1. Suppose that a sequence {a i i=1 is defined by a 1 = 2 and a i+1 = a i for i 1. (a) Write down R commands that compute a 10. (b) Write down R commands that compute the vector (a 1,..., a 10 ). Solutions. #(a) a <- 2 for (i in 1:9){ a <- sqrt(a) a #(b) a <- rep(2, 10) for (i in 1:9){ a[i+1] <- sqrt(a[i]) a 利用迴圈讓函數輸出向量. Example 2. 而 令 F 為 N(0, 1)CDF, 即 F (x) = x f(t)dt, x (, ) f(t) = e t2 /2 2π, t (, ). 執行以下 R 命令可得 F (0), F (1), F (2): f <- function(x){ exp(-x^2/2)/sqrt(2*pi) F <- function(x){ ans <- integrate(f, -Inf, x)$value #ans 為 f 從負無窮大到的 x 積分值 F(0);F(1);F(2) 1
但是 F(0:2) 無法執行. 寫下一個 R function F.fun, 當 input 為 (x 1,..., x n ) 時, F.fun 的 output 為 (F (x 1 ),..., F (x n )). 執行 F.fun(0:2) 以測試 F.fun. Solution. F.fun <- function(x){ ans <- x #define a vector to store the results n <- length(x) for (i in 1:n){ ans[i] <- F(x[i]) #update ans F.fun(0:2) 注意 : 使用函數繪圖指令 curve(f,a,b) 時, 函數 f 必須能輸出向量. curve(f,-1,1) #get error curve(f.fun,-1,1) Example 3. 寫一個函數 m.prod 計算矩陣相乘. m.prod <- function(a, B){ m <- dim(a)[1] n <- dim(b)[2] if (dim(b)[1]!=dim(a)[2]) return("error!") D <- matrix(0, m, n) ## 定義一個 m x n 的零矩陣 for (i in 1:m){ for (j in 1:n){ D[i,j] <- sum(a[i,]*b[,j]) return(d) 測試 A <- matrix(1:6, 3,2) B <- t(a) A %*% B m.prod(a, B) 畫 z = f(x, y) 的函數圖形, 其中 x {x 1,..., x m, y {y 1,..., y n. R 指令 : persp(xvec, yvec, z) xvec 為向量 (x 1,..., x m ) 2
yvec 為向量 (y 1,..., y n ) z 為 m n 矩陣, z 的第 (i, j) 位置為 f(x i, y j ). Example 4. Plot z = x 2 +y 2 for x { 0.9, 0.8,... 0.9, y { 2.20, 2.19,..., 2.20 using R command persp. x <- seq(-0.9, 0.9, by=0.1) y <- seq(-2.20, 2.20, by=0.01) m <- length(x) n <- length(y) z <- matrix(0, m, n) for (i in 1:m){ for (j in 1:n){ z[i,j] <- x[i]^2+y[j]^2 persp(x,y,z) #You can adjust the viewing angle #by changing the parameters theta and phi in persp. persp(x,y,z, theta=30, phi=15) persp(x,y,z, theta=0, phi=15) persp(x,y,z, theta=0, phi=30) 使用迴圈定義具有規律型態的矩陣 Example 5. Write down R commands that define a 100 100 matrix M, where M is 1 2 3... 100 0 1 2... 99 0 0 1... 98. 0 0 0... 1 Solution. n <- 100 M <- matrix(0, n, n) for (i in 1:n){ M[i,i:n] <- 1:(n-i+1) M R 中使用迴圈比較慢. 所以如果能避免使用迴圈就避免. 以下考慮二個計算絕對值並可輸出向量的函數 abs1.fun 和 abs2.fun. 3
函數 abs1.fun 使用迴圈和函數 abs0.fun, 而 abs0.fun 計算絕對值時只能輸出一個值 : abs0.fun <- function(x){ if (x<0) { return(-x) return(x) abs1.fun <- function(x){ ans <- x #define a vector to store the results n <- length(x) for (i in 1:n){ ans[i] <- abs0.fun(x[i]) #update ans 函數 abs2.fun 未使用迴圈也可計算絕對值並可輸出向量 : abs2.fun <- function(x){ ans <- x ans2 <- -x ans[x<0] <- ans2[x<0] 測試計算時間 x <- 1:5000000 t1 <- proc.time() y1 <- abs1.fun(x) t1 <- proc.time() - t1 t2 <- proc.time() y2 <- abs2.fun(x) t2 <- proc.time() - t2 t1 t2 在 R 中內建計算絕對值的函數叫做 abs, 可輸出向量. 當迴圈次數不固定, 而是在條件下停止執行時, 可使用 while 或 repeat. 假設當條件 A 為 TRUE 時要執行工作 B, 而在條件 A 為 FALSE 時停止執行. while 語法為 while(a){ B 4
假設要重複執行工作 B, 直到條件 A 為 TRUE 時停止執行. repeat 語法為 repeat{ B if (A) break Example 6. 重複執行 add 函數, 每次詢問使用者是否繼續, 直到使用者輸入 N 時停止執行. 先定義 add 函數和詢問字串 S. add <- function(){ x <- readline("please enter a number: \n") y <- readline("please enter another number: \n") qn <- paste("please enter the sum of", x,"and", y, ":\n", sep=" ") z <- readline(qn) x <- as.numeric(x) y <- as.numeric(y) z <- as.numeric(z) if (z==(x+y)) { res <- "Your answer is correct. Good job!" else { res <- paste("the sum of",x,"and", y,"should be",x+y, "!", sep=" ") cat(res, "\n") S <- "Do you want to play this game again? \n" S <- paste(s, "Enter N to stop or ", sep="") S <- paste(s, "enter something else to play again: \n", sep="") 用 repeat 加上 break 寫法如下 : repeat{ add() conti <- readline(prompt=s) if (conti=="n") break 用 while 寫法如下 : conti <- "Y" while( conti!="n" ){ add() conti <- readline(prompt=s) 5
Example 7. 以牛頓法計算 x 2 2 = 0 的正根. 取起始解為 1. 當本次疊代的解和前一次疊代的解相差 < 0.001 時停止疊代. 上述方法疊代公式如下. 令 x n 為第 n 次的疊代解, 則 x 1 = 1 且 用 while 寫法如下 : x n+1 = x n x2 n 2 2x n for n 1. x <- 1 d <- 100 #d can be any number >= 0.001 while ( d >= 0.001 ) { x.new <- x - (x^2-2)/(2*x) d <- abs(x.new - x) x <- x.new x #x 用 repeat 加上 break 寫法如下 : x <- 1 repeat { x.new <- x - (x^2-2)/(2*x) d <- abs(x.new - x) x <- x.new if (d < 0.001) break x #x 6