Develop metatree poisson

bayesml/poisson/_poisson.py

@@ -10,6 +10,7 @@
 from scipy.stats import poisson as ss_poisson
 from scipy.stats import gamma as ss_gamma 
 from scipy.stats import nbinom as ss_nbinom
+from scipy.special import gammaln
 import matplotlib.pyplot as plt
 
 from .. import base
@@ -218,6 +219,9 @@ def __init__(self,h0_alpha=1.0,h0_beta=1.0):
         self.p_r = 1.0
         self.p_theta = 0.5
 
+        #statistics
+        self._sum_log_factorial = 0.0
+
         self.set_h0_params(h0_alpha,h0_beta)
 
     def get_constants(self):
@@ -270,6 +274,7 @@ def set_hn_params(self,hn_alpha=None,hn_beta=None):
         hn_beta : float, optional
             a positive real number, by default None
         """
+        self._sum_log_factorial = 0.0
         if hn_alpha is not None:
             self.hn_alpha = _check.pos_float(hn_alpha,'hn_alpha',ParameterFormatError)
         if hn_beta is not None:
@@ -305,6 +310,7 @@ def update_posterior(self,x):
             self.hn_beta += x.size
         except:
             self.hn_beta += 1
+        self._sum_log_factorial += gammaln(x+1).sum()
         return self
 
     def _update_posterior(self,x):
@@ -479,3 +485,17 @@ def pred_and_update(self,x,loss="squared"):
         prediction = self.make_prediction(loss=loss)
         self.update_posterior(x)
         return prediction
+
+    def calc_log_marginal_likelihood(self):
+        """Calculate log marginal likelihood
+
+        Returns
+        -------
+        log_marginal_likelihood : float
+            The log marginal likelihood.
+        """
+        return (self.h0_alpha * np.log(self.h0_beta)
+                - gammaln(self.h0_alpha)
+                - self.hn_alpha * np.log(self.hn_beta)
+                + gammaln(self.hn_alpha)
+                - self._sum_log_factorial)