※ 引述《[email protected] (羊咩咩)》之銘言： > 請教一下如何利用R這個程式 > 使用accept-reject algorithm來產生常態分配的隨機變數 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 這是學習程式設計而作的練習嗎? 還是因有其他目的而需得 n 的常態分配的樣本? 若是後者, 在 R 用 rnorm(n) 就可以了, 甚至可用 rnorm(n, mean=3, sd=5) 去直接取得平均數為 3 標準差為 5 的常態樣本. > 演算法的步驟為 > 1. 產生 U~uniform(-1,1) > V~uniform(-1,1) > 2. let W=U^2 + V^2 > 3. If W>1 , 回到步驟一重新來一次 > 這樣的話就可以產生隨機變數 服從常態分配了 > 可是在寫程式的時候 要如何寫 > 才能夠在一個function內產生n個服從常態分配的樣本呢?? > 對於軟體方面還是新手 不熟悉該怎麼寫迴圈...>< > 不寫迴圈的話 一次只能產生一個樣本了...Orz ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 若是為了程式設計的練習 [但願不是], 下面這個被包裝 (wrap-up) 為 function 的可供參考: rNorm <- function(n=10) { z <- numeric(n) k <- 1 while(k <= n) { w <- sum(runif(2, -1, 1)^2) if(w <= 1) { z[k] <- w k <- k + 1 } } z } 例如: > rNorm(5) [1] 0.4193 0.9392 0.2459 0.7765 0.2454 > rNorm(10) [1] 0.5804 0.2529 0.7186 0.6608 0.1869 0.6637 0.3839 0.6172 0.5858 0.4746 只是這麼作非常沒有效率! 用下面這個小 function 去比較 nNorm() 與 R 內設提供的 rnorm() 在產生較大數量的樣本數時其耗費的時間差 (單位是秒) compN2n <- function(n) { system.time(rNorm(n))[3] - system.time(rnorm(n))[3] } > compN2n(1000) [1] 0.03 > compN2n(5000) [1] 0.13 > compN2n(10000) [1] 0.25 > compN2n(50000) [1] 1.2 若只需要產生一次大小為 n 的常態樣本, 那也就罷了; 但若這是為一個重複執行的模擬所需, 那就很沒效率了. Accept-reject algorithm 是個蠻重要的方法, 但是去練習這方法的題目或許應可更具有實用價值會比較好. -- 夫兵者不祥之器物或惡之故有道者不處君子居則貴左用兵則貴右兵者不祥之器非君子 之器不得已而用之恬淡為上勝而不美而美之者是樂殺人夫樂殺人者則不可得志於天下 矣吉事尚左凶事尚右偏將軍居左上將軍居右言以喪禮處之殺人之眾以哀悲泣之戰勝以 喪禮處之道常無名樸雖小天下莫能臣侯王若能守之萬物將自賓天地相合以降甘露民莫 之令而自均始制有名名亦既有夫亦將知止知止可以不殆譬道之在天下210.66.0.209海