[問題] 可靠度分析求 tp 值
[問題類型]:
程式諮詢(我想用 R 做某件事情，程式已經寫完了！但有一小部分覺得很奇怪)
[軟體熟悉度]:
入門(寫過其他程式，只是對語法不熟悉)
[問題敘述]:
我的目的是利用 for 迴圈，找出符合設定式子的 tp 值
p 可以是 0.01, 0.05, 0.1, 0.5，所以會有 t0.01, t0.05, t0.1, t0.5 四個結果
因為只有 t0.5 的部分會出問題，所以以下僅針對 t0.5 來做書寫
式子如下：
(其中 q, a, b 均為 MLE，且為已知)
tp0.5 <- NULL
for (tp in 1:1400000) {
if ( abs( pgamma(q, a*(tp*0.0000001)^b, rate = 1) -0.5 ) <= 0.000001 ) {
tp0.5 <- tp*0.0000001
}
}
tp 會設計成這樣，主要是把區間割的細一點，讓 pgamma 的值逐漸變大
大到與 0.5 的誤差小於等於 0.000001 時，輸出該 tp*0.0000001 值
就可以算出 t0.5 了！
以上單獨的求法是不會有問題的，可以算出 t0.5
但我真正想求的，是透過 1000 次 Bootstrap，來得到 t0.5 的信賴區間
每一次 Bootstrap 的步驟如下
Step1：用真實資料的 MLE 去生成一筆資料 (其中 MLE 已知)
Step2：再利用 Step1 生成的資料，透過我的模型，去估計 MLE*
Step3：利用 Step2 的 MLE* 去算 t0.5 (然後蓋到迴圈外的 NULL 向量去)
所以當 Bootstrap 跑完後，會得到 1000 個 t0.5
使用的程式與以上範例相同，僅修改迴圈內的名稱而已
但是！！！
最後發現...這 1000 個 t0.5 裡面，有些值是 NA (大約有 62 個)
也就是，那幾次迴圈內，並沒有求出 t0.5
我有將這些產生 NA 的 MLE* 抓出來看，都沒有問題，都有值
然後再利用這些 MLE* 重新在用上面提供的程式去計算 t0.5，發現還是求不出來
但奇怪的是，如果我不用迴圈條件去逼近 t0.5，而是單純直接給定 tp
卻都算得出來，因為我們想要的 t0.5 會在 0.11~0.13 之間
這區間內大大小小的值我都給過了，都滿 ok 的
但是在 Bootstrap 迴圈下就是求不出來！這讓我匪夷所思
[程式範例]:
如同以上例子提供，只是該程式會放在 Bootstrap 迴圈內
[環境敘述]:
與此問題無關
[關鍵字]:
tp, t0.5, Bootstrap, MLE
再麻煩各位版友抽空解惑，萬分感謝。
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 36.236.59.243
※ 文章網址: https://www.ptt.cc/bbs/R_Language/M.1493794193.A.5EA.html
→
05/04 00:16, , 1F
05/04 00:16, 1F
→
05/04 00:16, , 2F
05/04 00:16, 2F
→
05/04 00:16, , 3F
05/04 00:16, 3F
→
05/04 00:17, , 4F
05/04 00:17, 4F
→
05/04 00:17, , 5F
05/04 00:17, 5F
→
05/04 00:18, , 6F
05/04 00:18, 6F
→
05/04 00:19, , 7F
05/04 00:19, 7F
→
05/04 00:19, , 8F
05/04 00:19, 8F
→
05/04 00:21, , 9F
05/04 00:21, 9F
→
05/04 00:41, , 10F
05/04 00:41, 10F
因為程式很長，符號很多，所以我大概用文字描述一下
tp0.5 <- NULL ### 先設定一個 NULL 向量，給每一次迴圈的 t0.5 去蓋
for (bs in 1:1000) { ### 開始做 1000 次 Bootstrap
mlerg(h,21,13) ### mlerg 是用來生成資料的函數 (已在之前建構)
### h 是真實資料的 MLE，在前面已經估完了
### 總共生成 21 條 process，每條 13 個特徵值
optim(c(1,1,1,1), like) ### c(1,1,1,1) 為起始值
### like 是我模型的 likelihood function
### 這個 optim 在程式上會重複 12 次到收斂為止
### 收斂之後的 MLE 值會命名成 bh，共四個參數 (a, b, d, B)
for (t0.5 in 1:1400000) {
if ( abs( pgamma( ( log(1.6/0.9) )^( bh[3]+1 )/( (bh[3]+1)*bh[4] ),
bh[1]*(t0.5* 0.0000001)^bh[2]) - 0.5) <= 0.000001 ) {
tp0.5[bs] <- t0.5* 0.0000001
}
}
### 將時間分割成 1000000 塊，從 0.0000001 開始跑，每次增加 0.0000001
讓 pgamma 的結果慢慢逼近 0.5，但考慮到真實數字可能無法完全等於 0.5
所以讓與 0.5 的誤差小於等於 0.000001
### 所以每一次 BS 就是：
Step 1 ：用真實 MLE 生成資料
Step 2 ：用生成的資料求出 MLE*，所以每次 bs 迴圈就會有一組 MLE*
Step 3 ：用 MLE* 去求 tp0.5
### 真實的 MLE ：(a, b, d, B) = (676.77, 0.3329, -0.7199, 0.009083)
} ### bs 迴圈結束
程式如同以上敘述
問題就出在多數 tp0.5 是可以算出來的，1000 次裡面會有 62 次產生 NA
但是如果我直接用以下程式
pgamma( ( log(1.6/0.9) )^( bh[3]+1 )/( (bh[3]+1)*bh[4] ),
bh[1]*( 0.12 )^bh[2])
^
|
|
|
這裡直接帶 0.12 就可以，因為預期的結果 tp0.5 大概就是接近 0.12
也有試過 0.11，0.13, 0.115, ..., 很多值，發現都可以算出來
所以總結來說我覺得是迴圈的問題
不過 tp0.01, tp0.05, tp0.1 都可以算出來，都不會有 NA
唯獨 tp0.5 會出問題
如果是求 tp0.01，則程式如下：
for (t0.5 in 1:1400000) {
if ( abs( pgamma( ( log(1.6/0.9) )^( bh[3]+1 )/( (bh[3]+1)*bh[4] ),
bh[1]*(t0.5* 0.0000001)^bh[2]) - 0.99) <= 0.000001 ) {
tp0.5[bs] <- t0.5* 0.0000001 ^
} |
} |
|
|
|
只有這邊會不一樣
假設 p 為 0.01，這邊會用 1-p，是公式推導的結果
同理，t0.05 時那邊就會放 0.95、t0.1 時就會放 0.9
這 1000 個 tp0.5 是在做 Bootstrap 產生的
在迴圈外，我也會用真實資料的 MLE 去算真實的 tp0.5 (這邊不會有問題)
所以總結來說，感覺是 bs 迴圈內的某些因素導致 tp0.5 算不出來
我看過每一次 bs 的 MLE*，每次都有值，代表資料生成並沒有錯誤
而這些錯誤的 tp0.5 對應的 MLE*，相較其他沒錯的 MLE*，也沒有什麼顯著差別
所以我覺得會不會是程式運算上出現 bug？
※ 編輯: strp (36.236.59.243), 05/04/2017 11:02:29
推
05/04 13:29, , 11F
05/04 13:29, 11F
→
05/04 13:29, , 12F
05/04 13:29, 12F
→
05/04 13:30, , 13F
05/04 13:30, 13F
→
05/04 13:30, , 14F
05/04 13:30, 14F
→
05/04 13:30, , 15F
05/04 13:30, 15F
→
05/04 13:31, , 16F
05/04 13:31, 16F
→
05/04 13:31, , 17F
05/04 13:31, 17F
您好，我確實沒有 set.seed，因為我想說 Bootstrap 的概念
不就是抽一組樣本後 (這邊是選擇固定的 MLE) 反覆做同一件事情
所以我想每次 bs 生成的資料都不一樣、MLE* 也都不一樣好像也滿～～～合理的
因為如果生的資料都一樣，MLE* 不就會一模一樣了嗎？
如果您是說真實 MLE 值的話，在程式內我是直接給定的，因為在其他程式已經算完了
所以不會有真實 MLE 浮動的問題
會產生隨機效應的地方，只有「生成資料」與「估計 MLE*」這邊
如有誤會請見諒，因為我其實很少用 set.seed 這個 code
因為這個程式跑一次就要十幾個小時，所以也沒有詳細的去確認有幾個 NA
貼上來的這一次我確定是 62 個沒錯，因為我還有再把這 62 個抓出來再跑一次
結果還是全部都不行這樣 ...
您說的 set.seed 主要是要來確認是否每次都會是 62 個對嗎？
剛剛我試著把條件改為 abs( pgamma() - 0.5 ) < 0.00001 (原本為 0.000001)
發現這 62 個可以跑出來了！正在重跑程式中，但是我沒有 set.seed 不好意思
我有將 bh 做成一個 matrix，叫 allmlelook，發文之前有先確認了！
每一次的 bs 迴圈的 MLE* 值都有存在，所以才會覺得是 tp0.5 迴圈的問題
非常感謝您熱心的幫助，我看這一次跑不跑得出來，謝謝。
→
05/04 15:22, , 18F
05/04 15:22, 18F
→
05/04 15:27, , 19F
05/04 15:27, 19F
您的意思是找與 0.5 相減的絕對值中，最小的那一個嗎？
確實這樣好像就不用控制誤差，剛剛發現好像是誤差太小了導致程式找不到，才會 NA
※ 編輯: strp (36.236.59.243), 05/04/2017 16:17:54
→
05/04 16:44, , 20F
05/04 16:44, 20F
→
05/04 16:45, , 21F
05/04 16:45, 21F
嗯嗯謝謝，我先看看改了誤差值之後能不能跑出來，真的感謝您。
※ 編輯: strp (36.236.59.243), 05/04/2017 19:40:48
結果真的把誤差改大就可以跑了！感謝各位的幫忙！！！
※ 編輯: strp (36.236.59.243), 05/05/2017 09:36:22
R_Language 近期熱門文章
PTT數位生活區 即時熱門文章
