r - 在R中使用nls重新创建研究

Question

我正在学习如何在R中使用nls函数，并且遇到了一些问题。我只是试图重新创建目前在研究论文中发现的曲线。该模型对1987年股市崩盘前的股市走势拟合曲线。

我已经定义了一个函数func，如下所示:

func <- function(a,b,tc,t){
 a+b*log(tc-t)
}

我这样称呼nls:

nls1 <- nls(Y ~ func(a,b,tc,t), data2, start=list(a=0, b=1, tc=1466, t=1))

data2是一个由两列组成的数据框，一列是日期，另一列是一个值。有1466行。

head(data2)
 Date      Y
1  1/4/82 882.52
2  1/5/82 865.30
3  1/6/82 861.02
4  1/7/82 861.78
5  1/8/82 866.53
6 1/11/82 850.46

运行nls时收到以下消息，

Error in qr(.swts * attr(rhs, "gradient")) : 
  dims [product 4] do not match the length of object [1466]

In addition: Warning message:

In .swts * attr(rhs, "gradient") :
  longer object length is not a multiple of shorter object length

据我所知，这是数据帧设置方式的问题，但我找不到解决方案。

知道我该如何搬家吗？

非常感谢您的帮助。

Answer 1

基本问题是您尚未指定自变量。通过为start(...)指定a, b, tc, and t，您将告诉nls(...)这些都是模型的所有参数。

看起来您正在使用LPPL模型的简化版本，其中a, b, and tc是参数，而t是自变量。看来data2$Date包含时间变量。您需要确保data2$Date属于POSIXct类。所以你可以这样写:

df$Date <- as.POSIXct(df$Date, format="%m/%d/%y")
nls1 <- nls(Y~a+b*log(tc-Date), data=data2, start=list(a=0, b=1, tc=1466))

编辑:回应OP的评论

这是一个很好的问题，因为它说明了使用nls(...)的几个问题。您遇到的问题(现在已经正确指定了模型)是nls(...)没有收敛-令人痛苦的常见情况。基本上，除非您的起始参数估计值相对接近最终的拟合值(或者除非模型表现得非常“良好”)，否则nls将会失败。 [还请注意，您引用的论文提到b限制为b <0，而您以b = 1开头。]那么该怎么办？

软件包minpack.lm(...)中的minpack函数使用异常健壮的Levenberg-Marquardt算法进行非线性最小二乘估计。实际上，您引用的论文特别提到了L-M。 minpack.lm(...)的问题在于使用起来更加困难(您必须定义一个函数，该函数将在给定步骤返回残差，而不仅仅是定义适合的函数)。另外，minpack.lm(...)不会计算拟合度的统计信息。

因此解决方案是同时使用它们!使用minpack.lm(...)估计参数，然后将其用作nls(...)中的“起始值”。下面的代码可以做到这一点。使用nls(...)拟合模型将使生成拟合度，预测值，残差的统计信息以及将模型应用于新数据集变得更加容易。

# this section just grabs the DJIA for 1982 - 1987; you already have this
library(tseries)
library(zoo)
ts <- get.hist.quote(instrument="DJIA", 
                     start="1982-01-01", end="1987-08-01", 
                     quote="Close", provider="yahoo", origin="1970-01-01",
                     compression="d", retclass="zoo")
df <- data.frame(ts)
df <- data.frame(Date=as.Date(rownames(df)),Y=df$Close)
df <- df[!is.na(df$Y),]
# end of setup...
library(minpack.lm) # for nls.lm(...)
library(ggplot2)    # for ggplot
df$days <- as.numeric(df$Date - df[1,]$Date)
# model based on a list of parameters
f <- function(pars, xx) {pars$a + pars$b*log(pars$tc - xx)} 
# residual function
resids <- function(p, observed, xx) {df$Y - f(p,xx)}
# fit using Levenberg-Marquardt algorithm
nls.out <- nls.lm(par=list(a=1,b=-1,tc=5000), fn = resids, 
                  observed = df$Y, xx = df$days)
# use output of L-M algorithm as starting estimates in nls(...)
par <- nls.out$par
nls.final <- nls(Y~a+b*log(tc-days),data=df, 
                 start=c(a=par$a, b=par$b, tc=par$tc))
summary(nls.final)      # display statistics of the fit 
# append fitted values to df
df$pred <- predict(nls.final)
# plot the results
ggplot(df)+
  geom_line(aes(x=Date,y=Y),color="black")+
  geom_line(aes(x=Date,y=pred),color="blue",linetype=2)+
  labs(title="LPPL Model Applied to DJIA (1982 - 1987)",
       x="", y="DJIA (daily close)")+
  theme(plot.title=element_text(face="bold"))