Issues implementing recursive function in R - fisher scoring

user_12996569 0
Viewed 58

I am trying to implement Fisher Scoring on simulated i.i.d. Poisson data, but I am getting a stack overflow error. I did some simple prints from the function and found that the guess value is not changing after the first iteration.

fs_pois <- function(data, true, guess) {
if (abs(guess-true) < 0.01) return(guess)
else return(fs_pois(data, true, guess + (mean(data)*exp(-guess))-1))
}

data <- rpois(100,lambda=0.34)
fs_pois(data, 0.34, 0.36)

produces

Error: C stack usage  7970416 is too close to the limit

Does anyone know what is going on? Much appreciated

recursion r
2 Answers

You are updating guess with (mean(data)*exp(-guess))-1, what does not use true what you use in abs(guess-true) for the decision to exit. In your case the exit condition is not reached after 7970416 steps and typical the part (mean(data)*exp(-guess))-1 is getting to 0. Maybe you should use (mean(data)*exp(-guess))-1 for the exit decision.

fs_pois <- function(data, guess) {
  d <- (mean(data)*exp(-guess))-1
  if (abs(d) < 0.01) return(guess)
  fs_pois(data, guess + d)
}

set.seed(7)
data <- rpois(100,lambda=0.34)
fs_pois(data, 0.36)
#[1] -0.9326207
user_104885040

Thank you for your suggestions. After thinking about convergence, I looked back at the mathematical model and found that I had left out a log and in addition confused convergence to true parameter value with convergence to the ML estimate. My revised code is below:

fs_pois <- function(data, MLE, guess) {
  cat(guess,"\n")
  if (abs(guess-MLE) < 0.01) return(guess)
  else return(fs_pois(data, MLE, guess + (mean(data)*exp(-guess))-1))
}
set.seed(111)
data <- rpois(25,lambda=2)
fs_pois(data, log(mean(data)), -1)

As a general note to anyone reading: this is probably not how one would ever implement fisher scoring. This was simply for practicing recursion in R.

user_129965690
Related