R-Levenberg Marquardt中的非线性最小二乘拟合Heligman Pollard模型参数

我试图重现Kostakis的纸张解决方案。在本文中，使用de Heligman-Pollard模型将简略死亡率表扩展为完整的寿命表。该模型有8个参数必须安装。作者使用修改的高斯 - 牛顿算法;这个算法（E04FDF）是NAG计算机程序库的一部分。 Levenberg Marquardt不应该产生相同的一组参数吗？我的代码或LM算法的应用程序有什么问题？R-Levenberg Marquardt中的非线性最小二乘拟合Heligman Pollard模型参数

library(minpack.lm) 


## Heligman-Pollard is used to expand an abridged table. 
## nonlinear least squares algorithm is used to fit the parameters on nqx observed over 5 year intervals (5qx) 
AGE <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70) 
MORTALITY <- c(0.010384069, 0.001469140, 0.001309318, 0.003814265, 0.005378395, 0.005985625,  0.006741766, 0.009325056, 0.014149626, 0.021601755, 0.034271934, 0.053836246, 0.085287751, 0.136549522, 0.215953304) 

## The start parameters for de Heligman-Pollard Formula (Converged set a=0.0005893,b=0.0043836,c=0.0828424,d=0.000706,e=9.927863,f=22.197312,g=0.00004948,h=1.10003) 
## I modified a random parameter "a" in order to have a start values. The converged set is listed above. 
parStart <- list(a=0.0008893,b=0.0043836,c=0.0828424,d=0.000706,e=9.927863,f=22.197312,g=0.00004948,h=1.10003) 

## The Heligman-Pollard Formula (HP8) = qx/px = ...8 parameter equation 
HP8 <-function(parS,x) 
ifelse(x==0, parS$a^((x+parS$b)^parS$c) + parS$g*parS$h^x, 
      parS$a^((x+parS$b)^parS$c) + parS$d*exp(-parS$e*(log(x/parS$f))^2) + 
       parS$g*parS$h^x) 

## Define qx = HP8/(1+HP8) 
qxPred <- function(parS,x) HP8(parS,x)/(1+HP8(parS,x)) 

## Calculate nqx predicted by HP8 model (nqxPred(parStart,x)) 
nqxPred <- function(parS,x) 
(1 -(1-qxPred(parS,x)) * (1-qxPred(parS,x+1)) * 
    (1-qxPred(parS,x+2)) * (1-qxPred(parS,x+3)) * 
    (1-qxPred(parS,x+4))) 

##Define Residual Function, the relative squared distance is minimized 
ResidFun <- function(parS, Observed,x) (nqxPred(parS,x)/Observed-1)^2 

## Applying the nls.lm algo. 
nls.out <- nls.lm(par=parStart, fn = ResidFun, Observed = MORTALITY, x = AGE, 
        control = nls.lm.control(nprint=1, 
              ftol = .Machine$double.eps, 
              ptol = .Machine$double.eps, 
              maxfev=10000, maxiter = 500)) 

summary(nls.out) 


## The author used a modified Gauss-Newton algorithm, this alogorithm (E04FDF) is part of the NAG library of computer programs 
## Should not Levenberg Marquardt yield the same set of parameters

来源

2013-07-28 Arthur Rose

回车是您的朋友。 –

@HongOoi，不再。 –

“有四个参数，我可以适合一头大象，五头，我可以让他摆动他的后备箱。”（[John von Neumann]（http://en.wikiquote.org/wiki/John_von_Neumann））我相信这个是过度配合的严重情况。可能有许多当地的最低等级和其他nasties。制作一些诊断图来检查参数灵敏度。如果你有这样的问题，不同的算法会给出不同的结果。顺便说一下，你为什么不使用'nlsLM'前端？ – Roland

这里的底线是，@Roland是绝对正确的，这是一个非常病态问题，你不应该一定期望得到可靠的答案。下面，我

清理在一些小的方面的代码（这仅仅是审美）
改变了ResidFun返回残差，不残差平方。（前者是正确的，但这并没有太大区别。）
探索了几种不同优化器的结果。它实际上看起来像你得到的答案是更好比你上面列出的“融合参数”，我假设是来自原研究的参数（你能提供一个参考吗？）。

负载包：

library(minpack.lm)

数据，作为数据帧：

d <- data.frame(
    AGE = seq(0,70,by=5), 
    MORTALITY=c(0.010384069, 0.001469140, 0.001309318, 0.003814265, 
       0.005378395, 0.005985625, 0.006741766, 0.009325056, 
       0.014149626, 0.021601755, 0.034271934, 0.053836246, 
       0.085287751, 0.136549522, 0.215953304))

的数据的第一视图：

library(ggplot2) 
(g1 <- ggplot(d,aes(AGE,MORTALITY))+geom_point()) 
g1+geom_smooth() ## with loess fit

