作者strp (GlRoEvEeN)
看板R_Language
标题[问题] 可靠度分析求 tp 值
时间Wed May 3 14:49:50 2017
[问题类型]:
程式谘询(我想用 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://webptt.com/cn.aspx?n=bbs/R_Language/M.1493794193.A.5EA.html
1F：→ a78998042a: step1:看下来意思是假设资料是gamma，然後估计他若是 05/04 00:16
2F：→ a78998042a: gamma的参数，再用给定估计参计下的gamma以有母数 05/04 00:16
3F：→ a78998042a: booststrap生成资料 05/04 00:16
4F：→ a78998042a: step2:透过某个model估计出MLE(是gamma的likelihood ? 05/04 00:17
5F：→ a78998042a: step3:optim 05/04 00:17
6F：→ a78998042a: 先不论理解正确与否，但下面写说，这个过程没问题， 05/04 00:18
7F：→ a78998042a: 有问题的是程式的输出结果，不过没有 真实资料 跟 05/04 00:19
8F：→ a78998042a: 假定model应该没有办法重现你的问题? 05/04 00:19
9F：→ a78998042a: 还是请直接提供整个程式?能重现错误也不会理解错误? 05/04 00:21
10F：→ 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
11F：推 a78998042a: 你好，谢谢补充，因为没有足够重现现象的条件，而过程 05/04 13:29
12F：→ a78998042a: 看起来没有错误，只能猜测一下。 05/04 13:29
13F：→ a78998042a: 怀疑是，你生成的过程是否没有setseed(set.seed(100)) 05/04 13:30
14F：→ a78998042a: 导致你重覆验证时，其实是用不同的估计值，所以有了 05/04 13:30
15F：→ a78998042a: 错误的误解?(会这样怀疑是因为上面写，NA""大约""有62 05/04 13:30
16F：→ a78998042a: )，是否先将bh做成一个matrix，将所有结果存下来，再 05/04 13:31
17F：→ 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 回圈的问题
非常感谢您热心的帮助，我看这一次跑不跑得出来，谢谢。
18F：→ x88776544pc: 如果 code 都没问题，即回圈中没有一次符合 if 条件 05/04 15:22
19F：→ x88776544pc: 你可以试试用 min 取代 <0.000001 去验证 05/04 15:27
您的意思是找与 0.5 相减的绝对值中，最小的那一个吗？
确实这样好像就不用控制误差，刚刚发现好像是误差太小了导致程式找不到，才会 NA
※ 编辑: strp (36.236.59.243), 05/04/2017 16:17:54
20F：→ a78998042a: seed是加在整个loop外部，所以每次生成的资料仍不同， 05/04 16:44
21F：→ 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