improve numerical robustness in all PF algos

In auxiliary filter in particular, but presumably others, we exponentiate the weights and then sum and divide by the sum of weights. We also exponeniate and then take the log to get the log-likelihood. This is not robust to numerical underflow. We can presumably address this using the log-sum-exp trick to avoid exponentiating overly large magnitude values.

Here's a reprex from a user:

library(nimble)
library(nimbleSMC)

#### Model
SI <- nimbleCode({
  S[1] <- iniState[1]
  I[1] <- iniState[2]
  y[1] ~ dpois(pdet * I[1] + 1e-5)
  for (t in 2:nT) {
    lambda[t] <- beta * I[t-1] / N * S[t-1] + 1e-5
    negdS[t] ~ dpois(lambda[t])

    S[t] <- S[t-1] - negdS[t]
    I[t] <- I[t-1] + negdS[t]

    y[t] ~ dpois(pdet * I[t] + 1e-5)
    X[1:2,t] <- c(S[t], I[t])
  }

  beta ~ dunif(0, 0.3) #Transimissibility
  pdet ~ dunif(0, 1) #Probability of detection
})

states <- c(10000, 5)
constants <- list(
  N = sum(states),
  nT = 50,
  iniState = states
)

y <- c(0, 3, 2, 2, 8, 9, 3, 6, 7, 5, 6, 10, 14, 9, 8, 16, 
       15, 20, 27, 29, 37, 33, 32, 52, 55, 67, 56, 68, 85, 85, 111, 
       124, 145, 169, 206, 203, 260, 277, 320, 335, 414, 438, 539, 621, 
       681, 783, 932, 1015, 1099, 1222)

SIMod <- nimbleModel(code = SI,
                     constants = constants,
                     data = list(y = y),
                     inits = list(beta = 0.2, pdet = 0.5))

my_AuxF <- buildAuxiliaryFilter(SIMod, 'I',
                                control = list(saveAll = TRUE))
SIModC <- compileNimble(SIMod)
my_AuxFC <- compileNimble(my_AuxF, project = SIMod)
logLik <- my_AuxFC$run(m = 10000)
logLik # Always -Inflibrary(nimble)
library(nimbleSMC)

#### Model
SI <- nimbleCode({
  S[1] <- iniState[1]
  I[1] <- iniState[2]
  y[1] ~ dpois(pdet * I[1] + 1e-5)
  for (t in 2:nT) {
    lambda[t] <- beta * I[t-1] / N * S[t-1] + 1e-5
    negdS[t] ~ dpois(lambda[t])

    S[t] <- S[t-1] - negdS[t]
    I[t] <- I[t-1] + negdS[t]

    y[t] ~ dpois(pdet * I[t] + 1e-5)
    X[1:2,t] <- c(S[t], I[t])
  }

  beta ~ dunif(0, 0.3) #Transimissibility
  pdet ~ dunif(0, 1) #Probability of detection
})

states <- c(10000, 5)
constants <- list(
  N = sum(states),
  nT = 50,
  iniState = states
)

y <- c(0, 3, 2, 2, 8, 9, 3, 6, 7, 5, 6, 10, 14, 9, 8, 16, 
       15, 20, 27, 29, 37, 33, 32, 52, 55, 67, 56, 68, 85, 85, 111, 
       124, 145, 169, 206, 203, 260, 277, 320, 335, 414, 438, 539, 621, 
       681, 783, 932, 1015, 1099, 1222)

SIMod <- nimbleModel(code = SI,
                     constants = constants,
                     data = list(y = y),
                     inits = list(beta = 0.2, pdet = 0.5))

my_AuxF <- buildAuxiliaryFilter(SIMod, 'I',
                                control = list(saveAll = TRUE))
SIModC <- compileNimble(SIMod)
my_AuxFC <- compileNimble(my_AuxF, project = SIMod)
logLik <- my_AuxFC$run(m = 10000)
logLik # Always -Inf

nimble-dev / nimbleSMC

improve numerical robustness in all PF algos #3