Quick fix for stan-dev#1495 and the worst part of stan-dev#1496

martinmodrak · Jan 16, 2020 · 8978880 · 8978880
1 parent cfad0f7
commit 8978880
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 28 deletions.
diff --git a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp
@@ -9,6 +9,7 @@
 #include <stan/math/prim/fun/lgamma.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/prob/poisson_log_lpmf.hpp>
 #include <cmath>
 
 namespace stan {
@@ -84,19 +85,25 @@ return_type_t<T_log_location, T_precision> neg_binomial_2_log_lpmf(
   }
 
   for (size_t i = 0; i < max_size_seq_view; i++) {
-    if (include_summand<propto>::value) {
-      logp -= lgamma(n_vec[i] + 1.0);
+    if (phi__[i] > 1e10) {
+      // TODO(martinmodrak) This is wrong, but shouldn't brake most models.
+      // Will be adressed better in PR #1497
+      logp += poisson_log_lpmf(n_vec[i], eta__[i]);
+    } else {
+      if (include_summand<propto>::value) {
+        logp -= lgamma(n_vec[i] + 1.0);
+      }
+      if (include_summand<propto, T_precision>::value) {
+        logp += multiply_log(phi__[i], phi__[i]) - lgamma(phi__[i]);
+      }
+      if (include_summand<propto, T_log_location>::value) {
+        logp += n_vec[i] * eta__[i];
+      }
+      if (include_summand<propto, T_precision>::value) {
+        logp += lgamma(n_plus_phi[i]);
+      }
+      logp -= (n_plus_phi[i]) * logsumexp_eta_logphi[i];
     }
-    if (include_summand<propto, T_precision>::value) {
-      logp += multiply_log(phi__[i], phi__[i]) - lgamma(phi__[i]);
-    }
-    if (include_summand<propto, T_log_location>::value) {
-      logp += n_vec[i] * eta__[i];
-    }
-    if (include_summand<propto, T_precision>::value) {
-      logp += lgamma(n_plus_phi[i]);
-    }
-    logp -= (n_plus_phi[i]) * logsumexp_eta_logphi[i];
 
     if (!is_constant_all<T_log_location>::value) {
       ops_partials.edge1_.partials_[i]

diff --git a/stan/math/prim/prob/neg_binomial_2_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_lpmf.hpp
@@ -79,23 +79,25 @@ return_type_t<T_location, T_precision> neg_binomial_2_lpmf(
   }
 
   for (size_t i = 0; i < max_size_seq_view; i++) {
-    if (include_summand<propto>::value) {
-      logp -= lgamma(n_vec[i] + 1.0);
-    }
-    if (include_summand<propto, T_precision>::value) {
-      logp += multiply_log(phi__[i], phi__[i]) - lgamma(phi__[i]);
-    }
-    if (include_summand<propto, T_location>::value) {
-      logp += multiply_log(n_vec[i], mu__[i]);
-    }
-    if (include_summand<propto, T_precision>::value) {
-      logp += lgamma(n_plus_phi[i]);
-    }
-    logp -= (n_plus_phi[i]) * log_mu_plus_phi[i];
-
     // if phi is large we probably overflow, defer to Poisson:
-    if (phi__[i] > 1e5) {
-      logp = poisson_lpmf(n_vec[i], mu__[i]);
+    if (phi__[i] > 1e10) {
+      // TODO(martinmodrak) This is wrong, but shouldn't brake most models.
+      // Will be adressed better in PR #1497
+      logp += poisson_lpmf(n_vec[i], mu__[i]);
+    } else {
+      if (include_summand<propto>::value) {
+        logp -= lgamma(n_vec[i] + 1.0);
+      }
+      if (include_summand<propto, T_precision>::value) {
+        logp += multiply_log(phi__[i], phi__[i]) - lgamma(phi__[i]);
+      }
+      if (include_summand<propto, T_location>::value) {
+        logp += multiply_log(n_vec[i], mu__[i]);
+      }
+      if (include_summand<propto, T_precision>::value) {
+        logp += lgamma(n_plus_phi[i]);
+      }
+      logp -= (n_plus_phi[i]) * log_mu_plus_phi[i];
     }
 
     if (!is_constant_all<T_location>::value) {