Implement Poisson GLM model in PyMC

posterior_database/data/data/GLM_Poisson_Data.json

@@ -0,0 +1,5 @@
+{
+  "year": [-1.66802789939819, -1.58248800712136, -1.49694811484453, -1.4114082225677, -1.32586833029087, -1.24032843801404, -1.15478854573721, -1.06924865346038, -0.983708761183547, -0.898168868906717, -0.812628976629887, -0.727089084353056, -0.641549192076226, -0.556009299799396, -0.470469407522566, -0.384929515245736, -0.299389622968906, -0.213849730692075, -0.128309838415245, -0.0427699461384151, 0.0427699461384151, 0.128309838415245, 0.213849730692075, 0.299389622968906, 0.384929515245736, 0.470469407522566, 0.556009299799396, 0.641549192076226, 0.727089084353056, 0.812628976629887, 0.898168868906717, 0.983708761183547, 1.06924865346038, 1.15478854573721, 1.24032843801404, 1.32586833029087, 1.4114082225677, 1.49694811484453, 1.58248800712136, 1.66802789939819],
+  "C": [29, 36, 19, 28, 36, 29, 20, 19, 35, 32, 34, 34, 33, 47, 48, 46, 46, 49, 60, 64, 87, 85, 85, 95, 115, 127, 142, 168, 181, 194, 200, 210, 208, 235, 244, 272, 239, 263, 239, 245],
+  "n": 40
+}

posterior_database/models/info/GLM_Poisson_model.info.json

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

posterior_database/models/pymc/GLM_Poisson_model.py

@@ -0,0 +1,32 @@
+import numpy as np
+import pymc as pm
+
+
+def model(data):
+    years_data = np.array(data["year"])
+    counts_data = np.array(data["C"])
+    number_observations = data["n"]
+    coords = {
+        "features": ["years", "years**2", "years**3"],
+        "observation": np.arange(number_observations),
+    }
+    with pm.Model(coords=coords) as pymc_model:
+        alpha = pm.Uniform("alpha", -20, +20)
+        beta = pm.Uniform("beta", -10, +10, dims="features")
+        X = pm.Data(
+            "X",
+            np.column_stack([years_data, years_data**2, years_data**3]),
+            dims=["observation", "feature"],
+        )
+        y = pm.Data("y", counts_data, dims="observation")
+
+        log_lambda = pm.Deterministic(
+            "log_lambda", alpha + X @ beta, dims="observation"
+        )
+
+        counts = pm.Poisson(
+            "counts", mu=pm.math.exp(log_lambda), observed=y, dims="observation"
+        )
+
+    return pymc_model
+