Merge pull request #266 from fonnesbeck/M0_pymc

PyMC implementation of M0 model

Author: JTorgander

posterior_database/models/info/M0_model.info.json

@@ -11,6 +11,10 @@
     "stan": {
       "model_code": "models/stan/M0_model.stan",
       "stan_version": ">=2.26.0"
+    },
+    "pymc": {
+      "model_code": "models/pymc/M0_model.py",
+      "pymc_version": ">=5.16.0"
     }
   },
   "licence": "BSD3"

posterior_database/models/info/M0_model_discrete.info.json

@@ -0,0 +1,16 @@
+{
+  "name": "M0_model_discrete",
+  "keywords": ["BPA", "Ch.6", "Population", "Capture", "Recapture", "Individual", "Constant", "Discrete"],
+  "title": "Inferring population size with constant detection probability",
+  "description": "Detection probability of an individual is assumed constant over individuals and over time periods",
+  "urls": "https://github.com/stan-dev/example-models/blob/master/BPA/Ch.06",
+  "references": "kery2011population",
+  "added_by": "Chris Fonnesbeck",
+  "added_date": "2021-06-24",
+  "model_implementations": {
+    "pymc": {
+      "model_code": "models/pymc/M0_model_discrete.py",
+      "pymc_version": ">=5.16.0"
+    }
+  }
+}

posterior_database/models/pymc/M0_model.py

@@ -0,0 +1,35 @@
+import numpy as np
+import pymc as pm
+
+def model(data):
+    y = np.array(data["y"])  # capture history matrix
+    T = np.array(data["T"])  # time periods
+ 
+    coords = {"individual": np.arange(data["M"]),
+              "capture_period": np.arange(data["T"])}
+    with pm.Model(coords=coords) as pymc_model:
+        y_data = pm.Data("y", y, dims=("individual", "capture_period"))
+        s = pm.Deterministic("s", y_data.sum(axis=1), dims="individual")
+        is_observed = s > 0
+    
+        # Inclusion probability
+        omega = pm.Uniform("omega", 0, 1)
+        # Capture probability
+        p = pm.Uniform("p", 0, 1) 
+
+        # Defining bernoulli and binomial components
+        binom = pm.Binomial.dist(n=T, p=p)
+        log_omega = pm.math.log(omega)
+        log_one_minus_omega = pm.math.log(1-omega)
+        log_one_minus_p = pm.math.log(1-p)
+
+        # Computing marginalization mixture
+        logp_if_obs = log_omega + pm.logp(binom, s)
+        logp_zero_single = pm.math.logaddexp(log_omega + T * log_one_minus_p , log_one_minus_omega)
+        logp_each = pm.math.switch(is_observed, logp_if_obs, logp_zero_single)
+
+        pm.Potential("likelihood", logp_each.sum())
+        
+        omega_nd = pm.Deterministic("omega_nd", (omega * (1 -p)**T)  / (omega * (1 - p)**T  + (1 - omega)))
+    
+    return pymc_model

posterior_database/models/pymc/M0_model_discrete.py

@@ -0,0 +1,25 @@
+import numpy as np
+import pymc as pm
+
+def model(data):
+    y = np.array(data["y"])  # capture history matrix
+    T = np.array(data["T"])  # time periods
+    coords = {"individual": np.arange(data["M"]),
+              "capture_period": np.arange(data["T"])}
+    with pm.Model(coords=coords) as pymc_model:
+
+        # Inclusion probability
+        omega = pm.Uniform("omega", 0, 1)
+        # Capture probability
+        p = pm.Uniform("p", 0, 1) 
+
+        # Inclusion indicator
+        z = pm.Bernoulli("z", p=omega, dims="individual")
+    
+        y_obs = pm.Bernoulli("y_obs", p=z[:, None]*p, observed=y, dims=("individual", "capture_period"))
+
+        N = pm.Deterministic("N", pm.math.sum(z))
+
+        omega_nd = pm.Deterministic("omega_nd", (omega * (1 -p)**T)  / (omega * (1 - p)**T  + (1 - omega)))
+
+    return pymc_model

posterior_database/posteriors/M0_data-M0_model_discrete.json

@@ -0,0 +1,18 @@
+{
+  "keywords": ["bpa book", "capture-recapture"],
+  "urls": "https://github.com/stan-dev/example-models/blob/master/BPA/Ch.06",
+  "references": "kery2011population",
+  "dimensions": {
+    "omega": 1,
+    "p": 1,
+    "omega_nd": 1,
+    "N": 1,
+    "z": 237
+  },
+  "reference_posterior_name": null,
+  "added_date": "2026-01-16",
+  "added_by": "Chris Fonnesbeck",
+  "name": "M0_data-M0_model_discrete",
+  "model_name": "M0_model_discrete",
+  "data_name": "M0_data"
+}