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

作者: 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