Skip to content

Commit 3403d1a

Browse files
authored
SSM hurricane Notebooks update (#804)
* StateSpace Module Hurricane Case Study WIP rebased main to resolve conflicts * attempt to fix plot rendering on ReviewNB * attempt to fix plot rendering on ReviewNB try jupyterlab renderer * attempt to fix plot rendering on ReviewNB try notebook renderer * 1. Made changes in accordance to reviewer comments 2. Improved 24-hour ahead forecasts of simplest model 3. Added a section and model with how to add exogenous variables to the model * 1. added forecast to exogenous variables section 2. cleaned up typos and figure sizes/zoom * 1. updated the exogenous covariates section 2. Added a non-linear B-Splines section * 1. Added text to spline section 2. Added closing remarks section 3. Updated references bib file * 1. fixed typos, latex errors, and improved readability of text 2. Restructured function docstrings to comply with numpy's guidestyle 3. Added ODE derivation of Newtonian equations of motion 4. Added legend to origins scatterplot 5. Fixed issue with passing in future covariates into exogenous model when forecasting 6. Added model evaluations and comparisons 7. Fixed issue with diffuse initialization of P0 * updated notebook removing exogenous forecasts hack * made some figures static to reduce file size * 1. updated docstring indentation 2. updated generate forecast utility 3. updated how mode=JAX is specified 4. renamed simple model to Newtonian model 5. added explanations to exogenous data duplicating and mirroring 6. added note on how the newtonian model differs from typical parameteric models * attempt to fix latex rendering * missed spot in cleaning up latex * updated equation number to attempt to fix latex rendering * removed equation numbering * fixed broken reference in closing remarks * updated headers. Only one H1 header in notebook * updated authoring dates and fixed missed h3 header * made change to utilize conditional posteriors for predictive state means and covariances. Updated plotly setting to automatically select renderer * Revert "made change to utilize conditional posteriors for predictive state means and covariances. Updated plotly setting to automatically select renderer" This reverts commit 07355d2. Revert using conditional posterior to using save_kalman_filter_outputs_in_idata * try mimetype plotly renderer * try sphinx_gallery plotly renderer * try modifying javascript to render correctly * try disabling plotly mathjax * try iframe renderer * changed plotly to static rendering * fixed inconsistency where transition matrix was sometimes referred to with A and sometimes with T to only being referred to using T * try remove-cell tag to change thumbnail * updated the way in which we sample filter_outputs using the sample_filter_outputs method * updated to use new arviz * added arviz-variat style to arviz plots
1 parent 1d8a59e commit 3403d1a

2 files changed

Lines changed: 418 additions & 223 deletions

File tree

examples/case_studies/ssm_hurricane_tracking.ipynb

Lines changed: 370 additions & 199 deletions
Large diffs are not rendered by default.

examples/case_studies/ssm_hurricane_tracking.myst.md

Lines changed: 48 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ jupytext:
55
format_name: myst
66
format_version: 0.13
77
kernelspec:
8-
display_name: pymc_examples_dev
8+
display_name: pymc-examples
99
language: python
1010
name: python3
1111
myst:
@@ -91,8 +91,8 @@ import warnings
9191
9292
warnings.filterwarnings("ignore", message="The RandomType SharedVariables", category=UserWarning)
9393
94-
import arviz as az
9594
import arviz.labels as azl
95+
import arviz.preview as az
9696
import numpy as np
9797
import pymc as pm
9898
import pytensor.tensor as pt
@@ -119,6 +119,8 @@ from pymc_extras.statespace.utils.constants import (
119119
120120
# make all plotly figures static
121121
pio.renderers.default = "svg"
122+
# set arviz style
123+
az.style.use("arviz-variat")
122124
```
123125

124126
## Helper Functions
@@ -313,9 +315,7 @@ def generate_period_forecasts(
313315
longitude_cppc = az.extract(forecasts["forecast_observed"].sel(observed_state="x"))
314316
latitude_cppc = az.extract(forecasts["forecast_observed"].sel(observed_state="y"))
315317
cppc_var = forecasts["forecast_observed"].var(("chain", "draw"))
316-
cppc_covs = xr.cov(
317-
latitude_cppc["forecast_observed"], longitude_cppc["forecast_observed"], dim="sample"
318-
)
318+
cppc_covs = xr.cov(latitude_cppc, longitude_cppc, dim="sample")
319319
covs_list = []
320320
for i in range(cppc_covs.shape[0]):
321321
covs_list.append(
@@ -964,23 +964,32 @@ with pm.Model(coords=n_ssm.coords) as newtonian:
964964
965965
n_ssm.build_statespace_graph(
966966
data=fiona_df.select("longitude", "latitude").to_numpy(),
967-
save_kalman_filter_outputs_in_idata=True,
968967
)
969968
newtonian_idata = pm.sample(
970969
nuts_sampler="nutpie", nuts_sampler_kwargs={"backend": "jax", "gradient_backend": "jax"}
971970
)
972971
```
973972

974973
```{code-cell} ipython3
975-
az.summary(newtonian_idata, var_names="acceleration_innovations", kind="stats")
974+
az.summary(newtonian_idata, var_names="acceleration_innovations", kind="stats", round_to=4)
975+
```
976+
977+
```{code-cell} ipython3
978+
n_ssm_filter_outputs = n_ssm.sample_filter_outputs(
979+
newtonian_idata, filter_output_names=["predicted_covariances", "predicted_observed_states"]
980+
)
976981
```
977982

978983
```{code-cell} ipython3
979-
predicted_covs = newtonian_idata.posterior["predicted_covariance"].mean(("chain", "draw"))
984+
predicted_covs = n_ssm_filter_outputs.posterior_predictive["predicted_covariances"].mean(
985+
("chain", "draw")
986+
)
980987
```
981988

982989
```{code-cell} ipython3
983-
post_mean = newtonian_idata.posterior["predicted_observed_state"].mean(("chain", "draw"))
990+
post_mean = n_ssm_filter_outputs.posterior_predictive["predicted_observed_states"].mean(
991+
("chain", "draw")
992+
)
984993
```
985994

986995
Not bad for a model with only one parameter. We can see that the forecast gets wonky in the middle where the trajectory of the Hurricane changes directions over short time periods. Again, it is important to keep in mind that what we are plotting are the one-step/period ahead forecast. In our case, our periods are six hours apart. Unfortunately, a 6-hour ahead hurricane forecast is not very practical. Let's see what we get when we generate a 4-period (24-hour) ahead forecast.
@@ -1343,10 +1352,7 @@ with pm.Model(coords=exog_ssm.coords) as exogenous:
13431352
13441353
acceleration_innovations = pm.Gamma("acceleration_innovations", 0.1, 5, shape=(1,))
13451354
1346-
exog_ssm.build_statespace_graph(
1347-
data=fiona_df.select("longitude", "latitude").to_numpy(),
1348-
save_kalman_filter_outputs_in_idata=True,
1349-
)
1355+
exog_ssm.build_statespace_graph(data=fiona_df.select("longitude", "latitude").to_numpy())
13501356
exogenous_idata = pm.sample(
13511357
nuts_sampler="nutpie", nuts_sampler_kwargs={"backend": "jax", "gradient_backend": "jax"}
13521358
)
@@ -1355,7 +1361,7 @@ with pm.Model(coords=exog_ssm.coords) as exogenous:
13551361
Typically, the surface wind speed and the central pressure of a hurricane carry little information on the path the hurricane will take. The path of a hurricane is, generally, influenced by surrounding atmospheric conditions like pressure gradients. Knowing this, it makes sense to see that many of our beta parameters are close to zero, indicating little to no influence on the hurricanes' path.
13561362

13571363
```{code-cell} ipython3
1358-
az.plot_trace(exogenous_idata, var_names="acceleration_innovations");
1364+
az.plot_trace_dist(exogenous_idata, var_names="acceleration_innovations");
13591365
```
13601366

13611367
```{code-cell} ipython3
@@ -1367,8 +1373,18 @@ az.plot_forest(
13671373
### Make in-sample forecasts with new exogenous model
13681374

13691375
```{code-cell} ipython3
1370-
predicted_covs = exogenous_idata.posterior["predicted_covariance"].mean(("chain", "draw"))
1371-
post_mean = exogenous_idata.posterior["predicted_observed_state"].mean(("chain", "draw"))
1376+
exog_ssm_filter_outputs = exog_ssm.sample_filter_outputs(
1377+
exogenous_idata, filter_output_names=["predicted_covariances", "predicted_observed_states"]
1378+
)
1379+
```
1380+
1381+
```{code-cell} ipython3
1382+
predicted_covs = exog_ssm_filter_outputs.posterior_predictive["predicted_covariances"].mean(
1383+
("chain", "draw")
1384+
)
1385+
post_mean = exog_ssm_filter_outputs.posterior_predictive["predicted_observed_states"].mean(
1386+
("chain", "draw")
1387+
)
13721388
```
13731389

13741390
Our one-period ahead forecasts seem to be slightly worse than our Newtonian model. You will notice that at the end of the forecast we see that our trajectory is erroneously more north rather than north-east. Since the exogenous variables we added to the model don't carry additional information with respect to the hurricane's trajectory, this results are expected.
@@ -1805,10 +1821,7 @@ with pm.Model(coords=spline_ssm.coords) as spline_model:
18051821
18061822
acceleration_innovations = pm.Gamma("acceleration_innovations", 0.1, 5, shape=(1,))
18071823
1808-
spline_ssm.build_statespace_graph(
1809-
data=fiona_df.select("longitude", "latitude").to_numpy(),
1810-
save_kalman_filter_outputs_in_idata=True,
1811-
)
1824+
spline_ssm.build_statespace_graph(data=fiona_df.select("longitude", "latitude").to_numpy())
18121825
spline_idata = pm.sample(
18131826
nuts_sampler="nutpie", nuts_sampler_kwargs={"backend": "jax", "gradient_backend": "jax"}
18141827
)
@@ -1817,11 +1830,11 @@ with pm.Model(coords=spline_ssm.coords) as spline_model:
18171830
Most of our spline parameters are around zero, with a handful of exceptions. Let's take a look at how these effect our forecasts.
18181831

18191832
```{code-cell} ipython3
1820-
az.plot_trace(spline_idata, var_names="acceleration_innovations");
1833+
az.plot_trace_dist(spline_idata, var_names="acceleration_innovations");
18211834
```
18221835

18231836
```{code-cell} ipython3
1824-
az.plot_trace(spline_idata, var_names=["beta_exog"], compact=True, figsize=(20, 8));
1837+
az.plot_trace_dist(spline_idata, var_names=["beta_exog"], compact=True);
18251838
```
18261839

18271840
### Make in-sample forecasts with new spline model
@@ -1831,8 +1844,18 @@ az.plot_trace(spline_idata, var_names=["beta_exog"], compact=True, figsize=(20,
18311844
Our one-period ahead forecasts, look better than the ones we generated from the Exogenous covariates model, but worse than the original model that purely follows Newtonian kinematics.
18321845

18331846
```{code-cell} ipython3
1834-
predicted_covs = spline_idata.posterior["predicted_covariance"].mean(("chain", "draw"))
1835-
post_mean = spline_idata.posterior["predicted_observed_state"].mean(("chain", "draw"))
1847+
spline_ssm_filter_outputs = spline_ssm.sample_filter_outputs(
1848+
spline_idata, filter_output_names=["predicted_covariances", "predicted_observed_states"]
1849+
)
1850+
```
1851+
1852+
```{code-cell} ipython3
1853+
predicted_covs = spline_ssm_filter_outputs.posterior_predictive["predicted_covariances"].mean(
1854+
("chain", "draw")
1855+
)
1856+
post_mean = spline_ssm_filter_outputs.posterior_predictive["predicted_observed_states"].mean(
1857+
("chain", "draw")
1858+
)
18361859
```
18371860

18381861
```{code-cell} ipython3
@@ -1941,6 +1964,7 @@ fig.show(width=1000, renderer="png", config={"displayModeBar": False})
19411964
```
19421965

19431966
## Authors
1967+
* Updated by Jonathan Dekermanjian in August, 2025 to use the `sample_filter_outputs` method
19441968
* Authored by Jonathan Dekermanjian in June, 2025
19451969

19461970
+++

0 commit comments

Comments
 (0)