参数选择：

想必这些都是从原来的文件中的参数...

parConv <- c(a=0.0005893,b=0.0043836,c=0.0828424, 
      d=0.000706,e=9.927863,f=22.197312,g=0.00004948,h=1.10003)

摄动参数：

parStart <- parConv 
parStart["a"] <- parStart["a"]+3e-4

的公式：

HP8 <-function(parS,x) 
    with(as.list(parS), 
     ifelse(x==0, a^((x+b)^c) + g*h^x, 
       a^((x+b)^c) + d*exp(-e*(log(x/f))^2) + g*h^x)) 
## Define qx = HP8/(1+HP8) 
qxPred <- function(parS,x) { 
    h <- HP8(parS,x) 
    h/(1+h) 
} 
## Calculate nqx predicted by HP8 model (nqxPred(parStart,x)) 
nqxPred <- function(parS,x) 
    (1 -(1-qxPred(parS,x)) * (1-qxPred(parS,x+1)) * 
    (1-qxPred(parS,x+2)) * (1-qxPred(parS,x+3)) * 
    (1-qxPred(parS,x+4))) 
##Define Residual Function, the relative squared distance is minimized 
ResidFun <- function(parS, Observed,x) (nqxPred(parS,x)/Observed-1)

注：这从OP的版本略有改变; nls.lm想要残差，而不是平方残差。

求和的平方功能与其他优化利用：

ssqfun <- function(parS, Observed, x) { 
    sum(ResidFun(parS, Observed, x)^2) 
}

应用nls.lm。（不知道为什么ftol和ptol降低了从sqrt(.Machine$double.eps)到.Machine$double.eps - 前一般是实际极限精度...

nls.out <- nls.lm(par=parStart, fn = ResidFun, 
        Observed = d$MORTALITY, x = d$AGE, 
        control = nls.lm.control(nprint=0, 
              ftol = .Machine$double.eps, 
              ptol = .Machine$double.eps, 
              maxfev=10000, maxiter = 1000)) 

parNLS <- coef(nls.out) 

pred0 <- nqxPred(as.list(parConv),d$AGE) 
pred1 <- nqxPred(as.list(parNLS),d$AGE) 

dPred <- with(d,rbind(data.frame(AGE,MORTALITY=pred0,w="conv"), 
       data.frame(AGE,MORTALITY=pred1,w="nls"))) 

g1 + geom_line(data=dPred,aes(colour=w))

的线是不可区分的，但这些参数有一些大差异：

round(cbind(parNLS,parConv),5) 
##  parNLS parConv 
## a 1.00000 0.00059 
## b 50.46708 0.00438 
## c 3.56799 0.08284 
## d 0.00072 0.00071 
## e 6.05200 9.92786 
## f 21.82347 22.19731 
## g 0.00005 0.00005 
## h 1.10026 1.10003

d，F，G，H是接近，但A，B，C是不同数量级和e是50％不同。

综观原方程，这里发生了什么是a^((x+b)^c)是越来越设置为常数，因为a则接近1：一旦a约为1，b和c基本上是无关紧要的。

让我们来看看相关（我们需要一个广义逆，因为矩阵是如此强烈的相关性）：

obj <- nls.out 
vcov <- with(obj,deviance/(length(fvec) - length(par)) * 
       MASS::ginv(hessian)) 

cmat <- round(cov2cor(vcov),1) 
dimnames(cmat) <- list(letters[1:8],letters[1:8]) 

##  a b c d e f g h 
## a 1.0 0.0 0.0 0.0 0.0 0.0 -0.1 0.0 
## b 0.0 1.0 -1.0 1.0 -1.0 -1.0 -0.4 -1.0 
## c 0.0 -1.0 1.0 -1.0 1.0 1.0 0.4 1.0 
## d 0.0 1.0 -1.0 1.0 -1.0 -1.0 -0.4 -1.0 
## e 0.0 -1.0 1.0 -1.0 1.0 1.0 0.4 1.0 
## f 0.0 -1.0 1.0 -1.0 1.0 1.0 0.4 1.0 
## g -0.1 -0.4 0.4 -0.4 0.4 0.4 1.0 0.4 
## h 0.0 -1.0 1.0 -1.0 1.0 1.0 0.4 1.0

这实际上不是那么有用 - 它真的只是确认变量的地段是强烈的相关性...

library(optimx) 
mvec <- c('Nelder-Mead','BFGS','CG','L-BFGS-B', 
      'nlm','nlminb','spg','ucminf') 
opt1 <- optimx(par=parStart, fn = ssqfun, 
     Observed = d$MORTALITY, x = d$AGE, 
       itnmax=5000, 
       method=mvec,control=list(kkt=TRUE)) 
       ## control=list(all.methods=TRUE,kkt=TRUE)) ## Boom! 

