具有优化参数的函数与在 R 中使用 mle2 的数据不接近

Question

因此，我一直在尝试使用伽马误差分布优化 Michaelis-Menten 关系，以对我收集的一些数据的平均值进行建模。但是，无论我如何优化函数，我得到的最低 AIC 是用于甚至不接近数据的参数。有什么办法可以解决吗？

这是我的代码：

我首先创建一个最大似然函数：

MicNLL <- function(a,b){
  #a=150.6727
  #b=319.7007 optim val
  top <- a*x
  bot <- b+x
  Mic <- top/bot
  nll <- -sum(dgamma(y, shape=(Mic^2/var(x)), scale=(var(x)/Mic), log=TRUE))
  return(nll)
}

然后我使用bbmle包中的mle2()函数编写了优化函数：

MN <- mle2(minuslogl = MicNLL, parameters=list(a~Treatment, b~Treatment), start=list(a=100,b=260), data=list(x=NSug3$VolpulT, y=NSug3$SugarpugT), control=list(maxit=1e4), method="SANN", hessian=T)

MN 
AICMN <- (2*2)-(2*logLik(MN))
AICMN

虽然 a=100 和 b=260 的目测参数很适合我的数据，但它通常会将参数优化为 a=242 和 b=182，结果

Michealis <- function(a, b, x){
  top <- a*x
  bot <- b+x
  Mic <- top/bot
  return(Mic)
}
ggplot(NSug3, aes(x=VolpulT, y=SugarpugT))+
  geom_point(stat="identity", size=0.8)+
  theme_classic()+
  ggtitle("help")+
  ylab("Sugar concentration")+
  xlab("Volume per Extra floral nectary")+
  stat_function(fun= Michealis, args=c(a=100, b=260), colour="Orange", size=0.725)+
  stat_function(fun= Michealis, args=c(a=MN@coef[[1]], b=MN@coef[[2]]), colour="Red", size=0.725)

长话短说，我怎样才能确保我的优化模型真正贯穿我的数据？

Answer 1

对下面的脑残代码表示歉意...

我做了一个与你类似的可重复的例子，似乎给出了合理的结果。

• 您是否收到任何关于收敛失败/“达到最大迭代次数”的警告？
• 您的代码中似乎有一些关于处理的未使用/剩余内容；这是个好主意，但仅适用于公式界面（见下文）

一些辅助函数：

## Gamma parameterized by mean and variance
## m = a*s, v = a*s^2 -> s=v/m; a=m^2/v
rgamma2 <- function(n, m, v) {
    rgamma(n, shape=m^2/v, scale=v/m)
}
dgamma2 <- function(x, m, v, log=FALSE) {
    dgamma(x, shape=m^2/v, scale=v/m, log=log)
}
sgamma2<- function(m, v) {   ## for predict()
    list(title="Gamma", mean=m, sd=sqrt(v))
}
mm <- function(x, a=100, b=260) {
    a*x/(b+x)
}

模拟数据：

set.seed(101)
x <- rlnorm(100,meanlog=4,sdlog=1)
dd <- data.frame(x,y=rgamma2(100,m=mm(x), v= 100))

拟合（使用公式界面）：

library(bbmle)
m1 <-mle2(y~dgamma2(m=mm(x,a,b),v=exp(logv)),
     start=list(a=50,b=200,logv=0),
     data=dd,
     control=list(maxit=1000))

绘制结果：

plot(y~x,data=dd)
lines(sort(dd$x),mm(sort(dd$x)),col=2)     ## true
lines(sort(dd$x),sort(predict(m1)),col=3)  ## predicted