如何计算正在优化的参数的置信区间？

Question

我有一些温度数据，我已经优化了使用负对数似然的方程。

下面是数据：

structure(list(date = structure(c(16436, 16437, 16438, 16439, 
16440, 16441, 16442, 16443, 16444, 16445, 16446, 16447, 16448, 
16449, 16450, 16451, 16452, 16453, 16454, 16455, 16456, 16457, 
16458, 16459, 16460, 16461, 16462, 16463, 16464, 16465, 16466, 
16467, 16468, 16469, 16470, 16471, 16472, 16473, 16474, 16475, 
16476, 16477, 16478, 16479, 16480, 16481, 16482, 16483, 16484, 
16485, 16486, 16487, 16488, 16489, 16490, 16491, 16492, 16493, 
16494, 16495, 16496, 16497, 16498, 16499, 16500, 16501, 16502, 
16503, 16504, 16505, 16506, 16507, 16508, 16509, 16510, 16511, 
16512, 16513, 16514, 16515, 16516, 16517, 16518, 16519, 16520, 
16521, 16522, 16523, 16524, 16525, 16526, 16527, 16528, 16529, 
16530, 16531, 16532, 16533, 16534, 16535, 16536, 16537, 16538, 
16539, 16540, 16541, 16542, 16543, 16544, 16545, 16546, 16547, 
16548, 16549, 16550, 16551, 16552, 16553, 16554, 16555, 16556, 
16557, 16558, 16559, 16560, 16561, 16562, 16563, 16564, 16565, 
16566, 16567, 16568, 16569, 16570, 16571, 16572, 16573, 16574, 
16575, 16576, 16577, 16578, 16579, 16580, 16581, 16582, 16583, 
16584, 16585, 16586, 16587, 16588, 16589, 16590, 16591, 16592, 
16593, 16594, 16595, 16596, 16597, 16598, 16599, 16600, 16601, 
16602, 16603, 16604, 16605, 16606, 16607, 16608, 16609, 16610, 
16611, 16612, 16613, 16614, 16615, 16616, 16617, 16618, 16619, 
16620, 16621, 16622, 16623, 16624, 16625, 16626, 16627, 16628, 
16629, 16630, 16631, 16632, 16633, 16634, 16635, 16636, 16637, 
16638, 16639, 16640, 16641, 16642, 16643, 16644, 16645, 16646, 
16647, 16648, 16649, 16650, 16651, 16652, 16653, 16654, 16655, 
16656, 16657, 16658, 16659, 16660, 16661, 16662, 16663, 16664, 
16665, 16666, 16667, 16668, 16669, 16670, 16671, 16672, 16673, 
16674, 16675, 16676, 16677, 16678, 16679, 16680, 16681, 16682, 
16683, 16684, 16685, 16686, 16687, 16688, 16689, 16690, 16691, 
16692, 16693, 16694, 16695, 16696, 16697, 16698, 16699, 16700, 
16701, 16702, 16703, 16704, 16705, 16706, 16707, 16708, 16709, 
16710, 16711, 16712, 16713, 16714, 16715, 16716, 16717, 16718, 
16719, 16720, 16721, 16722, 16723, 16724, 16725, 16726, 16727, 
16728, 16729, 16730, 16731, 16732, 16733, 16734, 16735, 16736, 
16737, 16738, 16739, 16740, 16741, 16742, 16743, 16744, 16745, 
16746, 16747, 16748, 16749, 16750, 16751, 16752, 16753, 16754, 
16755, 16756, 16757, 16758, 16759, 16760, 16761, 16762, 16763, 
16764, 16765, 16766, 16767, 16768, 16769, 16770, 16771, 16772, 
16773, 16774, 16775, 16776, 16777, 16778, 16779, 16780, 16781, 
16782, 16783, 16784, 16785, 16786, 16787, 16788, 16789, 16790, 
16791, 16792, 16793, 16794, 16795, 16796, 16797, 16798, 16799, 
16800), class = "Date"), mint = c(-15.9, -16.5, -14.4, -11.2, 
-5.7, -2.4, -2.5, -3.2, -4.3, -4.6, -1.5, -1, -0.9, -6.3, -7, 
-5.7, -1.2, -0.9, 0.3, -2.7, -5.9, -10.1, -8.7, -7.3, -5.7, -3.5, 
-1.2, -0.4, -0.9, -0.7, -4.3, -4.3, -2.8, 1, 2.7, 3.1, 5.8, 6.2, 
3.8, 2.2, -0.7, -1.5, -0.9, -0.3, 1, 1, -1.6, -3.8, -3.9, -1.9, 
-0.6, -0.8, -3.8, -7, -8.8, -7, -2.2, -0.3, -1.1, -2.9, -5.1, 
-5.2, -9.2, -9.7, -6.9, -4.2, -3.1, -3.5, -3.8, -2.3, 3.5, 0.3, 
0.7, 5.8, 7, 7.4, 2.3, -0.6, -2.2, 0.7, 0.9, 1.6, 3.8, -0.9, 
-2.5, 1, 2.6, 1.8, -1.6, 2.3, -4.2, -6.6, -4.7, -4.2, -0.5, -1.4, 
-3, 0.3, -2.9, -2.3, 1.1, -0.4, -1.5, 0.5, -6.1, -7.3, -5, -0.5, 
0.6, 0.7, 1.2, 2.9, 4.3, 4.7, 2.1, 0.3, 0.5, 1.4, 3.4, 5, 4.9, 
4.2, 6.3, 6.7, 6, 6.3, 3.6, 3.5, 3.7, 1.1, 1.9, 4.9, 0.7, 1.2, 
5.8, 5.6, 4, 6.2, 8.3, 7, 6, 4.7, 7, 9.2, 8.1, 6.9, 7.9, 8.6, 
9.6, 9.4, 10.3, 10.4, 9.6, 8.2, 9.4, 9.8, 7.2, 9.4, 10.8, 12.4, 
14.5, 11.8, 11, 10.7, 11.3, 10.8, 9.7, 10.4, 10.6, 12.1, 10.3, 
10.5, 11.3, 10, 12.6, 13.6, 17.4, 19.9, 19.9, 18.9, 18.4, 18.9, 
20.1, 19, 17, 16.9, 14.8, 13.1, 14, 11.5, 10.6, 11.1, 12.7, 11.4, 
11.9, 12.5, 13.3, 13.6, 13.2, 11.8, 11.8, 12.6, 15, 11.4, 10, 
9.6, 9.3, 9.3, 8.2, 9.6, 9.7, 12, 14.3, 16.1, 16.5, 12.8, 13.7, 
11.3, 10.3, 12.2, 11.4, 11.8, 11.1, 10.9, 11.2, 13, 11.8, 9, 
9.7, 8.9, 10.1, 10, 11.5, 10.6, 12.2, 10.9, 12.6, 11.9, 11.9, 
13.1, 13.4, 11.4, 6.9, 6, 7.7, 9.7, 7.8, 2.2, 1.5, 0.9, 2.3, 
4.8, 6.3, 8.3, 10.4, 11.2, 12.8, 11, 7.5, 6.1, 5.5, 2.4, 3.5, 
5.8, 5.9, 6.2, 5.6, 6.1, 7.4, 9.9, 7.8, 6.4, 7.8, 11, 10.1, 4.8, 
3.5, 6.6, 4.6, 5.5, 5.9, 9.8, 8.3, 8.6, 6.4, 4.4, 6, 7.1, 6.9, 
7.5, 7.8, 6.9, 3.9, 1.8, 0.3, 0.3, -0.5, 3.2, 2.4, -0.3, 0.2, 
5.1, -1.5, -1.4, 4.7, 5.6, 1.6, -1.3, -3.8, -4.1, -4.6, -3.5, 
-0.8, -1.4, -6.5, -6, -5, -4.9, -3.9, -4.2, -6.1, -1.7, -0.2, 
-0.3, -3.6, -7.1, -6.4, -3.4, -5.2, -8.6, -9.6, -13.8, -16.3, 
-15.6, -14.5, -11.8, -4.6, 0, -7.6, -7.7, -1.3, 4.8, 4.6, 2.3, 
0.1, -2.2, -1.4, -2.6, -4.7, -9, -6.8, -4.4, -3.7, -3.9, -5.1, 
0, -1.8, -3.2, -9, -14.2, -17.4, -13, -8.2, -12.7, -17.5), day = c(1, 
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 
52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 
68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 
84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 
100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 
113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 
126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 
139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 
152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 
165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 
178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 
191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 
204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 
217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 
230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 
243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 
256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 
269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 
282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 
295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 
308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 
321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 
334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 
347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 
360, 361, 362, 363, 364, 365)), .Names = c("date", "mint", "day" 
), row.names = 8768:9132, class = "data.frame") 

