→ a78998042a: step1:看下來意思是假設資料是gamma，然後估計他若是 05/04 00:16
→ a78998042a: gamma的參數，再用給定估計參計下的gamma以有母數 05/04 00:16
→ a78998042a: booststrap生成資料 05/04 00:16
→ a78998042a: step2:透過某個model估計出MLE(是gamma的likelihood ? 05/04 00:17
→ a78998042a: step3:optim 05/04 00:17
→ a78998042a: 先不論理解正確與否，但下面寫說，這個過程沒問題， 05/04 00:18
→ a78998042a: 有問題的是程式的輸出結果，不過沒有 真實資料 跟 05/04 00:19
→ a78998042a: 假定model應該沒有辦法重現你的問題? 05/04 00:19
→ a78998042a: 還是請直接提供整個程式?能重現錯誤也不會理解錯誤? 05/04 00:21
→ x88776544pc: 從 NA 來猜的話，那 62 個 a 是負的嗎? 05/04 00:41
因為程式很長，符號很多，所以我大概用文字描述一下
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
推 a78998042a: 你好，謝謝補充，因為沒有足夠重現現象的條件，而過程 05/04 13:29
→ a78998042a: 看起來沒有錯誤，只能猜測一下。 05/04 13:29
→ a78998042a: 懷疑是，你生成的過程是否沒有setseed(set.seed(100)) 05/04 13:30
→ a78998042a: 導致你重覆驗證時，其實是用不同的估計值，所以有了 05/04 13:30
→ a78998042a: 錯誤的誤解?(會這樣懷疑是因為上面寫，NA""大約""有62 05/04 13:30
→ a78998042a: )，是否先將bh做成一個matrix，將所有結果存下來，再 05/04 13:31
→ a78998042a: 重新進行驗證? 05/04 13:31
您好，我確實沒有 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 迴圈的問題
非常感謝您熱心的幫助，我看這一次跑不跑得出來，謝謝。
→ x88776544pc: 如果 code 都沒問題，即迴圈中沒有一次符合 if 條件 05/04 15:22
→ x88776544pc: 你可以試試用 min 取代 <0.000001 去驗證 05/04 15:27
您的意思是找與 0.5 相減的絕對值中，最小的那一個嗎？
確實這樣好像就不用控制誤差，剛剛發現好像是誤差太小了導致程式找不到，才會 NA
※ 編輯: strp (36.236.59.243), 05/04/2017 16:17:54
→ a78998042a: seed是加在整個loop外部，所以每次生成的資料仍不同， 05/04 16:44
→ a78998042a: 目的是重現相同結果。 05/04 16:45
嗯嗯謝謝，我先看看改了誤差值之後能不能跑出來，真的感謝您。
※ 編輯: strp (36.236.59.243), 05/04/2017 19:40:48
結果真的把誤差改大就可以跑了！感謝各位的幫忙！！！
※ 編輯: strp (36.236.59.243), 05/05/2017 09:36:22