fix non-const phi

NEWS

@@ -8,6 +8,8 @@ bssm 1.1.6 (Release date: )
      ar1_lg model when predicting past or when using simulation smoother.
    * Fixed a bug which caused an error when predicting past values in 
      multivariate time series case.
+   * Fixed log-likelihood computation for gamma model with non-constant shape 
+     parameter when using (intermediate) Gaussian approximation.
    * Fixed sampling of negative binomial distribution in predict method, which 
      used std::negative_binomial which converts non-integer phi to integer. 
      Sampling now uses Gamma-Poisson mixture for simulation.

src/distr_consts.cpp

@@ -26,7 +26,7 @@ double negbin_log_const(double y, double u, double phi) {
 }
 
 double gamma_log_const(double y, double u, double phi) {
-  return phi * std::log(phi) - std::lgamma(phi) + (phi - 1) * std::log(y) - std::log(u);
+  return phi * std::log(phi) - std::lgamma(phi) + (phi - 1) * std::log(y) - phi * std::log(u);
 }
 
 
@@ -61,7 +61,7 @@ double negbin_log_const(const arma::vec&  y, const arma::vec& u, double phi) {
 double gamma_log_const(const arma::vec&  y, const arma::vec& u, double phi) {
   double res = 0;
   for(unsigned int i = 0; i < y.n_elem; i++) {
-    res += phi * std::log(phi) - std::lgamma(phi) + (phi - 1) * std::log(y(i)) - std::log(u(i));
+    res += phi * std::log(phi) - std::lgamma(phi) + (phi - 1) * std::log(y(i)) - phi * std::log(u(i));
   }
   return res;
 }

src/model_bsm_ng.cpp

@@ -91,7 +91,7 @@ double bsm_ng::log_prior_pdf(const arma::vec& x, const Rcpp::Function prior_fn)
 
   double log_prior = 0.0;
   arma::vec pars = x;
-  if (arma::accu(fixed) < 3 || noise) {
+  if (arma::accu(fixed) < 3 || noise || phi_est) {
     pars.subvec(0, pars.n_elem - xreg.n_cols - 1) = 
       arma::exp(pars.subvec(0, pars.n_elem - xreg.n_cols - 1));
     // add jacobian

src/model_ssm_mng.cpp

@@ -378,10 +378,10 @@ arma::vec ssm_mng::log_weights(const unsigned int t,  const arma::cube& alpha) c
             std::log(phi(j) + u(j,t) * std::exp(simsignal(j)));
           break;
         case 4 :
-          weights(i) += -phi(j) * simsignal(j) - (y(j,t) * phi(j) * exp(-simsignal(j)) / u(j,t));
+          weights(i) -= phi(j) * (simsignal(j) + (y(j,t) * exp(-simsignal(j)) / u(j,t)));
           break;
         case 5 :
-          weights(i) += -0.5 * std::pow((y(j,t) - simsignal(j)) / phi(j), 2.0);
+          weights(i) -= 0.5 * std::pow((y(j,t) - simsignal(j)) / phi(j), 2.0);
           break;
         }
         weights(i) += 
@@ -429,7 +429,7 @@ arma::vec ssm_mng::log_obs_density(const unsigned int t,
             std::log(phi(j) + u(j,t) * std::exp(simsignal(j)));
           break;
         case 4 :
-          weights(i) += -phi(j) * simsignal(j) - (y(j,t) * phi(j) * exp(-simsignal(j)) / u(j,t));
+          weights(i) += -phi(j) * (simsignal(j) + (y(j,t) * exp(-simsignal(j)) / u(j, t)));
           break;
         case 5 :
           weights(i) += -0.5 * std::pow((y(j,t) - simsignal(j)) / phi(j), 2.0);

src/model_ssm_ung.cpp

@@ -283,7 +283,7 @@ void ssm_ung::update_scales() {
   case 4  :
     for(unsigned int t = 0; t < n; t++) {
       if (arma::is_finite(y(t))) {
-        scales(t) = -phi * mode_estimate(t) - (y(t) * phi * exp(-mode_estimate(t)) / u(t)) +
+        scales(t) = -phi * (mode_estimate(t) + (y(t) * exp(-mode_estimate(t)) / u(t))) +
           0.5 * std::pow((approx_model.y(t) - mode_estimate(t)) / approx_model.H(t), 2.0);
       } 
     }
@@ -426,7 +426,7 @@ arma::vec ssm_ung::log_weights(
       for (unsigned int i = 0; i < alpha.n_slices; i++) {
         double simsignal = arma::as_scalar(D(t * Dtv) + Z.col(t * Ztv).t() *
           alpha.slice(i).col(t) + xbeta(t));
-        weights(i) = -phi * simsignal - (y(t) * phi * exp(-simsignal) / u(t)) +
+        weights(i) = -phi * (simsignal + (y(t) * exp(-simsignal) / u(t))) +
           0.5 * std::pow((approx_model.y(t) - simsignal) / approx_model.H(t), 2.0);
       }
       break;
@@ -483,7 +483,7 @@ arma::vec ssm_ung::log_obs_density(const unsigned int t,
       for (unsigned int i = 0; i < alpha.n_slices; i++) {
         double simsignal = arma::as_scalar(D(t * Dtv) + Z.col(t * Ztv).t() *
           alpha.slice(i).col(t) + xbeta(t));
-        weights(i) = -phi * simsignal - (y(t) * phi * exp(-simsignal) / u(t));
+        weights(i) = -phi * (simsignal + (y(t) * exp(-simsignal) / u(t)));
 
       }
       break;

tests/testthat/test_mcmc.R

@@ -189,6 +189,32 @@ test_that("MCMC results with psi-APF for Poisson model are correct", {
 })
 
 
+
+test_that("MCMC using SPDK for Gamma model works", {
+
+  set.seed(123)
+  n <- 20
+  u <- rgamma(n, 3, 1)
+  phi <- 5
+  x <- cumsum(rnorm(n, 0, 0.5))
+  y <- rgamma(n, shape = phi, scale = u * exp(x) / phi)
+  model_bssm <- bsm_ng(y, 
+    sd_level = gamma(0.1, 2, 10), u = u, phi = gamma(2, 2, 0.1),
+    distribution = "gamma", P1 = 2)
+
+  expect_error(mcmc_gamma <- run_mcmc(model_bssm, sampling_method = "spdk",
+    iter = 1000, particles = 5, seed = 42), NA)
+
+  expect_gt(mcmc_gamma$acceptance_rate, 0)
+  expect_gte(min(mcmc_gamma$theta), 0)
+  expect_lt(max(mcmc_gamma$theta), Inf)
+  expect_true(is.finite(sum(mcmc_gamma$alpha)))
+
+  expect_lt(sum(abs(summary(mcmc_gamma)[,"Mean"] - 
+      c(0.520146284042284, 2.17575390744017))), 0.3)
+
+})
+
 test_that("MCMC results for SV model using IS-correction are correct", {
   set.seed(123)
   expect_error(model_bssm <- svm(rnorm(10), rho = uniform(0.95, -0.999, 0.999),