拟合局部级别泊松（状态空间模型）

Question

我正在研究“使用指数平滑进行预测”。我被困在练习 16.4 的部分：

该数据集partx包含汽车零件的月销售历史记录。应用局部泊松模型。应该通过最大化似然或最小化平方误差总和来估计参数。

局部泊松模型定义为：

在哪里 和

我有以下代码，但似乎卡住了。优化总是返回接近起始值的东西。

我是否正确拟合了局部泊松模型？

library(expsmooth)
data("partx")
S <- function(x) {
  a <- x[1]
  if(a < 0 | a > 1)
    return(Inf)
  n <- length(partx)
  lambda <- numeric(n+1)
  error <- numeric(n)
  lambda[1] <- x[2]
  for(i in 1:n) {
    error[i] <- partx[i]-rpois(1,lambda[i])
    lambda[i+1] <-   (1-a)*lambda[i] + a*partx[i]
  }
  return(sum(error^2))
}

# returns a = 0.5153971 and lambda = 5.9282414
op1 <- optim(c(0.5,5L),S, control = list(trace = 1))
# returns a = 0.5999655 and lambda = 2.1000131
op2 <- optim(c(0.5,2L),S, control = list(trace = 1))

Answer 1

我知道这本书说你可以使用误差平方和或 MLE，但第一个选项似乎也是有线的，因为你必须对毒物分布进行采样，所以如果你修复参数，你每次都会得到不同的误差平方和。由于您没有说您已经尝试过 MLE 方法，所以我对其进行了编程。数学计算如下：

实现它的代码是

MLELocalPoisson = function(par,y){
  alpha = par[1]
  lambda_ini = par[2]
  n = length(y)
  vec_lambda = rep(NA, n)
  for(i in 1:n){
    if(i==1){
      vec_lambda[i] = (1-alpha)*lambda_ini+alpha*y[i]
    }else{
      vec_lambda[i] = (1-alpha)*vec_lambda[i-1]+alpha*y[i]
    }
  }
  vec_lambda = c(lambda_ini,vec_lambda[-n])
  sum_factorial = sum(sapply(y,function(x)log(factorial(x))))
  sum_lambda = sum(vec_lambda)
  sum_prod = sum(log(vec_lambda)*y)
  loglike = -sum_prod+sum_lambda+sum_factorial
  return(loglike)
}
optim(par = c(0.1,1),fn = MLELocalPoisson,y = partx, method = "L-BFGS-B",
      lower=c(1e-10,1e-10),upper = c(1,Inf),control = list(maxit = 10000))

较低的值集 a1e-10已完成，因此优化不会尝试c(0,0)，从而生成 的对数似然NaN。

编辑

查看泊松回归文献，通常定义 $\lambda = exp(x*\beta)$ 并将残差计算为 $y-exp(x*\beta)$ （看看）。因此，可以使用作者给出的 $\lambda$ 公式来解决这个问题，如下所示：

LocalPoisson = function(par,y,optim){
  alpha = par[1]
  lambda_ini = par[2]
  n = length(y)
  vec_lambda = rep(NA, n)
  y_hat = rep(NA, n)
  for(i in 1:n){
    if(i==1){
      vec_lambda[i] = (1-alpha)*lambda_ini+alpha*y[i]
    }else{
      vec_lambda[i] = (1-alpha)*vec_lambda[i-1]+alpha*y[i]
    }
  }
  if(optim){
    y_hat = c(lambda_ini,vec_lambda[-n])
    return(sum((y_hat-y)^2))
  } else {
    return(data.frame(y_hat = y_hat,y=y, lambda = vec_lambda))
  }

}
optim(par = c(0.1,1),fn = LocalPoisson,y = partx, optim =T,method = "L-BFGS-B",
                lower=c(1e-10,1e-10),upper = c(1,Inf),control = list(maxit = 10000))

它不会产生与 MLE 相同的结果（我对这个选项感觉更舒服，但它可能是估计参数的一种可能方法）。