代码可以写得简洁一些,用R编程序很多时候都可以避免显式的循环,比如:
f=function(){
x=15+2*rnorm(10)
m=mean(x); s=sd(x)
(15<m+.579*s) & (15>m-.579*s)
}
sum(replicate(100, f()))
不知道你是不是想画置信区间的图,试试这段代码:
<br />
f = function(n = 100, a = 15, b = 2, z = 0.579, rn = 10) {<br />
d = replicate(n, a + b * rnorm(rn))<br />
m = colMeans(d)<br />
s = apply(d, 2, sd)<br />
y0 = m - z * s<br />
y1 = m + z * s<br />
plot(1, xlim = c(0.5, n + 0.5), ylim = c(min(y0), max(y1)),<br />
type = "n", xlab = "", ylab = "")<br />
for (i in 1:n) {<br />
arrows(i, y0[i], i, y1[i], length = 0.1, angle = 90,<br />
code = 3, col = ifelse(a > y0[i] & a < y1[i], "blue",<br />
"red"))<br />
points(i, m[i])<br />
}<br />
abline(h = a, lty = 2)<br />
Sys.sleep(2)<br />
}<br />
<br />
replicate(20, f())
n是产生随机数(一组一组的)的组数,a+b*rnorm()就不说了,rn是每组随机数的个数,z是标准差前面乘的系数。这幅图用了箭头的一个小伎俩:把箭尾的角度设为90度,当然不一定用90度。Sys.sleep(2)停顿两秒是为了让你看清模拟的过程,你也可以在循环的过程中加上Sys.sleep()。