Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
587 changes: 492 additions & 95 deletions examples/case_studies/GEV.ipynb

Large diffs are not rendered by default.

45 changes: 22 additions & 23 deletions examples/case_studies/GEV.myst.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ jupytext:
format_name: myst
format_version: 0.13
kernelspec:
display_name: default
display_name: eabm
language: python
name: python3
---
Expand Down Expand Up @@ -38,14 +38,13 @@ Note that this parametrization of the shape parameter $\xi$ is opposite in sign
We will use the example of the Port Pirie annual maximum sea-level data used in {cite:t}`coles2001gev`, and compare with the frequentist results presented there.

```{code-cell} ipython3
import arviz as az
import arviz.preview as az
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm
import pymc_extras.distributions as pmx
import pytensor.tensor as pt

from arviz.plots import plot_utils as azpu
az.style.use("arviz-variat")
```

## Data
Expand Down Expand Up @@ -112,18 +111,13 @@ Let's get a feel for how well our selected priors cover the range of the data:

```{code-cell} ipython3
idata = pm.sample_prior_predictive(samples=1000, model=model)
az.plot_ppc(idata, group="prior", figsize=(12, 6))
ax = plt.gca()
ax.set_xlim([2, 6])
ax.set_ylim([0, 2]);
az.plot_ppc_dist(idata, group="prior_predictive", kind="ecdf")
```

And we can look at the sampled values of the parameters, using the `plot_posterior` function, but passing in the `idata` object and specifying the `group` to be `"prior"`:

```{code-cell} ipython3
az.plot_posterior(
idata, group="prior", var_names=["μ", "σ", "ξ"], hdi_prob="hide", point_estimate=None
);
az.plot_dist(idata, group="prior", var_names=["μ", "σ", "ξ"]);
```

## Inference
Expand All @@ -144,7 +138,7 @@ idata.extend(trace)
```

```{code-cell} ipython3
az.plot_trace(idata, var_names=["μ", "σ", "ξ"], figsize=(12, 12));
az.plot_trace_dist(idata, var_names=["μ", "σ", "ξ"]);
```

### Divergences
Expand All @@ -159,27 +153,32 @@ The trace exhibits divergences (usually). The HMC/NUTS sampler can have problems
The 95% credible interval range of the parameter estimates is:

```{code-cell} ipython3
az.hdi(idata, hdi_prob=0.95)
az.hdi(idata, prob=0.95)
```

And examine the prediction distribution, considering parameter variability (and without needing to assume normality):

```{code-cell} ipython3
az.plot_posterior(idata, hdi_prob=0.95, var_names=["z_p"], round_to=4);
az.plot_dist(
idata,
ci_prob=0.95,
var_names=["z_p"],
stats={"point_estimate": {"round_to": 2}},
);
```

And let's compare the prior and posterior predictions of $z_p$ to see how the data has influenced things:

```{code-cell} ipython3
az.plot_dist_comparison(idata, var_names=["z_p"]);
az.plot_prior_posterior(idata, var_names=["z_p"]);
```

## Comparison
To compare with the results given in {cite:t}`coles2001gev`, we approximate the maximum likelihood estimates (MLE) using the mode of the posterior distributions (the *maximum a posteriori* or MAP estimate). These are close when the prior is reasonably flat around the posterior estimate.

The MLE results given in {cite:t}`coles2001gev` are:

$$\left(\hat{\mu}, \hat{\sigma}, \hat{\xi} \right) = \left( 3.87, 0.198, -0.050 \right) $$
$$\left(\hat{\mu}, \hat{\sigma}, \hat{\xi} \right) = \left( 3.87, 0.198, -0.050 \right)$$


And the variance-covariance matrix of the estimates is:
Expand All @@ -189,13 +188,8 @@ $$ V = \left[ \begin{array} 0.000780 & 0.000197 & -0.00107 \\
-0.00107 & -0.000778 & 0.00965
\end{array} \right] $$


Note that extracting the MLE estimates from our inference involves accessing some of the Arviz back end functions to bash the xarray into something examinable:

```{code-cell} ipython3
_, vals = az.sel_utils.xarray_to_ndarray(idata["posterior"], var_names=["μ", "σ", "ξ"])
mle = [azpu.calculate_point_estimate("mode", val) for val in vals]
mle
az.mode(idata, var_names=["μ", "σ", "ξ"])
```

```{code-cell} ipython3
Expand All @@ -207,12 +201,17 @@ The results are a good match, but the benefit of doing this in a Bayesian settin
Finally, we examine the pairs plots and see where any difficulties in inference lie using the divergences

```{code-cell} ipython3
az.plot_pair(idata, var_names=["μ", "σ", "ξ"], kind="kde", marginals=True, divergences=True);
az.plot_pair(
idata,
var_names=["μ", "σ", "ξ"],
visuals={"divergence": True},
);
```

## Authors

