r - 拟合局部级别泊松(状态空间模型)

Question

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

The data set partx contains a history of monthly sales of an automobile part. Apply a local Poisson model. Parameters should be estimated by either maximizing the likelihood or minimizing the sum of squared errors.

局部泊松模型定义为:



哪里 和

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

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

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))

较低的值设置了 1e-10做完所以优化不要尝试 c(0,0)从而产生 NaN 的对数似然.

编辑

查看泊松回归文献，通常定义 $\lambda = exp(x*\beta)$ 并将残差计算为 $y-exp(x*\beta)$ ( have a look at )。因此，可以使用作者为 $\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 相同的结果(我对那个选项感觉更舒服，但它可能是一种估计参数的可能方法)。