fix error in ode notebook

biosustain · Apr 29, 2024 · 38adc5f · 38adc5f
1 parent fd034f5
commit 38adc5f
Show file tree

Hide file tree

Showing 22 changed files with 7,516 additions and 7,563 deletions.
diff --git a/docs/after_mcmc.html b/docs/after_mcmc.html
diff --git a/docs/after_mcmc.ipynb b/docs/after_mcmc.ipynb
diff --git a/docs/hierarchical_models.html b/docs/hierarchical_models.html
diff --git a/docs/hierarchical_models.ipynb b/docs/hierarchical_models.ipynb
diff --git a/docs/img/monod_posteriors.png b/docs/img/monod_posteriors.png
diff --git a/docs/img/monod_priors.png b/docs/img/monod_priors.png
diff --git a/docs/img/monod_simulated_data.png b/docs/img/monod_simulated_data.png
diff --git a/docs/index.ipynb b/docs/index.ipynb
@@ -205,7 +205,7 @@
         "Discussion).” *Bayesian Analysis* 16 (2): 667–718.\n",
         "<https://doi.org/10.1214/20-BA1221>."
       ],
-      "id": "e17fdecc-1fc7-4380-9ddc-8efa49cdb671"
+      "id": "3de702ca-c484-4274-b592-38e972460b2a"
     }
   ],
   "nbformat": 4,

diff --git a/docs/introduction_to_bayesian_inference.ipynb b/docs/introduction_to_bayesian_inference.ipynb
diff --git a/docs/mcmc_and_stan.ipynb b/docs/mcmc_and_stan.ipynb
diff --git a/docs/metropolis-hastings.html b/docs/metropolis-hastings.html
diff --git a/docs/metropolis-hastings.ipynb b/docs/metropolis-hastings.ipynb
diff --git a/docs/odes.html b/docs/odes.html
diff --git a/docs/odes.ipynb b/docs/odes.ipynb
diff --git a/docs/regression.html b/docs/regression.html
diff --git a/docs/regression.ipynb b/docs/regression.ipynb
diff --git a/docs/search.json b/docs/search.json
diff --git a/materials/img/monod_posteriors.png b/materials/img/monod_posteriors.png
diff --git a/materials/img/monod_priors.png b/materials/img/monod_priors.png
diff --git a/materials/img/monod_simulated_data.png b/materials/img/monod_simulated_data.png
diff --git a/materials/odes.qmd b/materials/odes.qmd
@@ -124,8 +124,10 @@ true_param_values = {
     "t_ks": 0.3,
     "t_gamma": 0.13,
     "species_zero": [
-        [np.exp(np.random.normal(-2.1, 0.1)), np.exp(np.random.normal(0.2, 0.1))]
-        for _ in range(N_tube)
+        [
+            np.exp(np.random.normal(-2.1, 0.05)), 
+            np.exp(np.random.normal(0.2, 0.05))
+        ] for _ in range(N_tube)
     ],
     "sigma_y": [0.08, 0.1],
     "ln_mu_max_z": np.random.normal(0, 1, size=N_strain).tolist(),
@@ -344,30 +346,22 @@ idata_prior
 We can find the prior intervals for the true abundance and plot them in the graph.
 
 ```{python}
-prior_abundances = (
-    idata_prior
-    .posterior["abundance"]
-    .to_dataframe()
-    .reset_index()
-    .join(timepoints, on="timepoint")
-)
+prior_abundances = idata_prior.posterior["abundance"]
+
 n_sample = 20
-chains_and_draws = rng.choice(prior_abundances[["chain", "draw"]].values, 10)
+chains = rng.choice(prior_abundances.coords["chain"].values, n_sample)
+draws = rng.choice(prior_abundances.coords["draw"].values, n_sample)
 f, axes = plot_sim(true_abundance, fake_measurements, species_to_ax)
-for tube_i in tubes:
-    for species_j in species:
-        abundance_sample = (
-            prior_abundances.loc[
-                lambda df: (df["tube"] == tube_i) & (df["species"] == species_j)
-            ]
-            .set_index(["chain", "draw"])
-            .loc[chains_and_draws.tolist()]
-            .reset_index()
-        )
-        axes[species_to_ax[species_j]].plot(
-            abundance_sample.set_index("time")["abundance"],
-            alpha=0.5, color="skyblue", zorder=-1
-        )
+
+for ax, species_i in zip(axes, species):
+    for tube_j in tubes:
+        for chain, draw in zip(chains, draws):
+            timeseries = prior_abundances.sel(chain=chain, draw=draw, tube=tube_j, species=species_i)
+            ax.plot(
+                timepoints.values, 
+                timeseries.values,
+                alpha=0.5, color="skyblue", zorder=-1
+            )
 f.savefig("img/monod_priors.png")
 ```
 
@@ -414,32 +408,22 @@ az.summary(idata_posterior)
 Show posterior intervals
 
 ```{python}
-posterior_abundances = (
-    idata_posterior
-    .posterior["abundance"]
-    .to_dataframe()
-    .reset_index()
-    .join(timepoints, on="timepoint")
-)
+prior_abundances = idata_posterior.posterior["abundance"]
+
 n_sample = 20
-chains_and_draws = rng.choice(
-    posterior_abundances[["chain", "draw"]].values, 10
-)
+chains = rng.choice(prior_abundances.coords["chain"].values, n_sample)
+draws = rng.choice(prior_abundances.coords["draw"].values, n_sample)
 f, axes = plot_sim(true_abundance, fake_measurements, species_to_ax)
-for tube_i in tubes:
-    for species_j in species:
-        abundance_sample = (
-            posterior_abundances.loc[
-                lambda df: (df["tube"] == tube_i) & (df["species"] == species_j)
-            ]
-            .set_index(["chain", "draw"])
-            .loc[chains_and_draws.tolist()]
-            .reset_index()
-        )
-        axes[species_to_ax[species_j]].plot(
-            abundance_sample.set_index("time")["abundance"],
-            alpha=0.5, color="skyblue", zorder=-1
-        )
+
+for ax, species_i in zip(axes, species):
+    for tube_j in tubes:
+        for chain, draw in zip(chains, draws):
+            timeseries = prior_abundances.sel(chain=chain, draw=draw, tube=tube_j, species=species_i)
+            ax.plot(
+                timepoints.values, 
+                timeseries.values,
+                alpha=0.5, color="skyblue", zorder=-1
+            )
 f.savefig("img/monod_posteriors.png")
 ```
 

diff --git a/src/stan/monod.stan b/src/stan/monod.stan
@@ -68,9 +68,9 @@ model {
   a_mu_max ~ normal(prior_a_mu_max[1], prior_a_mu_max[2]);
   a_ks ~ normal(prior_a_ks[1], prior_a_ks[2]);
   a_gamma ~ normal(prior_a_gamma[1], prior_a_gamma[2]);
-  t_mu_max ~ normal(prior_t_mu_max[1], prior_t_mu_max[2]);
-  t_ks ~ normal(prior_t_ks[1], prior_t_ks[2]);
-  t_gamma ~ normal(prior_t_gamma[1], prior_t_gamma[2]);
+  t_mu_max ~ lognormal(prior_t_mu_max[1], prior_t_mu_max[2]);
+  t_ks ~ lognormal(prior_t_ks[1], prior_t_ks[2]);
+  t_gamma ~ lognormal(prior_t_gamma[1], prior_t_gamma[2]);
   for (s in 1 : 2) {
     sigma_y[s] ~ lognormal(prior_sigma_y[s, 1], prior_sigma_y[s, 2]);
     for (t in 1 : N_tube){