0

晚上好
，我目前正在使用 R 中的 nlminb() 函数解决最大化问题。这是代码的一部分。这些是初始输入参数

w <- rnorm(12,0.5,1)
y <- 1:12
x <- rnorm(12,0,1)
h <- 0.25
n <- length(x)

g <- sd(x)

将评估参数的两个函数中的第一个

 Resid <- function(par, data) 
    {
      alpha.0 <- par[1] 
      alpha.1 <- par[2]
      beta.1  <- par[3]
      mu1     <- par[4]
      mu2     <- par[5]

  n <- length(x)
  sigma.sqs <- numeric(n) 
  epsilon <- numeric(n) 
  sigma.sqs[1] <- g 
  epsilon[1] = x[1] - mean(x)
  for(ii in c(1:(n-1))) {
    sigma.sqs[ii + 1] <- (
      alpha.0 +  
        alpha.1 * (epsilon[ii])^2 +
        beta.1 * sigma.sqs[ii])
  epsilon[ii+1] <- (x[ii+1]-mu1-mu2*x[ii])/sigma.sqs[ii]

    Ksum <- 0
    for(j in (1:(n-1))){
      Ksum <- Ksum + (((epsilon[ii]/(sigma.sqs[ii]^0.5))-w[j])/h)
    }
  }

  return(list(et = epsilon, ht = sigma.sqs, xt=Ksum)) 
}

第二部分，从函数 Resid 中​​获取 sigma 和 epsilons

 LogL <- function(par, data) {

      res <- Resid(par, data) 
      sigma.sqs <- res$ht
  epsilon <- res$et  
  f <- res$xt

 return( 1/n * sum( log(1/(n*h)*(1/((2*pi)^0.5))*exp(-0.5*(f)^2)) +    log(1/(sigma.sqs^0.5))))

}

最后最大化

 o <- nlminb(start=c(0.001,0.001,0.001,0.001,0.001), objective= LogL, lower=  0.0000001  ) 
    print(o)

代码运行，但它出现了 NaN。问题似乎出现在for循环中

        Ksum <- 0
        for(j in (1:(n-1))){
          Ksum <- Ksum + (((epsilon[ii]/(sigma.sqs[ii]^0.5))-w[j])/h)
        }

该循环应该计算 Ksum 的向量，每个 x 一个向量。我一直在努力找出问题所在，但我已经对解决方案视而不见。

有任何想法吗？

干杯

4

2 回答 2

1

实际上，您的问题似乎出在 LogL 函数中。那是第一个 NaN 值似乎出现的地方，它们似乎来自exp(-0.5*(f)^2)该值生成的术语Ksum。问题（至少在我运行它时）是 f 越来越小到 R 返回的程度exp(-0.5*(-313.5329)^2)=0，然后你记录了你得到 NaN 值的日志。

所以为了让它在数值上更稳定，我曾经log(a*b)=log(a)+log(b)重写过这个函数。你会想验证我写的东西在数学上是相同的，但它似乎不太可能产生越界问题。

return( 1/n * sum( 
    log( 1/(n*h) ) +  log( 1/(2*pi)^0.5 ) + -0.5*(f)^2  +
    log(1/(sigma.sqs^0.5))
))
于 2014-05-07T01:57:36.730 回答
0

您说您希望循环计算 Ksum 的向量，每个 x 一个向量。

您的循环实际上所做的是为每个 x 计算一次标量 Ksum，每次都覆盖它，然后将最终 x 生成的单个 Ksum 值传递给 LogL 函数。

看起来你需要改变

    Ksum <- 0

至

    Ksum <- numeric(n)

并将那行代码移到更大的 for 循环之外，for(ii in c(1:(n-1)))这样您就不会在 x 的每个新值中覆盖计算值。

您还需要更改此行

    Ksum <- Ksum + (((epsilon[ii]/(sigma.sqs[ii]^0.5))-w[j])/h)

引用您想要的任何索引，例如

    Ksum[ii] <- Ksum[ii] + (((epsilon[ii]/(sigma.sqs[ii]^0.5))-w[j])/h)

希望这会有所帮助。

于 2014-05-07T01:00:41.130 回答