##   fvalues  method fns grs itns conv KKT1 KKT2 xtimes 
## 2 8.988466e+307  BFGS NA NULL NULL 9999 NA NA  0 
## 3 8.988466e+307   CG NA NULL NULL 9999 NA NA  0 
## 4 8.988466e+307 L-BFGS-B NA NULL NULL 9999 NA NA  0 
## 5 8.988466e+307   nlm NA NA NA 9999 NA NA  0 
## 7  0.3400858   spg 1 NA 1 3 NA NA 0.064 
## 8  0.3400858  ucminf 1 1 NULL 0 NA NA 0.032 
## 1 0.06099295 Nelder-Mead 501 NA NULL 1 NA NA 0.252 
## 6 0.009275733  nlminb 200 1204 145 1 NA NA 0.708

此发出警告坏缩放，同时还发现了各种不同的答案：只有ucminf号称有converg ED，但nlminb得到了更好的答案 - 和itnmax参数似乎被忽略......

opt2 <- nlminb(start=parStart, objective = ssqfun, 
     Observed = d$MORTALITY, x = d$AGE,     
       control= list(eval.max=5000,iter.max=5000)) 

parNLM <- opt2$par

饰面，而是用假收敛警告...

round(cbind(parNLS,parConv,parNLM),5) 

##  parNLS parConv parNLM 
## a 1.00000 0.00059 1.00000 
## b 50.46708 0.00438 55.37270 
## c 3.56799 0.08284 3.89162 
## d 0.00072 0.00071 0.00072 
## e 6.05200 9.92786 6.04416 
## f 21.82347 22.19731 21.82292 
## g 0.00005 0.00005 0.00005 
## h 1.10026 1.10003 1.10026 

sapply(list(parNLS,parConv,parNLM), 
     ssqfun,Observed=d$MORTALITY,x=d$AGE) 
## [1] 0.006346250 0.049972367 0.006315034

它看起来像nlminb和minpack.lm得到了类似的答案，并且实际上做的是更好的比原来陈述的参数（相当多）：

pred2 <- nqxPred(as.list(parNLM),d$AGE) 

dPred <- with(d,rbind(dPred, 
       data.frame(AGE,MORTALITY=pred2,w="nlminb"))) 

g1 + geom_line(data=dPred,aes(colour=w)) 
ggsave("cmpplot.png")

enter image description here

ggplot(data=dPred,aes(x=AGE,y=MORTALITY-d$MORTALITY,colour=w))+ 
    geom_line()+geom_point(aes(shape=w),alpha=0.3) 
ggsave("residplot.png")

enter image description here

其他的事情一个可以尝试将是：

适当的比例 - 虽然这种快速测试似乎并没有帮助那么多
提供分析梯度
使用AD模型构建器
使用slice功能从bbmle探索新旧参数是否似乎代表不同极小的旧参数，或者是否只是虚假的收敛...
得到KKT（Karsh-从optimx或相关的包类似的检查

PS工作库恩 - 塔克）标准计算器：最大偏差（迄今为止）是最古老的年龄组，这可能也有小样本。从统计的角度来看，它可能是值得做一个合适的加权的个人点的精度...

来源

2013-07-28 15:52:42

我印象深刻。应该提到的是，如果您尝试适应如此多的参数，您应该拥有更多的数据。最好用独立的附加实验来估计一些参数。优选的结果应该用独立的数据进行验证，或者至少进行交叉验证。 – Roland

@BenBolker，谢谢你的回复。由于我无法附上论文，我通过电子邮件向您发送了论文。我从“官方统计杂志，第7卷，第3号，1991年，第311-323页”得到这篇论文。标题： –

@BenBolker标题：Heligman-Pollard公式作为扩展缩写生命表的工具作者：阿纳斯塔西娅科斯塔基链接：http://www.jos.nu/Articles/abstract.asp?article=73311 –

@BenBolker，拟合参数与整个数据集（基础qx）值。仍无法重现参数

library(minpack.lm) 

library(ggplot2) 

library(optimx) 

getwd() 

d <- data.frame(AGE = seq(0,74), MORTALITY=c(869,58,40,37,36,35,32,28,29,23,24,22,24,28, 
              33,52,57,77,93,103,103,109,105,114,108,112,119, 
              125,117,127,125,134,134,131,152,179,173,182,199, 
              203,232,245,296,315,335,356,405,438,445,535,594, 
              623,693,749,816,915,994,1128,1172,1294,1473, 
              1544,1721,1967,2129,2331,2559,2901,3203,3470, 
              3782,4348,4714,5245,5646)) 


d$MORTALITY <- d$MORTALITY/100000 

ggplot(d,aes(AGE,MORTALITY))+geom_point() 

##Not allowed to post Images 

g1 <- ggplot(d,aes(AGE,MORTALITY))+geom_point() 

