作者ritajen (asdfge)
看板R_Language
标题[问题] 关於蒙地卡罗统计模拟
时间Sat Dec 12 12:04:30 2015
[问题类型]
我想用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))
--
Sent from my Windows
--
※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 110.30.135.65
※ 文章网址: https://webptt.com/cn.aspx?n=bbs/R_Language/M.1449893072.A.5D8.html
1F：推 celestialgod: 你的seq出来是向量，你想存在一个的位置，当然出错 12/12 12:09
2F：→ celestialgod: 解决方法是对n回圈 12/12 12:10
3F：→ celestialgod: for (j in seq_along(n)){ 12/12 12:10
4F：推 celestialgod: final(n[j]) 12/12 12:11
5F：→ celestialgod: } 12/12 12:11
6F：→ celestialgod: l95,r95,l99,r99请放到function内 做return 12/12 12:12
7F：推 celestialgod: 不用做return....多打的...不过你最後两个mean要记 12/12 12:13
8F：推 celestialgod: 得用c合并再一起return 12/12 12:13
9F：→ ritajen: lt95 rt95 lt99 rt99这个放到function内，那for回圈还可 12/12 12:20
10F：→ ritajen: 以重复试验1000次吗? 还是"多加"在funtion中，然後最後 12/12 12:20
11F：→ ritajen: 两个mean要合并後return到function吗? 谢谢 12/12 12:20
12F：→ celestialgod: 有电脑再打详细程式码..... 12/12 12:27
13F：→ ritajen: ok 谢谢 12/12 12:32
15F：→ celestialgod: sted 12/12 12:32
16F：→ celestialgod: 没有排版就抱歉了.... 12/12 12:34
18F：→ ritajen: 点疑问，不晓得原因为何? 12/12 13:22
19F：→ celestialgod: 改好了 在上面的pastebin 12/12 13:39
20F：→ celestialgod: no.two不知道怎麽跑出来的，我明明打no.two的.... 12/12 13:40
21F：推 celestialgod: result[i,]才对，维度未考虑清楚，抱歉 12/12 13:42
22F：→ ritajen: 可以了，谢谢 可以请问一下，所以function自创函数的功 12/12 14:05
23F：→ ritajen: 能，就是可以同时跑出多笔结果吗? 12/12 14:05
24F：→ celestialgod: 自创函数是帮助你避免一再重复的程式码，并且增加程 12/12 14:12
25F：→ celestialgod: 式易读性 12/12 14:12
26F：→ ritajen: 了解，万分感谢^^ 12/12 15:43
27F：→ celestialgod: 干，为什麽rep还是变成two... 12/13 15:45
28F：→ ritajen: 对阿 依然会是two，自己改过就没问题了 12/13 22:21