Develop metatree

bayesml/bernoulli/_bernoulli.py

@@ -5,7 +5,7 @@
 import warnings
 import numpy as np
 from scipy.stats import beta as ss_beta
-# from scipy.stats import betabinom as ss_betabinom
+from scipy.special import gammaln
 import matplotlib.pyplot as plt
 
 from .. import base
@@ -486,3 +486,18 @@ 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 (gammaln(self.h0_alpha+self.h0_beta)
+                -gammaln(self.h0_alpha)
+                -gammaln(self.h0_beta)
+                -gammaln(self.hn_alpha+self.hn_beta)
+                +gammaln(self.hn_alpha)
+                +gammaln(self.hn_beta))

bayesml/categorical/_categorical.py

@@ -9,6 +9,7 @@
 import warnings
 import numpy as np
 from scipy.stats import dirichlet as ss_dirichlet
+from scipy.special import gammaln
 import matplotlib.pyplot as plt
 
 from .. import base
@@ -553,3 +554,16 @@ def pred_and_update(self, x, loss="squared",onehot=True):
         prediction = self.make_prediction(loss,onehot)
         self.update_posterior(x,onehot)
         return prediction
+
+    def calc_log_marginal_likelihood(self):
+        """Calculate log marginal likelihood
+
+        Returns
+        -------
+        log_marginal_likelihood : float
+            The log marginal likelihood.
+        """
+        return (gammaln(self.h0_alpha_vec.sum())
+                -gammaln(self.h0_alpha_vec).sum()
+                -gammaln(self.hn_alpha_vec.sum())
+                +gammaln(self.hn_alpha_vec).sum())

bayesml/exponential/_exponential.py

@@ -9,6 +9,7 @@
 import warnings
 import numpy as np
 from scipy.stats import expon as ss_expon, gamma as ss_gamma, lomax as ss_lomax
+from scipy.special import gammaln
 import matplotlib.pyplot as plt
 
 from .. import base
@@ -481,3 +482,16 @@ 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))

bayesml/linearregression/_linearregression.py

@@ -9,6 +9,7 @@
 from scipy.stats import gamma as ss_gamma
 from scipy.stats import multivariate_t as ss_multivariate_t
 from scipy.stats import t as ss_t
+from scipy.special import gammaln
 import matplotlib.pyplot as plt
 
 from .. import base
@@ -360,6 +361,9 @@ def __init__(
         self.p_lambda =  1.0
         self.p_nu = 2.0
 
+        # sample size
+        self._n = 0
+
         self.set_h0_params(
             h0_mu_vec,
             h0_lambda_mat,
@@ -462,6 +466,7 @@ def set_hn_params(
         hn_beta : float, optional
             a positive real number, by default None.
         """
+        self._n = 0
         if hn_mu_vec is not None:
             _check.float_vec(hn_mu_vec,'hn_mu_vec',ParameterFormatError)
             _check.shape_consistency(
@@ -539,6 +544,7 @@ def update_posterior(self, x, y):
         self.hn_alpha +=   x.shape[0]/2.0
         self.hn_beta += (-self.hn_mu_vec[np.newaxis,:] @ self.hn_lambda_mat @ self.hn_mu_vec[:,np.newaxis]
                          + y @ y + hn1_mu[np.newaxis,:] @ hn1_Lambda @ hn1_mu[:,np.newaxis])[0,0] /2.0
+        self._n += x.shape[0]
         return self
 
     def _update_posterior(self, x, y):
@@ -552,6 +558,7 @@ def _update_posterior(self, x, y):
         self.hn_alpha +=   x.shape[0]/2.0
         self.hn_beta += (-self.hn_mu_vec[np.newaxis,:] @ self.hn_lambda_mat @ self.hn_mu_vec[:,np.newaxis]
                          + y @ y + hn1_mu[np.newaxis,:] @ hn1_Lambda @ hn1_mu[:,np.newaxis])[0,0] /2.0
+        self._n += x.shape[0]
         return self
 
     def estimate_params(self,loss="squared",dict_out=False):
@@ -775,3 +782,23 @@ def pred_and_update(self,x,y,loss="squared"):
         prediction = self.make_prediction(loss=loss)
         self.update_posterior(x,y)
         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)
+            - self.hn_alpha * np.log(self.hn_beta)
+            + gammaln(self.hn_alpha)
+            - gammaln(self.h0_alpha)
+            + 0.5 * (
+                np.linalg.slogdet(self.h0_lambda_mat)[1]
+                - np.linalg.slogdet(self.hn_lambda_mat)[1]
+                - self._n * np.log(2*np.pi)
+            )
+        )