Implementing Fisher Scoring Recursively

[pocketlemon5]

I am trying to implement Fisher Scoring on simulated i.i.d. Poisson data recursively, but I am getting a stack overflow error. I did some simple tests 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  7969872 is too close to the limit

Does anyone know what is going on?
Much appreciated

Edit: Changing the stack limit does not resolve the issue because the issue is with convergence and not memory. (I.e. the guess value is not changing).

[technocrat]

Hi, and welcome!

Please see the FAQ: What's a reproducible example (`reprex`) and how do I create one? Using a reprex, complete with representative data will attract quicker and more answers.

Your question doesn't really need one (and due to the error, couldn't be produced, anyway) since there are no libraries to be loaded and no data. Still, it helps to enclose in triple backticks to make it easier to cut and past

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)

The code produces

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

This is an operating system constraint that this S/O post addresses for the *nices. Don't think I can help you on windows.

[pocketlemon5]

Thank you for the tips! I have applied them to the issue!

[nirgrahamuk]

As a general rule, I would advise, when performing calculations which rely on convergence, you should test for failure to converge

fs_pois <- function(data, true, guess) {
  cat(guess, " ",(mean(data)*exp(-guess)),"\n")
  if (abs(guess-true) < 0.01) return(guess)
  if (round((mean(data)*exp(-guess)),10) == 1) # i just simplifed the equation modifying guess
    cat("Failed to converge, ", guess,"\n")
  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)

example output:

> fs_pois(data, 0.34, 0.36)
0.36   0.2232564 
-0.4167436   0.4854443 
-0.9312993   0.8120974 
-1.119202   0.9799709 
-1.139231   0.9997967 
-1.139434   1 
-1.139434   1 
Failed to converge,  -1.139434

[pocketlemon5]

Thank you for the help! 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.