※ 引述《[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海