dynamichazard (version 0.6.4)

PF_get_score_n_hess: Approximate Negative Observation Matrix and Score Vector

Description

Returns a list of functions to approximate the negative observation matrix and score vector.

Usage

PF_get_score_n_hess(object, debug = FALSE, use_O_n_sq = FALSE)

Arguments

object

object of class PF_EM.

debug

TRUE if debug information should be printed to the console.

use_O_n_sq

TRUE if the method from Poyiadjis et al. (2011) should be used.

Value

A list with the following functions as elements

run_particle_filter

function to run particle filter as with PF_forward_filter.

set_parameters

function to set the parameters in the model. The first argument is a vectorized version of \(F\) matrix and \(Q\) matrix. The second argument is the fixed effect coefficients.

set_n_particles

sets the number of particles to use in run_particle_filter and get_get_score_n_hess when use_O_n_sq is TRUE.

get_get_score_n_hess

computes the approximate negative observation matrix and score vector. The argument toggles whether the approximate negative observation matrix should be computed. The last particle cloud from run_particle_filter is used when use_O_n_sq is FALSE.

Warning

The function is still under development so the output and API may change.

Details

The score vector and negative observed information matrix are computed with the (forward) particle filter. This comes at an \(O(d^2)\) variance where \(d\) is the number of periods. Thus, the approximation may be poor for long series. The score vector can be used to perform stochastic gradient descent.

If use_O_n_sq is TRUE then the method in Poyiadjis et al. (2011) is used. This may only have a variance which is linear in the number of time periods. However, the present implementation is \(O(N^2)\) where \(N\) is the number of particles. The method uses a particle filter as in Section 3.1 of Lin et al. (2005). There is no need to call run_particle_filter unless one wants a new approximation of the log-likelihood as a separate filter is run with get_get_score_n_hess when use_O_n_sq is TRUE.

References

Cappe, O. and Moulines, E. (2005) Recursive Computation of the Score and Observed Information Matrix in Hidden Markov Models. IEEE/SP 13th Workshop on Statistical Signal Processing.

Cappe, O., Moulines, E. and Ryden, T. (2005) Inference in Hidden Markov Models (Springer Series in Statistics). Springer-Verlag.

Doucet, A., and Tadi<U+0107>, V. B. (2003) Parameter Estimation in General State-Space Models Using Particle Methods. Annals of the Institute of Statistical Mathematics, 55(2), 409<U+2013>422.

Lin, M. T., Zhang, J. L., Cheng, Q. and Chen, R. (2005) Independent Particle Filters. Journal of the American Statistical Association, 100(472), 1412-1421.

Poyiadjis, G., Doucet, A. and Singh, S. S. (2011) Particle Approximations of the Score and Observed Information Matrix in State Space Models with Application to Parameter Estimation. Biometrika, 98(1), 65--80.

See Also

See the examples at https://github.com/boennecd/dynamichazard/tree/master/examples.

Examples

# NOT RUN {
library(dynamichazard)
.lung <- lung[!is.na(lung$ph.ecog), ]
# standardize
.lung$age <- scale(.lung$age)

set.seed(43588155)
pf_fit <- PF_EM(
  fixed = Surv(time, status == 2) ~ ph.ecog + age,
  random = ~ age, model = "exponential",
  data = .lung, by = 50, id = 1:nrow(.lung),
  Q_0 = diag(1, 2), Q = diag(.5^2, 2), type = "VAR",
  max_T = 800,
  control = PF_control(
    N_fw_n_bw = 250, N_first = 2000, N_smooth = 500, covar_fac = 1.1,
    nu = 6, n_max = 1000L, eps = 1e-5, est_a_0 = FALSE, averaging_start = 100L,
    n_threads = max(parallel::detectCores(logical = FALSE), 1)))

comp_obj <- PF_get_score_n_hess(pf_fit)
comp_obj$set_n_particles(N_fw = 10000L, N_first = 10000L)
comp_obj$run_particle_filter()
(o1 <- comp_obj$get_get_score_n_hess())

# O(N^2) method with lower variance
comp_obj <- PF_get_score_n_hess(pf_fit, use_O_n_sq = TRUE)
comp_obj$set_n_particles(N_fw = 2500L, N_first = 2500L)
(o2 <- comp_obj$get_get_score_n_hess())

# approximations may have large variance
o3 <- replicate(10L, {
  runif(1)
  pf_fit$seed <- .Random.seed
  comp_obj <- PF_get_score_n_hess(pf_fit)
  comp_obj$set_n_particles(N_fw = 10000L, N_first = 10000L)
  comp_obj$run_particle_filter()
  comp_obj$get_get_score_n_hess()
}, simplify = FALSE)
sapply(o3, function(x) x$observation$score)
sapply(o3, function(x) sqrt(diag(solve(x$observation$neg_obs_info))))
# }
# NOT RUN {
# }