编写一个程序来最小化递归指数函数的平方和

Question

这是我想在 R 中编写的函数，

i = 1,2,3,....j-1

a,b,c,f,g 由 nls 确定（起始值任意设置为 7,30,15,1,2）

S 和 Y 在数据集中

函数可以用计算更友好的递归方程来表示，

这是我对代码的尝试，但我无法让它收敛，

S=c(235,90,1775,960,965,1110,370,485,667,140,588,10,0,1340,600,0,930,1250,930,120,895,825,0,935,695,270,0,610,0,0,445,0,0,370,470,819,717,0,0,60,0,135,690,0,825,730,1250,370,1010,261,0,865,570,1425,150,1515,1143,0,675,1465,375,0,690,290,0,430,735,510,270,450,1044,0,928,60,95,105,60,950,0,1640,3960,1510,500,1135,0,0,0,181,568,60,1575,247,0,1270,870,290,510,0,540,455,120,580,420,90,525,1116,499,0,60,150,660,1080,1715,90,1090,840,975,280,850,633,30,1530,1765,880,150,225,77,1380,810,835,0,540,1017,1108,0,300,600,90,370,910,0,60,60,0,0,0,0,50,0,735,900)

Y=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,7.7,NA,NA,7.2,NA,NA,NA,NA,NA,NA,7.4,NA,NA,NA,NA,NA,NA,10.7,NA,NA,NA,NA,8.1,8.5,NA,NA,NA,NA,NA,9.9,NA,7.4,NA,NA,NA,9.5,NA,NA,9,NA,NA,NA,8.8,NA,NA,8.5,NA,NA,NA,6.9,NA,NA,7.9,NA,NA,NA,7.3,NA,7.9,8.3,NA,NA,NA,11.5,NA,NA,12.3,NA,NA,NA,6.1,NA,NA,9,NA,NA,NA,10.3,NA,NA,9.7,NA,NA,8.6,NA,9.1,NA,NA,11,NA,NA,12.4,11.1,10.1,NA,NA,NA,NA,11.7,NA,NA,9,NA,NA,NA,10.2,NA,NA,11.2,NA,NA,NA,11.8,NA,9.2,10,9.8,NA,9.5,11.3,10.3,9.5,10.2,10.6,NA,10.8,10.7,11.1,NA,NA,NA,NA,NA,NA,NA,NA,12.6,NA)

mydata = data.frame(Y,S)

f <- function(a,b,f,c,g,m) {

    model <- matrix(NA,nrow(m)+1,3)

    model[1,1]=0
    model[1,2]=0
    model[1,3]=a

    for (i in 2:nrow(model)){
        model[i,1]=exp(-1/c)*model[i-1,1] + m$S[i-1] 
        model[i,2]=exp(-1/g)*model[i-1,2] + m$S[i-1]
        model[i,3]=a+b*model[i,1]-f*model[i,2]
    }
    model <- as.data.frame(model)
    colnames(model) = c('l','m','Y')
    model$Y[which(m$Y>0)]
}

Y=mydata$Y
nls(Y ~ f(a,b,f,c,g,mydata), start=list(a=7,b=5.3651,f=5.3656,c=16.50329,g=16.5006),control=list(maxiter=1000,minFactor=1e-12))

我得到的错误取决于起始值：

nls 中的错误(Y ~ f(a, b, f, c, g, mydata), start = list(a = 7, :
迭代次数超过最大值 1000

nls 中的错误(Y ~ f(a, b, f, c, g, mydata), start = list(a = 7, :
奇异梯度

我被卡住了，不知道该怎么做，任何帮助将不胜感激。

Answer 1

试试这个：

ff <- function(a,b,f,c,g) {
   Y <- numeric(length(S))
   for(i in seq(from=2, to=length(S))) {
      j <- seq(length=i-1)
      Y[i] <- a + sum((b*exp(-(i-j)/c) - f*exp(-(i-j)/g))*S[j])
   }
   Y
}

S <- c(235,90,1775,960,965,1110,370,485,667,140,588,10,0,1340,600,0,930,1250,930,120,895,825,0,935,695,270,0,610,0,0,445,0,0,370,470,819,717,0,0,60,0,135,690,0,825,730,1250,370,1010,261,0,865,570,1425,150,1515,1143,0,675,1465,375,0,690,290,0,430,735,510,270,450,1044,0,928,60,95,105,60,950,0,1640,3960,1510,500,1135,0,0,0,181,568,60,1575,247,0,1270,870,290,510,0,540,455,120,580,420,90,525,1116,499,0,60,150,660,1080,1715,90,1090,840,975,280,850,633,30,1530,1765,880,150,225,77,1380,810,835,0,540,1017,1108,0,300,600,90,370,910,0,60,60,0,0,0,0,50,0,735,900)
Y <-  c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,7.7,NA,NA,7.2,NA,NA,NA,NA,NA,NA,7.4,NA,NA,NA,NA,NA,NA,10.7,NA,NA,NA,NA,8.1,8.5,NA,NA,NA,NA,NA,9.9,NA,7.4,NA,NA,NA,9.5,NA,NA,9,NA,NA,NA,8.8,NA,NA,8.5,NA,NA,NA,6.9,NA,NA,7.9,NA,NA,NA,7.3,NA,7.9,8.3,NA,NA,NA,11.5,NA,NA,12.3,NA,NA,NA,6.1,NA,NA,9,NA,NA,NA,10.3,NA,NA,9.7,NA,NA,8.6,NA,9.1,NA,NA,11,NA,NA,12.4,11.1,10.1,NA,NA,NA,NA,11.7,NA,NA,9,NA,NA,NA,10.2,NA,NA,11.2,NA,NA,NA,11.8,NA,9.2,10,9.8,NA,9.5,11.3,10.3,9.5,10.2,10.6,NA,10.8,10.7,11.1,NA,NA,NA,NA,NA,NA,NA,NA,12.6,NA)
nls(Y ~ f(a,b,f,c,g,mydata), start=list(a=7,b=5.3651,f=5.3656,c=16.50329,g=16.5006))

但我无法让 nls 在这里运行。您也可以尝试通用优化器。构造平方和函数（我们最大化它时的平方和）：

SS <- function(par) {
   a <- par[1]
   b <- par[2]
   f <- par[3]
   c <- par[4]
   g <- par[5]
  -sum((Y - ff(a,b,f,c,g))^2, na.rm=TRUE)
}

并最大化：

library(maxLik)
summary(a <- maxBFGS(SS, start=start))

它有效，但正如您所见，渐变仍然很大。如果我对 BFGS 的输出值重新运行 NR 优化器，我的梯度会变小：

summary(b <- maxNR(SS, start=coef(a)))

给出结果

Newton-Raphson maximisation 
Number of iterations: 1 
Return code: 2 
successive function values within tolerance limit 
Function value: -47.36338 
Estimates:
   estimate      gradient
a 10.584488  0.0016371615
b  6.954444 -0.0043306656
f  6.955095  0.0043327901
c 28.622035 -0.0005735572
g 28.619185  0.0003871179

我不知道这是否有意义。 nls 和其他优化器的问题表明您存在数值不稳定性，这与大数值或模型公式中的指数差异有关。

检查那里发生了什么:-)