r - NLS 逻辑回归的 optimx 问题
问题描述
我正在尝试根据历史数据(y)估计函数的一些参数。这是具有异步更新 (delta) 的逻辑回归。
这是有问题的系统:
其中(出于我们的目的)X[t] = y[t],n2 表示为 (1-n1),phi1 和 phi2 表示为 phi1+deltaphi(如 phi2>phi1),u[t] 是一个隐式错误由 y_hat 捕获的术语。
数据(日志值):
y <- c(-0.083522212, -0.080744273, -0.098003453, -0.090994700, -0.105010991, -0.112070623, -0.115681762, -0.143194134, -0.146458642, -0.139691305, -0.128929970, -0.118047088, -0.095065509, -0.082399946, -0.100997872, -0.092699501, -0.082550517, -0.050470001, 0.030390893, 0.131429122, 0.180958369, 0.212374498, 0.223668346, 0.226461144, 0.209141361, 0.195377626, 0.178487458,
0.201981948, 0.233653604, 0.245221474, 0.227886405, 0.141274238, 0.046683795, -0.047819717, -0.112630561, -0.203788442, -0.238529171, -0.211924261, -0.233738086, -0.241872522, -0.238041656, -0.230753558, -0.242931741, -0.231894162, -0.243119603, -0.233052377, -0.230820606, -0.225594126, -0.232095554, -0.244800121, -0.252265025, -0.241778694, -0.227898251, -0.242656156, -0.229516117, -0.216082812, -0.220941314, -0.211800617, -0.183642284, -0.165779424, -0.159285263, -0.147407410, -0.138607996, -0.130455753, -0.094857132, -0.039392141, -0.003361144, 0.076381508, 0.101627405, 0.103042608, 0.096997308, 0.100308333, 0.098658702, 0.083952591, 0.077534743, 0.064491677, 0.056002466, 0.082643906, 0.080460147, 0.090688462)
估计参数:
phi1, deltaphi, beta, delta
其他的东西:
T <- 80
y_hat <- mat.or.vec(nr = T, nc = 1)
R = 1.019656
alpha_bar = 0
n1 = mat.or.vec(nr = T, nc = 1)
功能:
f <- function(t, phi1, deltaphi, beta, delta, R, alpha_bar, y) {
n1[t] <<- (delta*n1[t-1])+(1-delta)*(1/(1+exp(-beta*((y[t-1]+alpha_bar-R*y[t-2])*((-deltaphi))*y[t-3]))))
y_hat = (((n1[t]*phi1+(1-n1[t])*(phi1+deltaphi))/((R+alpha_bar)))*y[t-1])
return((y[t]-y_hat)^2)
}
func <- function(par) sum(sapply(4:T,f, par[1],par[2], par[3],par[4], R, alpha_bar, y))
fit <- optimx(c(0.9, 1.05, 0.05, 0.6),
method = "nlm",
func,
hessian = TRUE)
我认为 n1 是这里的问题,我尝试了很多技巧来使其符合要求,因为它是一个没有起点的递归变量,因此是超级运算符。我有一些运气,但该功能总是在我身上崩溃,所以我认为该功能可能以某种方式错误指定。对 beta 和 delta 有一些限制 - beta 应该是非零的正数,而 delta 应该在 0 和 1 之间。
解决方案
我可以用不同的优化器得到一些答案(但不能用nlm
);它确实看起来您的目标函数可能被过度指定,不同的优化器会以不同的方式处理......我试图稍微简化目标函数。
设置常量
maxT <- 80
R <- 1.019656
alpha_bar <- 0
预测功能
我过去常常with(as.list(par) ...)
避免解包参数向量。这确实使调试变得更加困难,并且确实迫使我们<<-
在表达式中使用...
predfun <- function(par) {
y_hat <- n1 <- numeric(maxT) ## pre-allocate vectors
with(as.list(par), {
for (t in (4:maxT)) {
n1[t] <<- (delta*n1[t-1])+(1-delta)*
(1/(1+exp(-beta*((y[t-1]+alpha_bar-R*y[t-2])*
((-deltaphi))*y[t-3]))))
y_hat[t] <<- (((n1[t]*phi1+(1-n1[t])*(phi1+deltaphi))/
((R+alpha_bar)))*y[t-1])
}
})
return(y_hat)
}
目标函数
func <- function(par) {
return(sum((y[4:maxT]-predfun(par)[4:maxT])^2))
}
优化
p0 <- c(phi1=0.9, deltaphi=1.05, beta=0.05, delta=0.6)
func(p0) ## check to make sure the function works for starting values
library(optimx)
fit <- optimx(par=p0,
method = c("Nelder-Mead","BFGS","nlm"),
func,
hessian = TRUE)
结果
phi1 deltaphi beta delta value fevals
Nelder-Mead -0.4967063 3.026097 6.227032 -0.16820934 7.546323e-02 265
BFGS -0.4462906 2.903375 2.554713 -0.01802873 7.922255e-02 110
nlm NA NA NA NA 8.988466e+307 NA
gevals niter convcode kkt1 kkt2 xtime
Nelder-Mead NA NA 0 TRUE TRUE 0.143
BFGS 100 NA 1 FALSE FALSE 0.494
nlm NA NA 9999 NA NA 0.001
nlm 给出所有 NA 值;BFGS 不符合两个 KKT 标准;Nelder-Mead 看起来不错。Nelder-Mead 和 BFGS 给出了适度不同的参数估计,尽管两者都给出了明显小的 SSQ。
cc <- coef(fit)
plot(y)
lines(predfun(cc["Nelder-Mead",]), col=2)
lines(predfun(cc["BFGS",]), col=4)
拟合线(几乎)无法区分。
转换参数
如果您想在变换尺度上拟合参数,您可以在变换尺度上为优化函数提供初始参数,例如
p0 <- c(phi1=0.9, deltaphi=1.05, log_beta=log(0.05), logit_delta=qlogis(0.6))
然后在您的目标函数中,在使用它们之前对它们进行反向转换:
predfun <- function(par) {
y_hat <- n1 <- numeric(maxT) ## pre-allocate vectors
with(as.list(par), {
delta <- plogis(logit_delta)
beta <- exp(log_beta)
### ... then the rest of your objective function ...
推荐阅读
- wordpress - 按类别设置最小订单数量(Woocommerce)
- javascript - 使用 Ajax 调用刷新 JSP 页面上的 div 而没有响应
- hive - hive 独立元存储读取 avro 数据且架构不起作用
- amazon-web-services - 从 python 服务器订阅 aws-iot 主题
- javascript - 当我打电话给他时,我的 php 脚本返回 403 错误
- python - 下面的程序要改成多线程应该做些什么?
- javascript - Mongoose/Node:将元素推入数组类型的文档字段?
- javascript - Ajax 调用未在模式弹出窗口中触发:控制台中未显示错误
- python - pyshark - 如何在实时捕获期间打印目标 IP?
- python - Databricks 火花 UDF 不适用于过滤的数据帧