Я не сделал это точно так же, как вы, но я думаю, что ответ должен быть в порядке.
Настройка начальных условий
s0 <- list(b=0.1,c=0.1,d=0,Th=3)
X <- read.table("rogersdat.txt",header=TRUE)
Функция для вычисленного ожидаемого количества съеденного:
predfun <- function(b,c,d,Th,Te,N0,debug=FALSE) {
a <- (d+b*N0)/(1+c*N0)
r <- N0 - (1/(a*Th))*lambertW(a*Th*N0*exp(a*(Th*N0-Te)))
if (debug) cat(mean(a),b,c,d,Th,mean(r),"\n")
r
}
График данных (просто чтобы убедиться):
with(X,plot(Ne~N0))
## check starting value
lines(1:100,with(c(s0,X),predfun(b,c,d,Th,Te=72,N0=1:100)))
Fit:
library(emdbook)
n1 <- nls(Ne~predfun(b,c,d,Th,Te=72,N0),data=X,
lower=c(1e-6,1e-6,-Inf,1e-6),algorithm="port",
start=list(b=0.1,c=0.1,d=0.002,Th=3))
summary(n1)
Formula: Ne ~ predfun(b, c, d, Th, Te = 72, N0)
Parameters:
Estimate Std. Error t value Pr(>|t|)
b 0.0004155 0.0008745 0.475 0.63591
c 0.0000010 0.0657237 0.000 0.99999
d 0.0008318 0.0067374 0.123 0.90203
Th 4.0639937 1.4665102 2.771 0.00686 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.516 on 85 degrees of freedom
Algorithm "port", convergence message: relative convergence (4)
confint.default(n1)
2.5 % 97.5 %
b -0.001298523 0.002129547
c -0.128815083 0.128817083
d -0.012373153 0.014036799
Th 1.189686504 6.938300956
Оценки очень хорошо совпадают, доверительные интервалы не(эти типы доверительных интервалов на границе все же немного рискованны ...)
Я бы на самом деле предположил, что оценка максимального правдоподобия с биномиальными ошибками была бы немного лучше.
pdat <- data.frame(N0=1:100,Ne=predict(n1,newdata=data.frame(N0=1:100)))
with(pdat,lines(N0,Ne,col=2))
library(ggplot2)
ggplot(X,aes(x=N0,y=Ne))+stat_sum(aes(size=factor(..n..)),alpha=0.5)+theme_bw()+
geom_line(data=pdat,colour="red")
