[問題] 關於蒙地卡羅統計模擬

[ritajen (asdfge)]

[問題類型]
我想用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://www.ptt.cc/bbs/R_Language/M.1449893072.A.5D8.html
推 celestialgod: 你的seq出來是向量，你想存在一個的位置，當然出錯
→ celestialgod: 解決方法是對n迴圈
→ celestialgod: for (j in seq_along(n)){
推 celestialgod: final(n[j])
→ celestialgod: }
→ celestialgod: l95,r95,l99,r99請放到function內 做return
推 celestialgod: 不用做return....多打的...不過你最後兩個mean要記
推 celestialgod: 得用c合併再一起return
→ ritajen: lt95 rt95 lt99 rt99這個放到function內，那for迴圈還可
→ ritajen: 以重複試驗1000次嗎? 還是"多加"在funtion中，然後最後
→ ritajen: 兩個mean要合併後return到function嗎? 謝謝
→ celestialgod: 有電腦再打詳細程式碼.....
→ ritajen: ok 謝謝
→ celestialgod: http://pastebin.com/CY84tqmv 車上用手機打的 unte
→ celestialgod: sted
→ celestialgod: 沒有排版就抱歉了....
→ ritajen: 感謝您~http://pastebin.com/J2DgZyj0 裡面標記的部分有
→ ritajen: 點疑問，不曉得原因為何?
→ celestialgod: 改好了 在上面的pastebin
→ celestialgod: no.two不知道怎麼跑出來的，我明明打no.two的....
推 celestialgod: result[i,]才對，維度未考慮清楚，抱歉
→ ritajen: 可以了，謝謝  可以請問一下，所以function自創函數的功
→ ritajen: 能，就是可以同時跑出多筆結果嗎?
→ celestialgod: 自創函數是幫助你避免一再重複的程式碼，並且增加程
→ celestialgod: 式易讀性
→ ritajen: 了解，萬分感謝^^
→ celestialgod: 幹，為什麼rep還是變成two...
→ ritajen: 對阿 依然會是two，自己改過就沒問題了