All tests passed!

Improve bayes g lv examples

Author: Fontanapink

examples/gLV/examples-Rutter-Dekker.ipynb

@@ -279,7 +279,7 @@ def _run_inference_large(self, **kwargs) -> None:
             c2 = pm.InverseGamma("c2", 2, 8)
             tau = pm.HalfCauchy("tau", beta=tau0)
             lam = pm.HalfCauchy("lam", beta=1, shape=(ndim, ndim))
-            A = pm.Normal('A', mu=A_prior_mu, sigma=tau * lam *
+            A = pm.Normal('A', mu=A_prior_mu, sigma=tau * lam * \
                           at.sqrt(c2 / (c2 + tau**2 * lam**2)), shape=(ndim, ndim))
 
             # If noise covariance is provided, use it as a prior
@@ -438,14 +438,14 @@ def _run_inference_large_xs(self, **kwargs) -> None:
             c2_A = pm.InverseGamma("c2_A", 2, 1)
             tau_A = pm.HalfCauchy("tau_A", beta=tau0_A)
             lam_A = pm.HalfCauchy("lam_A", beta=1, shape=(nX, nX))
-            Ah = pm.Normal('Ah', mu=A_prior_mu, sigma=tau_A * lam_A *
+            Ah = pm.Normal('Ah', mu=A_prior_mu, sigma=tau_A * lam_A * \
                            at.sqrt(c2_A / (c2_A + tau_A**2 * lam_A**2)), shape=(nX, nX))
 
             tau0_B = (DB0 / (DB - DB0)) * 0.1 / np.sqrt(N)
             c2_B = pm.InverseGamma("c2_B", 2, 1)
             tau_B = pm.HalfCauchy("tau_B", beta=tau0_B)
             lam_B = pm.HalfCauchy("lam_B", beta=1, shape=(nS, nX))
-            Bh = pm.Normal('Bh', mu=0, sigma=tau_B * lam_B *
+            Bh = pm.Normal('Bh', mu=0, sigma=tau_B * lam_B * \
                            at.sqrt(c2_B / (c2_B + tau_B**2 * lam_B**2)), shape=(nS, nX))
 
             if noise_cov_prior is not None:

examples/gLV/examples-bayes-gLV.ipynb

@@ -2,27 +2,19 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pymc as pm
+import pandas as pd
+import seaborn as sns
 import pytensor.tensor as at
 import pickle
 import cloudpickle
+import os
+from typing import Optional, Union, List, Dict, Any
 
 from mimic.utilities import *
 from mimic.model_simulate.sim_gLV import *
 from mimic.model_infer.base_infer import BaseInfer
 
-from mimic.model_infer.base_infer import BaseInfer
-
-import os
-from typing import Optional, Union, List, Dict, Any
-
-
-import pandas as pd
-import numpy as np
-import seaborn as sns
-import matplotlib.pyplot as plt
-
 
-# Used in examples-Stein.ipynb
 def plot_params(mu_h, M_h, e_h, nsp):
     print("\ninferred params:")
     print("mu_hat/mu:")
@@ -69,8 +61,9 @@ def __init__(self,
                  prior_mu_sigma=None,
                  prior_Mii_mean=None,
                  prior_Mii_sigma=None,
-                 prior_Mij_sigma=None
-                 ):
+                 prior_Mij_sigma=None):
+
+        super().__init__()  # Call base class constructor
 
         # self.data = data  # data to do inference on
         self.X: Optional[np.ndarray] = X
@@ -217,7 +210,7 @@ def calculate_DA0(self, num_species, proportion=0.15):
         DA0 = int(round(expected_non_zero_elements))
         return max(DA0, 1)
 
-    def run_inference(self) -> None:
+    def run_inference(self, **kwargs) -> None:
         """
         This function infers the parameters for the Bayesian gLV model
 
@@ -314,9 +307,6 @@ def run_inference(self) -> None:
             # eg
             # print(f"mu_hat: {mu_hat.eval()}")
 
-            # initial_values = bayes_model.initial_point()
-            # print(f"Initial parameter values: {initial_values}")
-
             # Posterior distribution
             idata = pm.sample(
                 draws=draws,
@@ -327,7 +317,7 @@ def run_inference(self) -> None:
 
         return idata
 
-    def run_inference_shrinkage(self) -> None:
+    def run_inference_shrinkage(self, **kwargs) -> None:
         """
         This function infers the parameters for the Bayesian gLV model with Horseshoe prior for shrinkage
 
@@ -426,14 +416,15 @@ def run_inference_shrinkage(self) -> None:
             Y_obs = pm.Normal('Y_obs', mu=model_mean, sigma=sigma, observed=F)
 
             # For debugging:
+            # print if `debug` is set to 'high' or 'low'
+            if self.debug in ["high", "low"]:
+                initial_values = bayes_model.initial_point()
+                print(f"Initial parameter values: {initial_values}")
 
             # As tensor objects are symbolic, if needed print using .eval()
             # eg
             # print(f"mu_hat: {mu_hat.eval()}")
 
-            # initial_values = bayes_model.initial_point()
-            # print(f"Initial parameter values: {initial_values}")
-
             # Posterior distribution
             idata = pm.sample(
                 draws=draws,
@@ -443,7 +434,7 @@ def run_inference_shrinkage(self) -> None:
 
         return idata
 
-    def run_inference_shrinkage_pert(self) -> None:
+    def run_inference_shrinkage_pert(self, **kwargs) -> None:
         """
         This function infers the parameters for the Bayesian gLV model with Horseshoe prior for shrinkage
 
@@ -525,10 +516,8 @@ def run_inference_shrinkage_pert(self) -> None:
             lam = pm.HalfCauchy(
                 "lam", beta=1, shape=(
                     num_species, num_species - 1))
-            M_ij_hat = pm.Normal('M_ij_hat', mu=prior_Mij_sigma, sigma=tau * lam *
-                                 at.sqrt(c2 / (c2 + tau ** 2 * lam ** 2)),
-                                 shape=(num_species,
-                                        num_species - 1))
+            M_ij_hat = pm.Normal('M_ij_hat', mu=prior_Mij_sigma, sigma=tau * lam * at.sqrt(
+                c2 / (c2 + tau ** 2 * lam ** 2)), shape=(num_species, num_species - 1))
             # M_ij_hat = pm.Normal('M_ij_hat', mu=0, sigma=prior_Mij_sigma,
             # shape=(num_species, num_species - 1))  # different shape for
             # off-diagonal
@@ -553,14 +542,15 @@ def run_inference_shrinkage_pert(self) -> None:
             Y_obs = pm.Normal('Y_obs', mu=model_mean, sigma=sigma, observed=F)
 
             # For debugging:
+            # print if `debug` is set to 'high' or 'low'
+            if self.debug in ["high", "low"]:
+                initial_values = bayes_model.initial_point()
+                print(f"Initial parameter values: {initial_values}")
 
             # As tensor objects are symbolic, if needed print using .eval()
             # eg
             # print(f"mu_hat: {mu_hat.eval()}")
 
-            # initial_values = bayes_model.initial_point()
-            # print(f"Initial parameter values: {initial_values}")
-
             # Posterior distribution
             idata = pm.sample(
                 draws=draws,
@@ -587,17 +577,15 @@ def plot_posterior(self, idata):
         az.plot_posterior(
             idata,
             var_names=["mu_hat"],
-            ref_val=mu_hat_np.tolist()
-        )
+            ref_val=mu_hat_np.tolist())
         plt.savefig("plot-posterior-mu.pdf")
         plt.show()
         plt.close()
 
         az.plot_posterior(
             idata,
             var_names=["M_ii_hat"],
-            ref_val=np.diag(M_hat_np).tolist()
-        )
+            ref_val=np.diag(M_hat_np).tolist())
         plt.savefig("plot-posterior-Mii.pdf")
         plt.show()
         plt.close()
@@ -607,8 +595,7 @@ def plot_posterior(self, idata):
         az.plot_posterior(
             idata,
             var_names=["M_ij_hat"],
-            ref_val=M_ij.flatten().tolist()
-        )
+            ref_val=M_ij.flatten().tolist())
         plt.savefig("plot-posterior-Mij.pdf")
         plt.show()
         plt.close()
@@ -676,6 +663,11 @@ def plot_interaction_matrix(self, M, M_h):
                     color='white')
 
 
+######
+######
+######
+
+
 def param_data_compare(
         idata,
         F,
@@ -735,26 +727,24 @@ def curve_compare(idata, F, times, yobs, init_species_start, sim_gLV_class):
     # mu_h = idata.posterior['mu_hat'].mean(dim=('chain', 'draw')).values.flatten()
     # M_h= idata.posterior['M_hat'].mean(dim=('chain', 'draw')).values
 
-    predictor = sim_gLV(num_species=num_species,
-                        M=M_h.T,
-                        mu=mu_h
-                        )
+    predictor = sim_gLV(num_species=num_species, M=M_h.T, mu=mu_h)
     yobs_h, _, _, _, _ = predictor.simulate(
         times=times, init_species=init_species)
 
     plot_fit_gLV(yobs, yobs_h, times)
 
 
 def param_data_compare_pert(
-        idata,
-        F,
-        mu,
-        M,
-        epsilon,
-        num_perturbations,
-        times,
-        yobs,
-        init_species_start,
+    idata,
+    F,
+    mu,
+    M,
+    u,
+    epsilon,
+    num_perturbations,
+    times,
+    yobs,
+    init_species_start,
         sim_gLV_class):
     # az.to_netcdf(idata, 'model_posterior.nc')
     # Compare model parameters to the data
@@ -791,9 +781,9 @@ def param_data_compare_pert(
         epsilon=epsilon)
 
     yobs, init_species, mu, M, _ = simulator.simulate(
-        times=times, init_species=init_species, u=pert_fn)
+        times=times, init_species=init_species, u=u)
     yobs_h, _, _, _, _ = predictor.simulate(
-        times=times, init_species=init_species, u=pert_fn)
+        times=times, init_species=init_species, u=u)
 
     plot_fit_gLV(yobs, yobs_h, times)
 

mimic/model_infer/infer_VAR_bayes.py

@@ -20,11 +20,12 @@ def setup_data(request):
     X = np.random.randn(100, num_species+1)  # Random design matrix for 100 samples and 2 species
     F = np.random.randn(100, num_species)  # Random data matrix for 100 samples and 2 species
 
-    # Add a condition for different shapes
-    if 'test_run_inference' in request.node.nodeid or 'test_run_gLV_bayes_shrinkage' in request.node.nodeid:
-        X = np.random.randn(100, num_species + 1)  # Shape (100, num_species + 1)
-    elif 'test_run_bayes_gLV_shrinkage_pert' in request.node.nodeid or 'test_plot_posterior_pert' in request.node.nodeid:
-        X = np.random.randn(100, num_species + 2)  # Shape (100, num_species + 2)
+    # Add a condition for different shapes depending on number of parameters
+    test_name = request.node.nodeid.split("::")[-1] 
+    if test_name in ["test_run_inference", "test_run_inference_shrinkage"]:
+        X = np.random.randn(100, num_species + 1)
+    elif test_name in ["test_run_inference_shrinkage_pert", "test_plot_posterior_pert"]:
+        X = np.random.randn(100, num_species + 2)
 
 
     prior_mu_mean = 0
@@ -117,12 +118,12 @@ def test_run_inference(bayes_gLV_instance):
     assert len(idata.posterior["M_hat"]) > 0, "'M_hat' has no samples."
 
 
-def test_run_bayes_gLV_shrinkage(bayes_gLV_instance):
+def test_run_inference_shrinkage(bayes_gLV_instance):
     """
     Test the `run_bayes_gLV_shrinkage` function to check if it returns the correct output.
     """
     # Call the method to test
-    idata = bayes_gLV_instance.run_bayes_gLV_shrinkage()
+    idata = bayes_gLV_instance.run_inference_shrinkage()
 
     # Check that the output is an ArviZ InferenceData object
     assert isinstance(idata, az.InferenceData), "The output is not an InferenceData object."
@@ -132,12 +133,12 @@ def test_run_bayes_gLV_shrinkage(bayes_gLV_instance):
     assert "M_hat" in idata.posterior, "'M_hat' is not in the posterior."
 
 
-def test_run_bayes_gLV_shrinkage_pert(bayes_gLV_instance):
+def test_run_inference_shrinkage_pert(bayes_gLV_instance):
     """
     Test the `run_bayes_gLV_shrinkage_pert` function to check if it returns the correct output.
     """
     # Call the method to test
-    idata = bayes_gLV_instance.run_bayes_gLV_shrinkage_pert()
+    idata = bayes_gLV_instance.run_inference_shrinkage_pert()
 
     # Check that the output is an ArviZ InferenceData object
     assert isinstance(idata, az.InferenceData), "The output is not an InferenceData object."
@@ -153,7 +154,7 @@ def test_plot_posterior(bayes_gLV_instance):
     Test the `plot_posterior` function to ensure it runs without errors.
     """
     # Generate the InferenceData from the method
-    idata = bayes_gLV_instance.run_bayes_gLV_shrinkage()
+    idata = bayes_gLV_instance.run_inference_shrinkage()
 
     # Try to call the plot_posterior method
     try:
@@ -167,7 +168,7 @@ def test_plot_posterior_pert(bayes_gLV_instance):
     Test the `plot_posterior_pert` function to ensure it produces the correct plot.
     """
     # First, generate the InferenceData from the method
-    idata = bayes_gLV_instance.run_bayes_gLV_shrinkage_pert()
+    idata = bayes_gLV_instance.run_inference_shrinkage_pert()
 
     # Try to call the plot_posterior_pert method
     try: