用S作隨機模擬計算 作為統計工作者，我們除了可以用S迅速實現新的統計方法，還可以用S進行隨機模擬。隨機模擬可以驗證我們的演算法、比較不同演算法的的優缺點、發現改進統計方法的方向，是進行統計研究的最有力的計算工具之一。 隨機模擬最基本的需要是產生偽隨機數，S中已提供了大多數常用分佈的偽隨機數函式，可以返回一個偽隨機數序列向量。這些偽隨機數函式以字母r開頭，比如rnorm()是正態偽隨機數函式，runif()是均勻分佈偽隨機數函式，其第一個自變數是偽隨機數序列長度n。關於這些函式可以參見第14節以及系統幫助檔案。下例產生1000個標準正態偽隨機數： > y <- rnorm(1000) 這些偽隨機數函式也可以指定與分佈有關的引數，比如下例產生1000個均值為150、標準差為100的正態偽隨機數： > y <- rnorm(1000, mean=150, sd=100) 產生偽隨機數序列是不重複的，實際上，S在產生偽隨機數時從一個種子出發，不斷迭代更新種子，所以產生若干隨機數後內部的隨機數種子就已經改變了。有時我們需要模擬結果是可重複的，這隻要我們儲存當前的隨機數種子，然後在每次產生偽隨機數序列之前把隨機數種子置為儲存值即可： > the.seed <- .Random.seed > …………… > .Random.seed <- the.seed > y <- rnorm(1000) 作為例子，我們來產生服從一個簡單的線性迴歸的資料。 # 簡單線性迴歸的模擬 lm.simu <- function(n){ # 先生成自變數。假設自變數x的取值範圍在150到180之間，大致服從正態分佈。 x <- rnorm(n, mean=165, sd=7.5) # 再生成模型誤差。假設服從N(0, 1.2)分佈 eps <- rnorm(n, 0, 1.2) # 用模型生成因變數 y <- 0.8 * x + eps return(data.frame(y,x)) } S沒有提供多元隨機變數的模擬程式，這裡給出一個進行三元正態隨機變數模擬的例子。假設要三元正態隨機向量 的 n個獨立觀測，可以先產生n個服從三元標準正態分佈的觀測，放在一個 n行3列的矩陣中： U <- matrix(rnorm(3*n), ncol=3, byrow=T) 可以認為矩陣U的每一行是一個標準的三元正態分佈的觀測。設矩陣 的Choleski分解為 ， A為上三角矩陣，若隨機向量 ，則 。因此， 作為一個三行 n列的矩陣每一行都是服從 分佈的，且各行之間獨立。經過轉置，產生的 X X <- matrix(rep(mu,n), ncol=3, byrow=T) + U %*% A 是一個 n行三列的矩陣。 有時模擬需要的計算量很大，多的時候甚至要計算幾天的時間。對於這種問題我們要善於把問題拆分成可以單獨計算的小問題，然後單獨計算每個小問題，把結果儲存在S物件中或文字檔案中，最後綜合儲存的結果得到最終結果。 如果某一個問題需要的計算時間比較長，我們在程式設計時可以採用以下的技巧：每隔一定時間就顯示一下任務的進度，以免計算已經出錯或進入死迴圈還不知道；應該把中間結果每隔一段時間就記錄到一個文字檔案中（cat()函式可以帶一個file引數和append引數，對這種記錄方法提供了支援），如果需要中斷程式，中間結果可能是有用的，有些情況下還可以根據記錄的中間結果從程式中斷的地方繼續執行。 參考文獻： http://www.math.pku.edu.cn/teachers/lidf/docs/statsoft/html/s/13.html
用S作隨機模擬計算 作為統計工作者，我們除了可以用S迅速實現新的統計方法，還可以用S進行隨機模擬。隨機模擬可以驗證我們的演算法、比較不同演算法的的優缺點、發現改進統計方法的方向，是進行統計研究的最有力的計算工具之一。 隨機模擬最基本的需要是產生偽隨機數，S中已提供了大多數常用分佈的偽隨機數函式，可以返回一個偽隨機數序列向量。這些偽隨機數函式以字母r開頭，比如rnorm()是正態偽隨機數函式，runif()是均勻分佈偽隨機數函式，其第一個自變數是偽隨機數序列長度n。關於這些函式可以參見第14節以及系統幫助檔案。下例產生1000個標準正態偽隨機數： > y <- rnorm(1000) 這些偽隨機數函式也可以指定與分佈有關的引數，比如下例產生1000個均值為150、標準差為100的正態偽隨機數： > y <- rnorm(1000, mean=150, sd=100) 產生偽隨機數序列是不重複的，實際上，S在產生偽隨機數時從一個種子出發，不斷迭代更新種子，所以產生若干隨機數後內部的隨機數種子就已經改變了。有時我們需要模擬結果是可重複的，這隻要我們儲存當前的隨機數種子，然後在每次產生偽隨機數序列之前把隨機數種子置為儲存值即可： > the.seed <- .Random.seed > …………… > .Random.seed <- the.seed > y <- rnorm(1000) 作為例子，我們來產生服從一個簡單的線性迴歸的資料。 # 簡單線性迴歸的模擬 lm.simu <- function(n){ # 先生成自變數。假設自變數x的取值範圍在150到180之間，大致服從正態分佈。 x <- rnorm(n, mean=165, sd=7.5) # 再生成模型誤差。假設服從N(0, 1.2)分佈 eps <- rnorm(n, 0, 1.2) # 用模型生成因變數 y <- 0.8 * x + eps return(data.frame(y,x)) } S沒有提供多元隨機變數的模擬程式，這裡給出一個進行三元正態隨機變數模擬的例子。假設要三元正態隨機向量 的 n個獨立觀測，可以先產生n個服從三元標準正態分佈的觀測，放在一個 n行3列的矩陣中： U <- matrix(rnorm(3*n), ncol=3, byrow=T) 可以認為矩陣U的每一行是一個標準的三元正態分佈的觀測。設矩陣 的Choleski分解為 ， A為上三角矩陣，若隨機向量 ，則 。因此， 作為一個三行 n列的矩陣每一行都是服從 分佈的，且各行之間獨立。經過轉置，產生的 X X <- matrix(rep(mu,n), ncol=3, byrow=T) + U %*% A 是一個 n行三列的矩陣。 有時模擬需要的計算量很大，多的時候甚至要計算幾天的時間。對於這種問題我們要善於把問題拆分成可以單獨計算的小問題，然後單獨計算每個小問題，把結果儲存在S物件中或文字檔案中，最後綜合儲存的結果得到最終結果。 如果某一個問題需要的計算時間比較長，我們在程式設計時可以採用以下的技巧：每隔一定時間就顯示一下任務的進度，以免計算已經出錯或進入死迴圈還不知道；應該把中間結果每隔一段時間就記錄到一個文字檔案中（cat()函式可以帶一個file引數和append引數，對這種記錄方法提供了支援），如果需要中斷程式，中間結果可能是有用的，有些情況下還可以根據記錄的中間結果從程式中斷的地方繼續執行。 參考文獻： http://www.math.pku.edu.cn/teachers/lidf/docs/statsoft/html/s/13.html