g1+geom_smooth()## with loess fit

报告的参数：

parConv <- c(a=0.0005893,b=0.0043836,c=0.0828424,d=0.000706,e=9.927863,f=22.197312, 
      g=0.00004948,h=1.10003) 

parStart <- parConv 

parStart["a"] <- parStart["a"]+3e-4 


## Define qx = HP8/(1+HP8) 

HP8 <-function(parS,x) 
with(as.list(parS), 
ifelse(x==0, a^((x+b)^c) + g*h^x, a^((x+b)^c) + d*exp(-e*(log(x/f))^2) + g*h^x)) 



qxPred <- function(parS,x) { 
    h <- HP8(parS,x) 
    h/(1+h) 
} 



##Define Residual Function, the relative squared distance is minimized, 
ResidFun <- function(parS, Observed,x) (qxPred(parS,x)/Observed-1) 

ssqfun <- function(parS, Observed, x) { 
    sum(ResidFun(parS, Observed, x)^2) 
} 

nls.out <- nls.lm(par=parStart, fn = ResidFun, Observed = d$MORTALITY, x = d$AGE, 
        control = nls.lm.control(nprint=1, ftol = sqrt(.Machine$double.eps), 
        ptol = sqrt(.Machine$double.eps), maxfev=1000, maxiter=1000)) 


parNLS <- coef(nls.out) 

pred0 <- qxPred(as.list(parConv),d$AGE) 
pred1 <- qxPred(as.list(parNLS),d$AGE) 


#Binds Row wise the dataframes from pred0 and pred1 
dPred <- with(d,rbind(data.frame(AGE,MORTALITY=pred0,w="conv"), 
     data.frame(AGE,MORTALITY=pred1,w="nls"))) 


g1 + geom_line(data=dPred,aes(colour=w)) 

round(cbind(parNLS,parConv),7) 

mvec <- c('Nelder-Mead','BFGS','CG','L-BFGS-B','nlm','nlminb','spg','ucminf') 
opt1 <- optimx(par=parStart, fn = ssqfun, 
    Observed = d$MORTALITY, x = d$AGE, 
    itnmax=5000, 
    method=mvec, control=list(all.methods=TRUE,kkt=TRUE,) 
## control=list(all.methods=TRUE,kkt=TRUE)) ## Boom 

get.result(opt1, attribute= c("fvalues","method", "grs", "itns", 
      "conv", "KKT1", "KKT2", "xtimes")) 

##  method  fvalues grs itns conv KKT1 KKT2 xtimes 
##5   nlm 8.988466e+307 NA NA 9999 NA NA  0 
##4 L-BFGS-B 8.988466e+307 NULL NULL 9999 NA NA  0 
##2   CG 8.988466e+307 NULL NULL 9999 NA NA 0.02 
##1  BFGS 8.988466e+307 NULL NULL 9999 NA NA  0 
##3 Nelder-Mead  0.5673864 NA NULL 0 NA NA 0.42 
##6  nlminb  0.4127198 546 62 0 NA NA 0.17 


opt2 <- nlminb(start=parStart, objective = ssqfun, 
    Observed = d$MORTALITY, x = d$AGE, 
    control= list(eval.max=5000,iter.max=5000)) 

parNLM <- opt2$par

检查参数：

round(cbind(parNLS,parConv,parNLM),5) 

## parNLS parConv parNLM 
##a 0.00058 0.00059 0.00058 
##b 0.00369 0.00438 0.00369 
##c 0.08065 0.08284 0.08065 
##d 0.00070 0.00071 0.00070 
##e 9.30948 9.92786 9.30970 
##f 22.30769 22.19731 22.30769 
##g 0.00005 0.00005 0.00005 
##h 1.10084 1.10003 1.10084

SSE评论：

sapply(list(parNLS,parConv,parNLM), 
    ssqfun,Observed=d$MORTALITY,x=d$AGE) 

##[1] 0.4127198 0.4169513 0.4127198

无法上传的绘图，但代码是这里。当使用完整的死亡率数据（未删减或子集）时，仍然显示文章中找到的参数不是最合适的。

##pred2 <- qxPred(as.list(parNLM),d$AGE) 

##dPred <- with(d,rbind(dPred, 
    data.frame(AGE,MORTALITY=pred2,w="nlminb"))) 

##g1 + geom_line(data=dPred,aes(colour=w)) 
ggplot(data=dPred,aes(x=AGE,y=MORTALITY-d$MORTALITY,colour=w)) 
     + geom_line()+geom_point(aes(shape=w),alpha=0.3)

来源

2013-07-30 13:13:56

@BenStolker，我提供给你int他回答完整的数据集。 –

R-Levenberg Marquardt中的非线性最小二乘拟合Heligman Pollard模型参数

回答

相关问题