Skip to content

Commit 2ead8fe

Browse files
committed
1. Added text to spline section
2. Added closing remarks section 3. Updated references bib file
1 parent af68743 commit 2ead8fe

File tree

3 files changed

+221
-109
lines changed

3 files changed

+221
-109
lines changed

examples/case_studies/ssm_hurricane_tracking.ipynb

+172-106
Large diffs are not rendered by default.

examples/case_studies/ssm_hurricane_tracking.myst.md

+43-3
Original file line numberDiff line numberDiff line change
@@ -1129,7 +1129,7 @@ f_mean, cppc_vcov = generate_period_forecasts(
11291129
)
11301130
```
11311131

1132-
Compared to the previous 24-hour forecast, it looks like the uncertainty of the model has decreased, generally. We can also see that the mid-section and even the mid-to-end section have slightly improved. However, it still seems that our finally trajectory is more north than north-east, but is less north than our one-period ahead forecast.
1132+
Compared to the previous 24-hour forecast, it looks like the uncertainty of the model has decreased, generally. We can also see that the mid-section and even the mid-to-end section have slightly improved. However, it still seems that our finally trajectory is more north than north-east, but is less north than our last one-period ahead forecast.
11331133

11341134
```{code-cell} ipython3
11351135
fig = plot_hurricane_path(
@@ -1139,7 +1139,11 @@ fig.show(config={"displayModeBar": False})
11391139
```
11401140

11411141
# Add B-Splines
1142-
Still need to add text to this section
1142+
In the previous section, we tried adding an interaction term between the maximum sustained surface wind speed and the minumum central pressure. However, our estimated parameters were not too far off from zero. In this section we are going to attempt to model the non-linear complexities of the path, particularyly in the mid-section, using cubic B-splines.
1143+
1144+
To do this we first need to define what variables we are going to model as a smooth function. In our case, we are going to model the longitude values as a smooth function of the latitude values and vice versa.
1145+
1146+
To keep things simple, we are going to define a constant number of knots that are equal in both variables (same #knots for both longitude and latitude) and we are going to place the knots using a quantile function over the variable space. You can see the knots plotted out below for each variable.
11431147

11441148
```{code-cell} ipython3
11451149
num_knots = 15
@@ -1154,7 +1158,7 @@ fig.add_traces(
11541158
go.Scatter(
11551159
x=fiona_df["longitude"],
11561160
y=fiona_df["latitude"],
1157-
name="actuals",
1161+
name="Actuals",
11581162
mode="lines+markers",
11591163
line=dict(color="black"),
11601164
hovertemplate=[
@@ -1211,6 +1215,8 @@ fig.update_layout(
12111215
fig.show(config={"displayModeBar": False})
12121216
```
12131217

1218+
Next we need to create the basis functions over the defined variable space knot locations for each variable
1219+
12141220
```{code-cell} ipython3
12151221
B_longitude = dmatrix(
12161222
"bs(longitude, knots=longitude_knots, degree=3, include_intercept=True) - 1",
@@ -1227,6 +1233,25 @@ B_latitude = dmatrix(
12271233
exog_data = np.column_stack((np.asarray(B_longitude, order="F"), np.asarray(B_latitude, order="F")))
12281234
```
12291235

1236+
Our new models' structure is going to be similar to our last model that had exogenous variables. However, in this case our data are going to be the basis functions we created earlier. These will be inserted into our design matrix ($Z$) and the beta parameters corresponding to each spline will be added to our state vector ($x_{t}$). Again, these parameters will be constant (non-time varying). We will also have to re-adjust our transition matrix ($T$) and selection matrix ($R$) similar to how we did previously. We will now have:
1237+
1238+
$$\hat{y}_{longitude_{t+1}} = longitude_{t} + longitude\_velocity_{t}\Delta t + \frac{longitude\_acceleration_{t}\Delta t^{2}}{2} + \beta_{spline\_params_{longitude}} longitude\_spline\_basis\_functions$$
1239+
1240+
$$\hat{y}_{latitude_{t+1}} = latitude_{t} + latitude\_velocity_{t}\Delta t + \frac{latitude\_acceleration_{t}\Delta t^{2}}{2} + \beta_{spline\_params_{latitude}} latitude\_spline\_basis\_functions$$
1241+
1242+
our design matrix will be: $$ Z' = \begin{bmatrix} Z & X_{basis\_functions} \end{bmatrix} $$ Where $Z$ is our previously defined design matrix and $X_{exogenous}$ are the basis functions we defined earlier.
1243+
1244+
Our transition matrix will be $$T' = \begin{bmatrix} T & F \\ 0 & I_{n\_spline\_params} \end{bmatrix} $$
1245+
1246+
Where $T$ is the transition matrix defined for the Newtonian kinematics (the top-left 6x6 block in our previous model)
1247+
and $$ F = \begin{bmatrix} 1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0 \\ 0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1&0&1 \\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0 \\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0 \\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0 \\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0 \end{bmatrix} $$
1248+
and the 0 in the matrix above is a matrix of 0s of shape (number of spline parameters by number of spline parameters)
1249+
1250+
Finally, we have $$R' = \begin{bmatrix} R \\ 0 \end{bmatrix} $$
1251+
1252+
Where R is the selectiom matrix over our endogenous states (identity matrix of shape (number of states))
1253+
and again the 0 in the matrix is a matrix of 0s with shape (number of spline parameters by number of states)
1254+
12301255
```{code-cell} ipython3
12311256
class SplineSSM(PyMCStateSpace):
12321257
def __init__(self, k_exog: int = None):
@@ -1522,6 +1547,8 @@ with pm.Model(coords=spline_ssm.coords) as spline_model:
15221547
spline_idata = pm.sample(nuts_sampler="numpyro", target_accept=0.95)
15231548
```
15241549

1550+
Most of our spline parameters are around zero, with a handful of exceptions. Let's take a look at how these effect our forecasts.
1551+
15251552
```{code-cell} ipython3
15261553
az.plot_trace(spline_idata, var_names="acceleration_noise");
15271554
```
@@ -1532,6 +1559,10 @@ az.plot_trace(spline_idata, var_names=["beta_exog"], compact=True, figsize=(20,
15321559

15331560
# Make in-sample forecasts with new spline model
15341561

1562+
+++
1563+
1564+
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.
1565+
15351566
```{code-cell} ipython3
15361567
predicted_covs = spline_idata.posterior["predicted_covariance"].mean(("chain", "draw"))
15371568
post_mean = spline_idata.posterior["predicted_observed_state"].mean(("chain", "draw"))
@@ -1547,6 +1578,8 @@ fig = plot_hurricane_path(
15471578
fig.show(config={"displayModeBar": False})
15481579
```
15491580

1581+
Our 24-hour (4-period) forecasts, look pretty good. So far, this follows the true trajectory during the mid-section the best.
1582+
15501583
```{code-cell} ipython3
15511584
# our last i will equal 60 so we need to pad our exogenous data to len 64
15521585
exog_data_padded = np.concatenate(
@@ -1572,6 +1605,11 @@ fig = plot_hurricane_path(
15721605
fig.show(config={"displayModeBar": False})
15731606
```
15741607

1608+
# Closing Remarks
1609+
In this case study we looked at how we can track a hurricane in two-dimensional space using a state space representation of Newtonian kinematics. We proceeded to expand on the pure Newtonian model and added exogenous variables that may hold information pertintent to the Hurricane's track. We then expanded our model by modeling our variables as smooth functions using cubic B-splines. Throughout, the case study we saw that the pure Newtonian kinematics model's one-period ahead forecast performed the best. However, adding exogenous variables allowed our 24-hour (four-period) ahead forecasts to produce better forecasts compared to the pure Newtonian model. A possible explanation is that for a short period forecast the Newtonian kinematics contain all the information needed to produce the best next forecast. However, for longer periods the Newtonian kinematics are too deterministic whereaas the exogenous variables may act to regularize the output from the Newtonian kinematics. In turn producing better longer period forecasts.
1610+
1611+
+++
1612+
15751613
# Authors
15761614

15771615
+++
@@ -1580,6 +1618,8 @@ fig.show(config={"displayModeBar": False})
15801618

15811619
:::{bibliography}
15821620
:filter: docname in docnames
1621+
1622+
becker2023kalman
15831623
:::
15841624

15851625
+++

examples/references.bib

+6
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,12 @@ @article{bauer2005probing
5353
year = {2005},
5454
publisher = {Taylor \& Francis}
5555
}
56+
@book{becker2023kalman,
57+
title = {Kalman Filter from the Ground Up},
58+
author = {Becker, Alex},
59+
year = {2023},
60+
publisher = {Alex Becker}
61+
}
5662
@book{berry1996statistics,
5763
title = {Statistics: a Bayesian perspective},
5864
author = {Berry, Donald A},

0 commit comments

Comments
 (0)