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

Have you tried putting the debugger on the Resid code and stepping through it to see what the problem is? Or at least narrow it down? You can call `debugonce(Resid)` and run the code. Then when it hits the Resid function, you can step though it and see how the variables are changing. See `?debug` for more info. — MrFlick, May 07 '14 at 00:23

score 1 · Answer 1 · answered May 07 '14 at 01:57

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

Sean Murphy · Answer 2 · 2014-05-07T01:30:04.733

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.

For loop not running as intended in R

2 Answers2