Development

CITATION.cff

@@ -3,10 +3,16 @@ message: "If you use posteriordb, please cite it as below."
 authors:
   - family-names: Magnusson
     given-names: Måns
+  - family-names: Torgander
+    given-names: Jakob
   - family-names: Bürkner
     given-names: Paul
+  - family-names: Zhang
+    given-names: Lu
+  - family-names: Carpenter
+    given-names: Bob
   - family-names: Vehtari
     given-names: Aki
-title: "posteriordb: a set of posteriors for Bayesian inference and probabilistic programming"
-version: 0.5
-date-released: 2023-10-15
+title: "posteriordb: Testing, Benchmarking and Developing Bayesian Inference Algorithms"
+version: 1.0
+date-released: 2025-04-29

README.Rmd

@@ -45,7 +45,15 @@ Most model code are using the the same BSD-3 licence. Although, some models and
 To simplify the use of `posteriordb`, there are convenience functions both in Python and in R. To use R, see the [posteriordb-r](https://github.com/stan-dev/posteriordb-r) repository, and to use Python, see the [posteriordb-python](https://github.com/stan-dev/posteriordb-python) repository.
 
 ## Citing `posteriordb`
-Developing and maintaining open-source software is an important yet often underappreciated contribution to scientific progress. Thus, please make sure to cite it appropriately so that developers get credit for their work. Information on how to cite `posteriordb` can be found in the [CITATION.cff](https://github.com/stan-dev/posteriordb/blob/master/CITATION.cff) file. Use the "cite this repository" button under "About" to get a simple BibTeX or APA snippet.
+Developing and maintaining open-source software is an important yet often underappreciated contribution to scientific progress. Thus, please make sure to cite it appropriately so that developers get credit for their work. 
+
+To cite the actual posterior database on github, information can be found in the [CITATION.cff](https://github.com/stan-dev/posteriordb/blob/master/CITATION.cff) file. Use the "cite this repository" button under "About" to get a simple BibTeX or APA snippet.
+
+We have also released a paper describing `posteriordb` that we encourage to cite when using `posteriordb`:
+
+Magnusson, M., Torgander, J., Bürkner, P., Zhang, L., Carpenter, B., Vehtari, A.. (2025). posteriordb: Testing, Benchmarking and Developing Bayesian Inference Algorithms. *Proceedings of The 28th International Conference on Artificial Intelligence and Statistics*, in *Proceedings of Machine Learning Research* 258:1198-1206 Available from (https://proceedings.mlr.press/v258/magnusson25a.html)[https://proceedings.mlr.press/v258/magnusson25a.html].
+
+If you do extensive testing of individual posterior, also please consider citing the paper that introduced the posterior by citing the origin of the posterior. This information can be found in the meta data associated with the posterior, data an/or model.
 
 As `posteriordb` rely heavily on Stan, so please consider also to cite Stan:
 

README.md

@@ -73,11 +73,29 @@ repository.
 Developing and maintaining open-source software is an important yet
 often underappreciated contribution to scientific progress. Thus, please
 make sure to cite it appropriately so that developers get credit for
-their work. Information on how to cite `posteriordb` can be found in the
+their work.
+
+To cite the actual posterior database on github, information can be
+found in the
 [CITATION.cff](https://github.com/stan-dev/posteriordb/blob/master/CITATION.cff)
 file. Use the “cite this repository” button under “About” to get a
 simple BibTeX or APA snippet.
 
+We have also released a paper describing `posteriordb` that we encourage
+to cite when using `posteriordb`:
+
+Magnusson, M., Torgander, J., Bürkner, P., Zhang, L., Carpenter, B.,
+Vehtari, A.. (2025). posteriordb: Testing, Benchmarking and Developing
+Bayesian Inference Algorithms. *Proceedings of The 28th International
+Conference on Artificial Intelligence and Statistics*, in *Proceedings
+of Machine Learning Research* 258:1198-1206 Available from
+(<https://proceedings.mlr.press/v258/magnusson25a.html>)\[<https://proceedings.mlr.press/v258/magnusson25a.html>\].
+
+If you do extensive testing of individual posterior, also please
+consider citing the paper that introduced the posterior by citing the
+origin of the posterior. This information can be found in the meta data
+associated with the posterior, data an/or model.
+
 As `posteriordb` rely heavily on Stan, so please consider also to cite
 Stan:
 

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"    
     }
   },
   "licence": "BSD3"

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",
+      "pymc_version": ">=5.16.2"
     }
   },
   "licence": "BSD3"

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/info/earn_height.info.json

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

posterior_database/models/info/eight_schools_noncentered.info.json

@@ -13,8 +13,9 @@
       "model_code": "models/stan/eight_schools_noncentered.stan",
       "stan_version": ">=2.26.0"
     },
-    "pymc3": {
-      "model_code": "models/pymc3/eight_schools_noncentered.py"
+    "pymc": {
+      "model_code": "models/pymc/eight_schools_noncentered.py",
+      "pymc_version": ">=5.16.2" 
     }
   },
   "added_by": "Mans Magnusson",

posterior_database/models/info/rats_model.info.json

@@ -1,7 +1,7 @@
 {
   "name": "rats_model",
-  "keywords": ["Heirarchical", "Normal", "Rats", "Random Effects", "Linear Growth", "Linear"],
-  "title": "Normal Heirarchical Model to Model Rats' Weight Over Time",
+  "keywords": ["Hierarchical", "Normal", "Rats", "Random Effects", "Linear Growth", "Linear"],
+  "title": "Normal Hierarchical Model to Model Rats' Weight Over Time",
   "description": "The model is essentially a random effects linear growth curve.",
   "urls": "https://github.com/stan-dev/example-models/tree/master/bugs_examples/vol1/rats",
   "references": "gelfand1990illustration",

posterior_database/models/info/seeds_centered_model.info.json

@@ -1,7 +1,7 @@
 {
   "name": "seeds_centered_model",
   "keywords": ["Random", "Effect", "Logistic", "Regression", "Stanified", "Centered"],
-  "title": "Normal Heirarchical Model to Model Rats' Weight Over Time",
+  "title": "Random Effect Logistic Regression for Seed Germination Proportion",
   "description": "The model is essentially a random effects logistic,\n          allowing for over-dispersion. Models the proportion of germination\n          of different seeds in different root extracts. 'Stanified' & centered:\n          tau replaced by sigma, direct estimation using \n          narrower semi-informative priors, coefficients are centered",
   "urls": "http://www.mrc-bsu.cam.ac.uk/wp-content/uploads/WinBUGS_Vol1.pdf",
   "references": "crowder1978beta",

posterior_database/models/info/seeds_stanified_model.info.json

@@ -1,7 +1,7 @@
 {
   "name": "seeds_stanified_model",
   "keywords": ["Random", "Effect", "Logistic", "Regression", "Stanified"],
-  "title": "Normal Heirarchical Model to Model Rats' Weight Over Time",
+  "title": "Random Effect Logistic Regression for Seed Germination Proportion",
   "description": "The model is essentially a random effects logistic,\n          allowing for over-dispersion. Models the proportion of germination\n          of different seeds in different root extracts. 'Stanified':\n          tau replaced by sigma, direct estimation using \n          narrower semi-informative priors",
   "urls": "http://www.mrc-bsu.cam.ac.uk/wp-content/uploads/WinBUGS_Vol1.pdf",
   "references": "crowder1978beta",

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

posterior_database/models/pymc/GLM_Poisson_model.py

@@ -0,0 +1,36 @@
+import numpy as np
+import pymc as pm
+import pytensor.tensor as pt
+
+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)
+        beta1 = pm.Uniform("beta1", -10, +10)
+        beta2 = pm.Uniform("beta2", -10, +10)
+        beta3 = pm.Uniform("beta3", -10, +10)
+
+        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")
+
+        beta = pm.Deterministic("beta", pt.stack([beta1, beta2, beta3]), dims = "feature")
+
+        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

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/models/pymc/arma11.py

@@ -0,0 +1,43 @@
+import pymc as pm
+import pytensor
+import pytensor.tensor as pt
+def model(data):
+    y_obs = data["y"]  # Observed series
+    up_to = data["T"]  # Timesteps
+    with pm.Model() as pymc_model:
+        mu = pm.Normal(
+            "mu", mu=0, sigma=10
+        )  # mean coefficient
+        phi = pm.Normal(
+            "phi", mu=0, sigma=2
+        ) # autoregressive coef
+        theta = pm.Normal(
+            "theta", mu=0, sigma=2
+        ) # moving average
+        sigma = pm.HalfCauchy(
+            "sigma", beta=2.5
+        ) # noise scale
+
+        # scan variables
+        y = pt.as_tensor_variable(y_obs)
+        y_hat0 = mu + phi * mu
+        err0 = y[0] - y_hat0
+        outputs_info = [y_hat0, err0, y[0]]
+
+        def step(curr_obs, prev_pred, prev_err, prev_obs, mu, phi, theta):
+            y_hat = mu + phi * prev_obs + theta * prev_err
+            new_err = curr_obs - y_hat
+            return [y_hat, new_err, curr_obs]
+
+        [predictions, _, _] , _ = pytensor.scan(fn=step,
+                                              outputs_info=outputs_info,    
+                                              sequences=y[1:],              
+                                              non_sequences=[mu, phi, theta])
+
+        # concatenating predictions
+        y_hat0_vec = y_hat0.dimshuffle("x")
+        final_predictions = pt.concatenate([y_hat0_vec, predictions])
+
+        pm.Normal("y", mu=final_predictions, sigma=sigma, observed=y)
+
+    return pymc_model

posterior_database/models/pymc/earn_height.py

@@ -0,0 +1,23 @@
+import numpy as np
+import pymc as pm
+
+def model(data):
+
+    # Define PyMC Model
+    coords = {
+        'obs_idx': np.arange(data["N"]),
+        'feature': ['intercept', 'height']
+    }
+
+    with pm.Model(coords=coords) as earn_height:
+        N = pm.Data("N", data["N"])
+        earn = pm.Data('earn', data['earn'], dims=['obs_idx'])
+        height = pm.Data('height', data['height'], dims=['obs_idx'])
+
+        beta = pm.Flat('beta', dims=['feature'])
+        mu = beta[0] + beta[1] * height
+
+        sigma = pm.HalfFlat('sigma')
+        earn_hat = pm.Normal('earn_hat', mu=mu, sigma=sigma, observed=earn, dims=['obs_idx'])
+
+    return earn_height

posterior_database/models/pymc/eight_schools_noncentered.py

@@ -0,0 +1,18 @@
+import numpy as np
+import pymc as pm
+
+
+def model(data):
+    y_obs = np.array(data["y"])  # estimated treatment
+    sigma = np.array(data["sigma"])  # std of estimated effect
+    coords = {"school": np.arange(data["J"])}
+    with pm.Model(coords=coords) as pymc_model:
+
+        mu = pm.Normal(
+            "mu", mu=0, sigma=5
+        )  # hyper-parameter of mean, non-informative prior
+        tau = pm.HalfCauchy("tau", beta=5)  # hyper-parameter of sigma
+        theta_trans = pm.Normal("theta_trans", mu=0, sigma=1, dims="school")
+        theta = mu + tau * theta_trans
+        y = pm.Normal("y", mu=theta, sigma=sigma, observed=y_obs)
+    return pymc_model