* Authored by [Colin Caprani](https://github.com/ccaprani), October 2021
* Updated by Osvaldo Martin, January 2026

+++

Expand Down
313 changes: 107 additions & 206 deletions examples/case_studies/factor_analysis.ipynb

Large diffs are not rendered by default.

67 changes: 31 additions & 36 deletions examples/case_studies/factor_analysis.myst.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ jupytext:
format_name: myst
format_version: 0.13
kernelspec:
display_name: Python 3 (ipykernel)
display_name: eabm
language: python
name: python3
myst:
Expand All @@ -33,7 +33,7 @@ Factor analysis is a widely used probabilistic model for identifying low-rank st
:::

```{code-cell} ipython3
import arviz as az
import arviz.preview as az
import numpy as np
import pymc as pm
import pytensor.tensor as pt
Expand All @@ -42,7 +42,6 @@ import seaborn as sns
import xarray as xr

from matplotlib import pyplot as plt
from matplotlib.lines import Line2D
from numpy.random import default_rng
from xarray_einstats import linalg
from xarray_einstats.stats import XrContinuousRV
Expand All @@ -52,7 +51,7 @@ print(f"Running on PyMC v{pm.__version__}")

```{code-cell} ipython3
%config InlineBackend.figure_format = 'retina'
az.style.use("arviz-darkgrid")
az.style.use("arviz-variat")

np.set_printoptions(precision=3, suppress=True)
RANDOM_SEED = 31415
Expand Down Expand Up @@ -128,11 +127,13 @@ with pm.Model(coords=coords) as PPCA:
At this point, there are already several warnings regarding failed convergence checks. We can see further problems in the trace plot below. This plot shows the path taken by each sampler chain for a single entry in the matrix $W$ as well as the average evaluated over samples for each chain.

```{code-cell} ipython3
for i in trace.posterior.chain.values:
samples = trace.posterior["W"].sel(chain=i, observed_columns=3, latent_columns=1)
plt.plot(samples, label=f"Chain {i + 1}")
plt.axhline(samples.mean(), color=f"C{i}")
plt.legend(ncol=4, loc="upper center", fontsize=12, frameon=True), plt.xlabel("Sample");
az.plot_trace(
trace,
var_names="W",
coords={"observed_columns": 3, "latent_columns": 1},
sample_dims=["draw"],
figure_kwargs={"sharey": True},
);
```

Each chain appears to have a different sample mean and we can also see that there is a great deal of autocorrelation across chains, manifest as long-range trends over sampling iterations.
Expand Down Expand Up @@ -194,13 +195,7 @@ with pm.Model(coords=coords) as PPCA_identified:
F = pm.Normal("F", dims=("latent_columns", "rows"))
sigma = pm.HalfNormal("sigma", 1.0)
X = pm.Normal("X", mu=W @ F, sigma=sigma, observed=Y, dims=("observed_columns", "rows"))
trace = pm.sample(tune=2000, random_seed=rng) # target_accept=0.9

for i in range(4):
samples = trace.posterior["W"].sel(chain=i, observed_columns=3, latent_columns=1)
plt.plot(samples, label=f"Chain {i + 1}")

plt.legend(ncol=4, loc="lower center", fontsize=8), plt.xlabel("Sample");
trace = pm.sample(tune=2000, random_seed=rng, target_accept=0.9)
```

$W$ (and $F$!) now have entries with identical posterior distributions as compared between sampler chains, although it's apparent that some autocorrelation remains.
Expand Down Expand Up @@ -251,29 +246,28 @@ When we compare the posteriors calculated using MCMC and VI, we find that (for a
```{code-cell} ipython3
col_selection = dict(observed_columns=3, latent_columns=1)

ax = az.plot_kde(
trace.posterior["W"].sel(**col_selection).values,
label=f"MCMC posterior for the explicit model",
plot_kwargs={"color": f"C{1}"},
)

az.plot_kde(
trace_amortized.posterior["W"].sel(**col_selection).values,
label="MCMC posterior for amortized inference",
plot_kwargs={"color": f"C{2}", "linestyle": "--"},
dt = az.from_dict(
{
"posterior": {
"MCMC_explicit": trace.posterior["W"].sel(**col_selection),
"MCMC_amortized": trace_amortized.posterior["W"].sel(**col_selection),
"FR-ADVI_amortized": trace_vi.posterior["W"].sel(**col_selection),
}
}
)


az.plot_kde(
trace_vi.posterior["W"].sel(**col_selection).squeeze().values,
label="FR-ADVI posterior for amortized inference",
plot_kwargs={"alpha": 0},
fill_kwargs={"alpha": 0.5, "color": f"C{0}"},
pc = az.plot_dist(
dt,
cols=None,
aes={"color": ["__variable__"]},
visuals={
"title": False,
"point_estimate_text": False,
"point_estimate": False,
"credible_interval": False,
},
)


ax.set_title(rf"PDFs of $W$ estimate at {col_selection}")
ax.legend(loc="upper left", fontsize=10);
pc.add_legend("__variable__")
```

### Post-hoc identification of F
Expand Down Expand Up @@ -389,6 +383,7 @@ We find that our model does a decent job of capturing the variation in the origi
* Updated by [Christopher Krapu](https://github.com/ckrapu) on April 4, 2021
* Updated by Oriol Abril-Pla to use PyMC v4 and xarray-einstats on March, 2022
* Updated by Erik Werner on Dec, 2023 ([pymc-examples#612](https://github.com/pymc-devs/pymc-examples/pull/612))
* Updated by Osvaldo Martin on January, 2026

+++

Expand Down
Loading