[问题] 关于蒙地卡罗统计模拟

楼主: ritajen (asdfge)   2015-12-12 12:04:30
[问题类型]
我想用R做统计模拟,看看重复试验后的信赖区间是否名符其实
[软件熟悉程度]
新手
[问题叙述]
已知期望值、标准差随机抽取n个样本后,重复1000,想检查其95%的覆蓋率是否属实,这可以使用for循环得到解,我的问题是如果我的抽取样本数变成不是固定的,如:
5,10,15,20,…,95,100个样本,这样我是可以利用"function"得到结果吗? 如果是以下是我目前的程式,但结果输出后系统出现警示
"In r95[i] <- mean(x) + qnorm(0.975) * sqrt(sigma^2/n) :
被替换的项目不是替换值长度的倍数"
我猜测是function有问题,不过不晓得应该如何解决?
再者,我如果要将结果画成图横轴为sample size,纵轴为覆蓋率,是否应该利用plot的方式进行?
谢谢
[程式范例]
rm(list = ls())
mu <- 7; sigma <- 2; n <- seq(from=5,to=100,by=5); no.rep <- 1000
l95 <- rep(NA,no.rep)
r95 <- rep(NA,no.rep)
l99 <- rep(NA,no.rep)
r99 <- rep(NA,no.rep)
final=function(n){
for(i in 1:no.rep){ #重复1000次
print(i)
set.seed(i)
x <- rnorm(n,mu,sigma)
l95[i] <- mean(x)-qnorm(0.975)*sqrt(sigma^2/n)
r95[i] <- mean(x)+qnorm(0.975)*sqrt(sigma^2/n)
l99[i] <- mean(x)-qnorm(0.995)*sqrt(sigma^2/n)
r99[i] <- mean(x)+qnorm(0.995)*sqrt(sigma^2/n)
}
mean((l95<=mu) & (mu<=r95)) # 检查覆蓋率(coverage)
mean((l99<=mu) & (mu<=r99))
}
final(seq(from=5,to=100,by=5))
作者: celestialgod (天)   2015-12-12 12:09:00
你的seq出来是向量,你想存在一个的位置,当然出错解决方法是对n循环for (j in seq_along(n)){final(n[j])}l95,r95,l99,r99请放到function内 做return不用做return....多打的...不过你最后两个mean要记得用c合并再一起return
楼主: ritajen (asdfge)   2015-12-12 12:20:00
lt95 rt95 lt99 rt99这个放到function内,那for循环还可以重复试验1000次吗? 还是"多加"在funtion中,然后最后两个mean要合并后return到function吗? 谢谢
作者: celestialgod (天)   2015-12-12 12:27:00
有电脑再打详细程式码.....
楼主: ritajen (asdfge)   2015-12-12 12:32:00
ok 谢谢
作者: celestialgod (天)   2015-12-12 12:32:00
http://pastebin.com/CY84tqmv 车上用手机打的 untested没有排版就抱歉了....
楼主: ritajen (asdfge)   2015-12-12 13:22:00
感谢您~http://pastebin.com/J2DgZyj0 里面标记的部分有点疑问,不晓得原因为何?
作者: celestialgod (天)   2015-12-12 13:39:00
改好了 在上面的pastebinno.two不知道怎么跑出来的,我明明打no.two的....result[i,]才对,维度未考虑清楚,抱歉
楼主: ritajen (asdfge)   2015-12-12 14:05:00
可以了,谢谢 可以请问一下,所以function自创函数的功能,就是可以同时跑出多笔结果吗?
作者: celestialgod (天)   2015-12-12 14:12:00
自创函数是帮助你避免一再重复的程式码,并且增加程式易读性
楼主: ritajen (asdfge)   2015-12-12 15:43:00
了解,万分感谢^^
作者: celestialgod (天)   2015-12-13 15:45:00
干,为什么rep还是变成two...
楼主: ritajen (asdfge)   2015-12-13 22:21:00
对阿 依然会是two,自己改过就没问题了

Links booklink

Contact Us: admin [ a t ] ucptt.com