真是夠烏龍的 ... ※ 引述《cjlu (Joe)》之銘言： > ※ 引述《[email protected] (羊咩咩)》之銘言： > > 演算法的步驟為 > > 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 是個蠻重要的方法, > 但是去練習這方法的題目或許應可更具有實用價值會比較好. -- 夫兵者不祥之器物或惡之故有道者不處君子居則貴左用兵則貴右兵者不祥之器非君子 之器不得已而用之恬淡為上勝而不美而美之者是樂殺人夫樂殺人者則不可得志於天下 矣吉事尚左凶事尚右偏將軍居左上將軍居右言以喪禮處之殺人之眾以哀悲泣之戰勝以 喪禮處之道常無名樸雖小天下莫能臣侯王若能守之萬物將自賓天地相合以降甘露民莫 之令而自均始制有名名亦既有夫亦將知止知止可以不殆譬道之在 140.116.152.108海