diff --git a/posterior_database/models/info/eight_schools_noncentered.info.json b/posterior_database/models/info/eight_schools_noncentered.info.json
index b83367f9..9799a1eb 100644
--- a/posterior_database/models/info/eight_schools_noncentered.info.json
+++ b/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",
diff --git a/posterior_database/models/pymc/eight_schools_noncentered.py b/posterior_database/models/pymc/eight_schools_noncentered.py
new file mode 100644
index 00000000..b7cc0e27
--- /dev/null
+++ b/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
diff --git a/posterior_database/models/pymc3/eight_schools_noncentered.py b/posterior_database/models/pymc3/eight_schools_noncentered.py
deleted file mode 100644
index 91b6c611..00000000
--- a/posterior_database/models/pymc3/eight_schools_noncentered.py
+++ /dev/null
@@ -1,18 +0,0 @@
-import numpy as np
-import pymc3 as pm3
-
-
-def model(data):
-    J = data["J"]  # number of schools
-    y_obs = np.array(data["y"])  # estimated treatment
-    sigma = np.array(data["sigma"])  # std of estimated effect
-    with pm3.Model() as pymc_model:
-
-        mu = pm3.Normal(
-            "mu", mu=0, sd=5
-        )  # hyper-parameter of mean, non-informative prior
-        tau = pm3.Cauchy("tau", alpha=0, beta=5)  # hyper-parameter of sd
-        theta_trans = pm3.Normal("theta_trans", mu=0, sd=1, shape=J)
-        theta = mu + tau * theta_trans
-        y = pm3.Normal("y", mu=theta, sd=sigma, observed=y_obs)
-    return pymc_model