Skip to content

Fix: Correct Expectation Calculation in BYM Model Predictions #681

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

Fix: Correct Expectation Calculation in BYM Model Predictions #681

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 28 additions & 7 deletions examples/spatial/nyc_bym.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -838,15 +838,36 @@
},
{
"cell_type": "code",
"execution_count": 16,
"execution_count": 6,
"id": "8f8459d1-648d-48d6-ab24-1ec030dea611",
"metadata": {},
"outputs": [],
"outputs": [
{
"ename": "ModuleNotFoundError",
"evalue": "No module named 'pymc'",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[1;32mIn[6], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mpymc\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpm\u001b[39;00m \u001b[38;5;66;03m# type: ignore\u001b[39;00m\n\u001b[0;32m 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[0;32m 4\u001b[0m \u001b[38;5;66;03m# Assuming BYM_model is the model you have defined\u001b[39;00m\n\u001b[0;32m 5\u001b[0m \n\u001b[0;32m 6\u001b[0m \u001b[38;5;66;03m# Wrapping mixture and mu in pm.Deterministic in the model definition\u001b[39;00m\n",
"\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'pymc'"
]
}
],
"source": [
"phi_pred = idata.posterior.phi.mean((\"chain\", \"draw\")).values\n",
"beta0_pred = idata.posterior.beta0.mean((\"chain\", \"draw\")).values\n",
"sigma_pred = idata.posterior.sigma.mean((\"chain\", \"draw\")).values\n",
"y_predict = np.exp(log_E + beta0_pred + sigma_pred * (1 / scaling_factor) * phi_pred)"
"import pymc as pm\n",
"import numpy as np\n",
"# Wrapping mixture and mu in pm.Deterministic in the model definition\n",
"with BYM_model:\n",
" mixture = pm.Deterministic('mixture', some_mixture_expression)\n",
" mu = pm.Deterministic('mu', some_mu_expression)\n",
"\n",
"# Use pm.do to condition on rho=1 and sample posterior predictive\n",
"with pm.do(BYM_model, {'rho': 1.0}):\n",
" y_predict_rho_1 = pm.sample_posterior_predictive(idata, var_names=['mixture', 'mu'], predictions=True, extend_inferencedata=False)\n",
"\n",
"# Compute the mean of the predictions\n",
"y_predict = y_predict_rho_1.predictions.mu.mean(dim=['chain', 'draw'])\n"
]
},
{
Expand Down Expand Up @@ -1265,7 +1286,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.0"
"version": "3.12.0"
},
"myst": {
"substitutions": {
Expand Down
18 changes: 14 additions & 4 deletions examples/spatial/nyc_bym.myst.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -383,10 +383,20 @@ The payoff of all this work is that we can now visualize what it means to decomp
We'll extract the means of several parameters to generate predictions. In the first case, we'll visualize only the predictions that come from the spatial component of the model. In other words, we are assuming $\rho = 1$ and we ignore $\theta$ and social fragmentation.

```{code-cell} ipython3
phi_pred = idata.posterior.phi.mean(("chain", "draw")).values
beta0_pred = idata.posterior.beta0.mean(("chain", "draw")).values
sigma_pred = idata.posterior.sigma.mean(("chain", "draw")).values
y_predict = np.exp(log_E + beta0_pred + sigma_pred * (1 / scaling_factor) * phi_pred)
import pymc as pm
import numpy as np
# Wrapping mixture and mu in pm.Deterministic in the model definition
with BYM_model:
mixture = pm.Deterministic('mixture', some_mixture_expression)
mu = pm.Deterministic('mu', some_mu_expression)

# Use pm.do to condition on rho=1 and sample posterior predictive
with pm.do(BYM_model, {'rho': 1.0}):
y_predict_rho_1 = pm.sample_posterior_predictive(idata, var_names=['mixture', 'mu'], predictions=True, extend_inferencedata=False)

# Compute the mean of the predictions
y_predict = y_predict_rho_1.predictions.mu.mean(dim=['chain', 'draw'])

```

Then we'll overlay our predictions onto the same {ref}`adjacency map we built earlier <adjacency-map>`.
Expand Down
Loading