这里是我的优化：

TMIN <- function(a,b,c,x){a*sin(b*x+c)} 

plot(test$mint~test$day) 

normTMIN<-function(params,k,x){ 
    a=params[1] 
    b=params[2] 
    c=params[3] 
    d=params[4] 
    Mean<-TMIN(a,b,c,x) 
    -sum(dnorm(k,mean=Mean,sd=d,log=TRUE)) #shape= Mean(a,b)/scale 
} 

# qnorm(0.95,mean=TMIN(a,b,c,c(1:400)),sd=optTMIN$par['d']) 
# chi test 

optTMIN <- optim(par=c(a=20,b=0.017,c=11,d=1),k=test$mint,x=test$day,fn=normTMIN, control=list(trace = TRUE)) 

plot(test$mint~test$day) 
curve(TMIN(optTMIN$par[1],optTMIN$par[2],optTMIN$par[3],x),add=TRUE,col="blue") #,add=TRUE 

我的目标是建立一个模型，然后我就可以从模拟数据。我想知道如何获得optim估计值的置信区间。

Answer 1

假设你的目标函数是正确的（有疑惑），你可以引导残差这样的：

test$fitted <- TMIN(optTMIN$par[1],optTMIN$par[2],optTMIN$par[3],test$day) 
test$residuals <- test$mint - TMIN(optTMIN$par[1],optTMIN$par[2],optTMIN$par[3],test$day) 

library(boot) 
set.seed(42) 
bootres <- boot(test, function(DF, i) { 
    newy <- test$fitted + test[i, "residuals"] 
    optTMIN <- optim(par=c(a=20,b=0.017,c=11,d=1),k=newy,x=test$day,fn=normTMIN) 
    optTMIN$par 
}, R = 1000) 

#Bootstrap Statistics : 
#  original  bias std. error 
#t1* 16.7813787 -1.1483371368 1.639690649 
#t2* 0.0129091 -0.0006805658 0.007442933 
#t3* 11.6789783 -0.3390739337 0.332904632 
#t4* 6.5193201 -1.1309324860 0.520732955 

apply(bootres$t, 2, quantile, probs = c(0.025, 0.5, 0.975)) 
#   [,1]  [,2]  [,3]  [,4] 
#2.5% 13.47923 -0.01408214 10.46399 4.715159 
#50% 15.32346 0.01414520 11.43496 5.305173 
#97.5% 19.58337 0.01589570 11.57207 6.470475