<jupyter_start><jupyter_text>Stan Implementation of the modelIn the paper, the maximum marginal likelihood is used to select the scale for the random walk scale. This notebook uses a model where a hyperprior is placed on the random walk scale instead:$$\sigma_{rw} \sim \mathcal{N}^+(5)$$For inference, I use NUTS instead of SMC. 1. Preprocess the dataI'll just use the same example as in `quick_start.ipynb`.<jupyter_code>import pandas as pd

ltla = "Thurrock"
report_dt = "2020-12-14"

df = pd.read_csv(f"../data/cases-{report_dt}.csv"
                ).query(f"`Area name` == '{ltla}' & `Area type` == 'ltla'")
df = df.sort_values('Specimen date')

ys_sequence = df['Daily lab-confirmed cases'].values[-14:]
is_weekend = (pd.to_datetime(df["Specimen date"]).dt.weekday >= 5).astype(int).values[-14:]
import sys
sys.path.append("..")

from data_processing.data import get_data, moment_match_theta_priors
from datetime import timedelta

last_report_dt = pd.to_datetime(report_dt)
report_dates = pd.date_range(start=last_report_dt - timedelta(14), end=last_report_dt)

historical_df = get_data(ltla, report_dates=report_dates)
alphas, betas = moment_match_theta_priors(historical_df, n_lags=len(ys_sequence))<jupyter_output>../data_processing/data.py:95: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  row[clcc] = np.nan<jupyter_text>2. Fit the modelStan grumbles a bit that $\lambda$ can become negative sometimes because the AR(1) process can easily walk in that direction - a future version of this model might do this on the log scale instead.<jupyter_code>import pystan

model = pystan.StanModel("../stan_models/nowcast.stan")
data = {"n_days": len(ys_sequence), "y": ys_sequence, "a": alphas, "b": betas, "is_weekend": is_weekend}
fit = model.sampling(data)<jupyter_output>WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated
WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed
WARNING:pystan:2105 of 4000 iterations saturated the maximum tree depth of 10 (52.6 %)
WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation<jupyter_text>The Rhat warning here is because some elements of `z` have zero variance because they are deterministically set to 1, meaning Rhat is undefined.  3. Plots<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

def plot_posterior(samples, name=None):
    # plot the posterior mean and 95% credible interval
    mean = samples.mean(axis=0)
    lo = np.percentile(samples, 2.5, axis=0)
    hi = np.percentile(samples, 97.5, axis=0)
    
    fig, ax = plt.subplots(figsize=(7, 5))
    ax.fill_between(np.arange(samples.shape[1]), lo, hi, color="xkcd:sky blue")
    ax.plot(mean, c="xkcd:denim blue", lw=3, label="" or name)
    
    return ax 
ax = plot_posterior(fit["lambda"], name="$\lambda$")
ax.plot(ys_sequence, c="r", ls="--", lw=3, label="reported")
ax.legend(loc="upper left");
import seaborn as sns

fig, ax = plt.subplots(figsize=(7, 5))
sns.distplot(fit["sigma_rw"], ax=ax)
ax.set_xlabel("$\sigma_{rw}$", size=18);
ax = plot_posterior(fit["kappa"])
ax.set_ylabel("$\kappa$", size=18);
ax = plot_posterior(fit["x"], name="$x$")
ax.plot(ys_sequence, c="r", ls="--", lw=3, label="reported")
ax.legend(loc="upper left");<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13045912653975364, 'fraction_numerical': 0.019036954087346025, 'mean_word_length': 5.403225806451613, 'pattern_counts': {'":': 5, '<': 11, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '479197', 'license': 'MIT', 'n_tokens_mistral': 1294, 'n_tokens_neox': 1206, 'n_words': 350, 'path': 'notebooks/stan_version.ipynb', 'repo_name': 'anguswilliams91/jbc-turing-rss-nowcasting', 'size': '3572'}