For loop not running as intended in R

Question

Good evening
I am currently working on a maximization problem using the nlminb() function in R. Heres part of the code. These are the initial input parameters

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

g <- sd(x)

The first of two functions from which the parameters will be evaluated

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

The second part, taking the sigma and epsilons from the function Resid

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

}

and finally maximizing

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

The code runs, but its coming up with NaNs. The problem seems to be arising in the for-loop regaring

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

The loop is supposed to calculate a vector of Ksum's, one for every x. I have been banging my head trying to find out what is wrong, but i have become blind for the solution.

Any ideas?

Cheers

Answer 1

Actually, it looks like your problem is in the LogL function. That's where the first NaN values seem to appear and they seem to come from the exp(-0.5*(f)^2) term which is generated by the Ksum value. The problem (at least when I ran it) was that f was getting small to the point that R returned exp(-0.5*(-313.5329)^2)=0 and then you took the log of that you got a NaN value.

So to perhaps make it more numerically stable, I used log(a*b)=log(a)+log(b) to re-write the function. You'll want to verify what I wrote is mathematically identical, but it seems less likely to produce the out-of-bounds problem.

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

Answer 2

You said you want the loop to calculate a vector of Ksum's, one for every x.

What your loop actually does is calculate a scalar Ksum once for each x, over-writing it each time, and then passing the single Ksum value generated by the final x to the LogL function.

It looks like you need to change

    Ksum <- 0

to

    Ksum <- numeric(n)

and move that line of code outside the larger for loop for(ii in c(1:(n-1))) so that you're not overwriting your calculated values in each new value of x.

You'll also need to change this line

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

to reference whatever index you want e.g.

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

Hopefully that helps.