fix issues with 'PF_get_score_n_hess' with 'use_O_n_sq = TRUE'

R/PF.R

@@ -1806,7 +1806,6 @@ PF_get_score_n_hess <- function(object, debug = FALSE, use_O_n_sq = FALSE){
     n_rng <- NCOL(Q)
     drng <- 2L * n_rng * n_rng
     org_dim <- dfix + drng
-    stopifnot(length(cpp_res$score) == org_dim)
 
     # output dimension and matrix to multiply objects by
     out_dim <- dfix + n_rng * n_rng + (n_rng * (n_rng + 1L)) / 2L

README.Rmd

@@ -34,7 +34,7 @@ devtools::install_github("boennecd/dynamichazard")
 You can also download the package from CRAN by calling:
 
 ```R
-installed.packages("dynamichazard")
+install.packages("dynamichazard")
 ```
 
 ## Example - ddhazard

examples/firstOrderVAR.Rmd

@@ -543,9 +543,9 @@ This is important for long time series.
 
 ```{r O_N_sq_sgd, cache = 1, dependson = "sgd"}
 # number of iterations
-n_runs <- 600L
+n_runs <- 400L
 # number of particles to use
-ns <- ceiling(exp(seq(log(100), log(500), length.out = n_runs)))
+ns <- ceiling(exp(seq(log(200), log(1000), length.out = n_runs)))
 
 out <- sgd(n_runs = n_runs, ns = ns, use_O_n_sq = TRUE,  object = pf_start, 
            lr = .005)

src/PF/PF_score_n_Hess.cpp

@@ -110,7 +110,7 @@ class score_n_hess : public score_n_hess_base {
   const arma::mat &get_hess_terms() const override {
     return hess_terms;
   }
-  const double get_weight() const override {
+  double get_weight() const override {
     return weight;
   }
 
@@ -417,7 +417,7 @@ namespace {
     const arma::mat &get_hess_terms() const override {
       return hess_terms;
     }
-    const double get_weight() const override {
+    double get_weight() const override {
       return weight;
     }
 
@@ -471,11 +471,6 @@ score_n_hess_O_N_sq::score_n_hess_O_N_sq(
             &score_dim);
       }
     }
-
-    if(!only_score)
-      F77_CALL(dsyr)(
-        &C_U, &dfixd, &D_ONE, score.memptr(), &I_ONE, hess_terms.memptr(),
-        &score_dim);
   }
 
   /* take a copy of the gradient w.r.t. the terms from the observation equation */
@@ -515,7 +510,8 @@ score_n_hess_O_N_sq::score_n_hess_O_N_sq(
       if(!is_first_it){
         score_terms(obs_span)  = old_res->get_score()(obs_span);
         score_terms(sta_span) += old_res->get_score()(sta_span);
-      }
+      } else
+        score_terms(obs_span).zeros();
 
       /* add terms to score */
       score += w_i * score_terms;
@@ -532,8 +528,8 @@ score_n_hess_O_N_sq::score_n_hess_O_N_sq(
           arma::kron(
             dat.K, w_i * ((- innovation_Qinv) * innovation_Qinv.t() + dat.K_1_2));
 
-        /* add outer product of score terms from state equation and old the
-         * state's score vector */
+        /* add outer product of score terms from this pair */
+        score_terms(obs_span) += obs_score_term;
         F77_CALL(dsyr)(
             &C_U, &score_dim, &w_i, score_terms.memptr(), &I_ONE,
             hess_terms.memptr(), &score_dim);
@@ -545,18 +541,6 @@ score_n_hess_O_N_sq::score_n_hess_O_N_sq(
   }
 
   if(!only_score){
-    /* we lag one outer product in the Hessian term between the score terms
-     * from the observation equation and the sum of the weighted score terms
-     * from the previous state and the state equation */
-    {
-      arma::vec tmp = score;
-      tmp(obs_span) -= obs_score_term;
-      R_BLAS_LAPACK::dger(
-        &dfixd, &score_dim, &D_ONE, obs_score_term.memptr(),
-        &I_ONE, tmp.memptr(), &I_ONE, hess_terms.memptr(), &score_dim);
-
-    }
-
    /* subtract outer product of score */
    F77_CALL(dsyr)(
      &C_U, &score_dim, &D_NEG_ONE, score.memptr(), &I_ONE, hess_terms.memptr(),
@@ -618,7 +602,7 @@ std::vector<std::unique_ptr<score_n_hess_base> > PF_get_score_n_hess_O_N_sq_comp
 
   /* run particle filter and compute smoothed functionals recursively */
   for(unsigned t = 1;  t <= (unsigned)data.d; ++t, ++r_obj){
-    if((t + 1L) % 3L == 0L)
+    if(t % 5L == 0L)
       Rcpp::checkUserInterrupt();
 
     std::shared_ptr<PF_cdist> y_dist = dens_calc.get_y_dist(t);

src/PF/PF_utils.h

@@ -246,7 +246,7 @@ class score_n_hess_base {
   virtual const arma::vec &get_score() const = 0;
   /* the betas in https://doi.org/10.1093/biomet/asq062 */
   virtual const arma::mat &get_hess_terms() const = 0;
-  virtual const double get_weight() const = 0;
+  virtual double get_weight() const = 0;
 
   virtual ~score_n_hess_base() = default;
 };

src/PF/est_params_dens.cpp

@@ -7,7 +7,7 @@ using pair_iterator = std::vector<const particle_pairs*>::const_iterator;
 class update_parameters_data {
   /* function to create a list of pointers to be used in loops later */
   std::vector<const smoother_output::particle_pairs*>
-  set_pairs(){ /* TODO: Just store two start and end iterators */
+  set_pairs(){ /* TODO: Just store a start and end iterators */
     std::vector<const particle_pairs*>::size_type total_n_pairs = 0L;
     for(auto i = tr->begin(); i != tr->end(); ++i)
       total_n_pairs += i->size();
@@ -16,7 +16,7 @@ class update_parameters_data {
     auto ptr = out.begin();
     for(auto i = tr->begin(); i != tr->end(); ++i)
       for(auto j = i->begin(); j != i->end(); ++j, ++ptr)
-        *ptr = &(*j);
+        *ptr = &*j;
 
     return out;
   }
@@ -46,12 +46,12 @@ class generator_dens final : public qr_data_generator {
   const update_parameters_data &dat;
   pair_iterator i_start;
   pair_iterator i_end;
-  const unsigned long int n_pairs;
+  const std::size_t n_pairs;
 
 public:
   generator_dens
     (const update_parameters_data &dat, pair_iterator i_start,
-     pair_iterator i_end, const unsigned long int n_pairs):
+     pair_iterator i_end, const std::size_t n_pairs):
     dat(dat), i_start(i_start), i_end(i_end), n_pairs(n_pairs) { }
 
   qr_work_chunk get_chunk() const override {
@@ -104,7 +104,7 @@ PF_parameters
 
     /* setup generators */
     const unsigned long max_per_block =
-        max_bytes / ((dat.n_elem_X + dat.n_elem_Y) * 8L) + 1L;
+        max_bytes / (dat.n_elem_X + dat.n_elem_Y) / 8L + 1L;
 
     if(debug)
       Rcpp::Rcout << "Running `est_params_dens` with `max_per_block` "
@@ -159,7 +159,7 @@ PF_parameters
     arma::inplace_trans(out.R_top_F);
 
     /* use that
-     *    Q = Y^\top Y - F^\top F
+     *    Q = (Y^\top Y - F^\top F) / d
      *
      * where F is the output from qr_parallel.compute().
      */