Использование mle2 для оценки параметров с ошибкой и прогнозированием - PullRequest
Новая
       108

Использование mle2 для оценки параметров с ошибкой и прогнозированием

0 голосов
/ 31 октября 2019

Я использую mle2 для оценки параметров нелинейной модели, и мне нужны оценки погрешности вокруг оценок параметров (стандартная ошибка). Кроме того, я хотел бы использовать модель для прогнозирования с использованием новых данных, и у меня возникают проблемы (ошибки) с несколькими шагами в этом процессе.

Вот данные: 

table<- "ac.temp performance
1      2.17   47.357923
4      2.17  234.255317
7      2.17  138.002633
10     2.17  227.545902
13     2.17   28.118072
16     9.95  175.638448
19     9.95  167.392218
22     9.95  118.162747
25     9.95  102.770622
28     9.95  191.874867
31    16.07  206.897159
34    16.07   74.741619
37    16.07  127.219884
40    16.07  208.231559
42    16.07   89.544612
44    20.14  314.946107
47    20.14  290.994063
50    20.14  243.322497
53    20.14  192.335133
56    20.14  133.841776
58    23.83  139.746673
61    23.83  224.135993
64    23.83  126.726493
67    23.83  246.443386
70    23.83  163.019896
83    28.04    4.614154
84    28.04    2.851866
85    28.04    2.935584
86    28.04  153.868415
87    28.04  103.884295
88    30.60    0.000000
89    29.60    0.000000
90    30.30    0.000000
91    29.90    0.000000
92    30.80    0.000000
93    28.90    0.000000
94    30.00    0.000000
95    30.20    0.000000
96    30.40    0.000000
97    30.70    0.000000
98    27.90    0.000000
99    28.60    0.000000
100   28.60    0.000000
101   30.40    0.000000
102   30.60    0.000000
103   29.70    0.000000
104   30.50    0.000000
105   29.30    0.000000
106   28.90    0.000000
107   30.40    0.000000
108   30.20    0.000000
109   30.10    0.000000
110   29.50    0.000000
111   31.00    0.000000
112   30.60    0.000000
113   30.90    0.000000
114   31.10    0.000000"

perfdat<- read.table(text=table, header=T)

Сначала я должен установить пару фиксированных параметров для моей нелинейной модели поведения животного по температуре 

pi = mean(subset(perfdat, ac.temp<5)$performance)
ti = min(perfdat$ac.temp)

определение переменной x (температура) 

t = perfdat$ac.temp

создание функции для нелинейной модели 

tpc = function(tm, topt, pmax) {
    perfmu = pi+(pmax-pi)*(1+((topt-t)/(topt-tm)))*(((t-ti)/(topt-ti))^((tm-ti)/(topt-tm)))
    perfmu[perfmu<0] = 0.00001
    return(perfmu)
}

создание функции отрицательного логарифмического правдоподобия 

LL1 = function(tm, topt, pmax, performance=perfdat$performance) {
    perfmu = tpc(tm=tm, topt=topt, pmax=pmax)
    loglike = -sum(dnorm(x=performance, mean=perfmu, log=TRUE))
    return(loglike)
}

производительность модели с использованием mle 2 - maximum likelihood estimator

m1<- mle2(minuslogl = LL1, start=list(tm=15, topt=20, pmax=150), control=list(maxit=5000))

summary(m1)

Это дает мне оценки параметров, но не оценки ошибок (стандартная ошибка) с warning message: In sqrt(diag(object@vcov)) : NaNs produced. Тем не менее, оценки параметров являются хорошими и дают мне прогнозы, которые имеют смысл.

График нелинейной кривой с использованием оценок параметров

Я пробовал много разных оптимизаторов и методови получить ту же ошибку о невозможности рассчитать стандартное отклонение. ошибка, обычно с предупреждением о невозможности инвертировать гессиан. ИЛИ я получаю действительно сомнительные оценки своих параметров, которые не имеют смысла.

ЕСЛИ я использую: 

confint(m1)

Я получаю 95% интервалы для каждого из моих параметров, но я не могу включить их в метод прогнозирования, который я мог бы использовать для создания графика, подобногониже, который я сделал, используя модель nls и predict():

нелинейную модель с отображенной ошибкой

Если я воссоздаю свою модель mle2() повстраивание формулы модели в функцию mle2: 

tpcfun<- function(t, tm.f, topt.f, pmax.f) {
    perfmu = pi+(pmax.f-pi)*(1+((topt.f-t)/(topt.f-tm)))*(((t-ti)/(topt.f-ti))^((tm.f-ti)/(topt.f-tm.f)))
    perfmu[perfmu<0] = 0.00001
    return(perfmu)
}

m2<- mle2(performance ~ dnorm(mean=-sum(log(tpcfun(t=ac.temp, tm.f, topt.f, pmax.f))), sd=1),  method='L-BFGS-B', lower=c(tm.f=1, topt.f=1, pmax.f=1), control=list(maxit=5000, parscale=c(tm.f=1e1, topt.f=1e1, pmax.f=1e2)), start=list(tm.f=15, topt.f=20, pmax.f=150), data=perfdat)

summary(m2)

Я получаю бессмысленные оценки для своих параметров и все еще не получаю оценки для ошибки.

Мой вопрос таков:Кто-нибудь может увидеть что-то не так с одной из моих моделей (функции модели и функции правдоподобия) или что-то еще, что я делаю неправильно? У меня есть ощущение, что я могу написать неправильную функцию вероятности, но я пробовал все виды дистрибутивов и разными способами, но, возможно, я все испортил.

ИЛИ есть ли способ получить оценки погрешности вокруг моих параметров, чтобы я мог использовать их для визуализации погрешности прогнозирования моей модели в графиках?

Спасибо,

Rylee

PS. Я решил построить график с точечными оценками, а затем с трендовой линией без ошибок вокруг него, но я хотел поставить столбцы для 95% -ного доверительного интервала на каждой из точечных оценок, но confint () дает мне повторно малые CI, которые не даютдаже не отображаются на графике, потому что они меньше, чем точечный символ, который я использую, ха.

...