maxbiostat / truncation_tests

Preliminary code to test truncation algorithms for infinite series and their applications in Statistics.
GNU General Public License v3.0
2 stars 1 forks source link

Doubling algorithm is sensitive to initial guess on number of iterations #5

Open maxbiostat opened 3 years ago

maxbiostat commented 3 years ago

This script shows this.

wellington36 commented 7 months ago

The script mentioned above was deleted in commit https://github.com/maxbiostat/truncation_tests/commit/75032857db5b71debf4bac421b47977789a68ea3 (from Apr 16, 2021).

maxbiostat commented 6 months ago

The script mentioned above was deleted in commit 7503285 (from Apr 16, 2021).

True. Went back and recovered this


negativeBinomial_marginalised <- function(k, theta){
  mu <- theta[1]
  phi <- theta[2]
  eta <- theta[3]
  x <- theta[4]
  ans <- ifelse(k < x,  -Inf, 
                dnbinom(k, mu = mu, size = phi, log = TRUE) + dbinom(x = x, size = k, prob = eta, log = TRUE))
  return(ans)
}

Mu <- 50
Phi <- .5
Eta <- .008
obsX <- 20
Theta <- c(Mu, Phi, Eta, obsX)
Eps <- 1E-16
mit <- 1E5
(TrueValue <- dnbinom(x = obsX, mu = Eta*Mu, size = Phi, log = TRUE) )

d20 <-  approx_doubling(negativeBinomial_marginalised, theta = Theta, N_start = 20,
                        n0 = obsX, epsilon = Eps, max_iter = mit)
d100 <- approx_doubling(negativeBinomial_marginalised, theta = Theta, N_start = 100,
                        n0 = obsX, epsilon = Eps, max_iter = mit)
d500 <- approx_doubling(negativeBinomial_marginalised, theta = Theta, N_start = 500,
                      n0 = obsX, epsilon = Eps, max_iter = mit)

d20$iter
d100$iter
d500$iter

robust_difference(d20$sum, TrueValue)
robust_difference(d100$sum, TrueValue)
robust_difference(d500$sum, TrueValue)

all.equal(d20$sum, TrueValue, tolerance = 0)
all.equal(d100$sum, TrueValue, tolerance = 0)
all.equal(d500$sum, TrueValue, tolerance = 0)

Not sure how relevant this is anymore. But may be interesting to investigate whether guessing too low in the initial thing for doubling has an impact.