Implement Binomial GLM model in PyMC

.gitignore

@@ -16,3 +16,4 @@ TODO.tdl
 *.hpp
 *.dll
 tmp*
+.ipynb_checkpoints/*

CITATION.cff

@@ -8,5 +8,5 @@ authors:
   - family-names: Vehtari
     given-names: Aki
 title: "posteriordb: a set of posteriors for Bayesian inference and probabilistic programming"
-version: 0.4
+version: 0.5
 date-released: 2023-10-15

doc/sampling_pymc_models.md

@@ -0,0 +1,30 @@
+## Sampling PyMC models
+
+Steps to sample PyMC models: 
+
+1. load the model json data from the respctive zipfile as a Python dictionary. 
+
+For example, for the GLM binomial model:
+
+```python
+import json
+import zipfile
+
+# Load zipfile
+zfile = zipfile.ZipFile(r"posterior_database/data/data/GLM_Binomial_data.json.zip", 'r')
+data = json.loads(zfile.read("GLM_Binomial_data.json"))
+```
+
+2. Call the model defining function, with the data dictionary and, in a `with` context, call `pm.sample`:
+
+```python
+import pymc as pm
+
+from posterior_database.models.pymc.GLM_Binomial_model import model
+
+# Call the function that creates the PyMC model and sample it
+with model(data):
+    trace = pm.sample()
+
+print(pm.stats.summary(trace))
+```

posterior_database/models/info/GLM_Binomial_model.info.json

@@ -11,6 +11,10 @@
     "stan": {
       "model_code": "models/stan/GLM_Binomial_model.stan",
       "stan_version": ">=2.26.0"
+    },
+    "pymc": {
+        "model_code": "models/pymc/GLM_Binomial_model.py",
+        "pymc_version": ">=5.16.2"    
     }
   }
 }

posterior_database/models/pymc/GLM_Binomial_model.py

@@ -0,0 +1,21 @@
+import numpy as np
+import pymc as pm
+
+def model(data):
+
+    with pm.Model(coords={"year": range(40)}) as pymc_model:    
+        year = pm.Data("year_data", np.array(data["year"]), dims="year")
+        N = pm.Data("N_data", np.array(data["N"]), dims="year")
+        C = pm.Data("C_data", np.array(data["C"]), dims="year")
+
+        alpha = pm.Normal("alpha", mu=0, sigma=100)
+        beta1 = pm.Normal("beta1", mu=0, sigma=100)
+        beta2 = pm.Normal("beta2", mu=0, sigma=100)
+
+        year_squared = year ** 2
+        logit_p = alpha + beta1 * year + beta2 * year_squared        
+        pm.Binomial("C", n=N, logit_p=logit_p, observed=C, dims="year")
+
+        p = pm.Deterministic("p", pm.math.invlogit(logit_p), dims="year")
+
+    return pymc_model