<jupyter_start><jupyter_text>Data The data consits of COVID-19 positive test counts (`Daily lab-confirmed cases`) by Lower tier local authority (LTLA) as reported on a given date (`Report date`) for all preceding test dates (`Specimen date`).<jupyter_code>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')<jupyter_output><empty_output><jupyter_text>We model the last two weeks of data.<jupyter_code>ys_sequence = df['Daily lab-confirmed cases'].values[-14:]
specimen_dates = df['Specimen date'].unique()[-14:]<jupyter_output><empty_output><jupyter_text>The dashed line in the plot below indicates that counts for the most recent dates do not reflect the final count yet and are subject to change in future reports due to the lag in reporting (it takes about 4 days for all tests to be processed and positive counts reported).<jupyter_code>fig, ax = plt.subplots(1, 1, figsize=(14,4))

plt.plot(specimen_dates[:-3], ys_sequence[:-3], "xkcd:denim blue")
plt.plot(specimen_dates[-4:], ys_sequence[-4:], color="xkcd:denim blue", linestyle="--")

plt.title(f"{ltla}, {report_dt} report", fontsize=16)
plt.ylabel('Daily lab-confirmed cases', fontsize=16)
plt.xlabel('Specimen date', fontsize=16)
plt.setp(ax.get_xticklabels(), fontsize=14)
plt.setp(ax.get_yticklabels(), fontsize=14)
plt.setp(ax.xaxis.get_majorticklabels(), rotation=70)

plt.show()<jupyter_output><empty_output><jupyter_text>Reporting lag priors ($\theta_j$)<jupyter_code>from data_processing.data import get_data, moment_match_theta_priors<jupyter_output><empty_output><jupyter_text>A core assumption of the model is that the reporting lag is consistent and predictable. We load two weeks of prior reports to learn the most recent reporting behaviour from the data. Specifically we look at the proportion of the final count that is reported (`prop_reported`) at each `lag` (where `lag` is the difference in days between `Specimen date` and `Report date`). We denote the reporting proportion parameters $\theta_j$.<jupyter_code>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)
historical_df[['Specimen date', 'Report date', 'Daily lab-confirmed cases', 'prop_reported', 'lag']
             ].query("(`Specimen date` == '2020-12-03' | `Specimen date` == '2020-12-04') & lag <= 5"
                    ).sort_values(['Specimen date','Report date'])<jupyter_output><empty_output><jupyter_text>For each `lag` $j$ we propose that $\theta_j \sim Beta(\alpha_j, \beta_j)$ and we estimate  the $\alpha_j$ and $\beta_j$ shape parameters by moment matching to the data (`prop_reported`).<jupyter_code>alphas, betas = moment_match_theta_priors(historical_df, n_lags=len(ys_sequence))<jupyter_output><empty_output><jupyter_text>Below we plot the reporting data and the corresponding Beta prior for the first 5 lags.<jupyter_code>fig, ax = plt.subplots(1, 5, figsize=(14, 4))
x = np.linspace(0, 1, 1002)[1:-1]

for i in range(1, 6):
    curr_ax = ax[i-1]
    
    empirical = historical_df[historical_df['lag']==i][:-4]
    curr_ax.hist(empirical['prop_reported'], range=(0,1), bins=15, label='Data', alpha=0.25, density=True)

    dist = beta(alphas[-i], betas[-i])
    curr_ax.plot(x, dist.pdf(x), label='Prior', color="#0e1c2e")

    curr_ax.set_ylim(0,20)

    if i == 0:
        curr_ax.set_ylabel(f"{l}", fontsize=16)
    curr_ax.set_yticklabels([])

    curr_ax.set_title(f"Lag {i}", fontsize=16)

    if i==3:
        curr_ax.set_xlabel('$\theta$', fontsize=24)

    plt.setp(curr_ax.get_xticklabels(), fontsize=14)
    plt.setp(curr_ax.get_yticklabels(), fontsize=14)
    
ax[0].legend(fontsize=16)
plt.suptitle(f"", fontsize=16)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Model<jupyter_code>from models.inference import Nowcaster<jupyter_output><empty_output><jupyter_text>To instantiate the Nowcaster model class we have to specify (i) the reporting lag priors learned in the previous step and (ii) the scale ($\sigma$) of the temporal smoothing applied by the random walk prior. In practice we would run the model with a number of different scales and choose the model that maximises the model evidence.<jupyter_code>model = Nowcaster(
            random_walk_scale=4,
            prior_theta_shape_a = alphas,
            prior_theta_shape_b = betas
)<jupyter_output><empty_output><jupyter_text>Once the model is instantiated, we pass the reported counts and dates to the `model.fit()` method. Including the dates is required for the model to account for weekend effects.<jupyter_code>model.fit(
    observations=ys_sequence,
    dates=specimen_dates
)<jupyter_output><empty_output><jupyter_text>Model outputsThe Nowcaster class comes with methods to access and plot the model outputs once fit to data.<jupyter_code>model.plot_lambda_smoothing()
model.plot_x_smoothing()
model.plot_drift_smoothing()<jupyter_output><empty_output><jupyter_text>The samples plotted above can also be accessed directly as properties of the Nowcaster class:- `model.lambda_smoothing_samples`- `model.x_smoothing_samples`- `model.kappa_smoothing_samples`<jupyter_code>[x.mean() for x in model.x_smoothing_samples]<jupyter_output><empty_output><jupyter_text>The Nowcaster class also returns the model evidence.<jupyter_code>sum(model.log_evidence)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12151898734177215, 'fraction_numerical': 0.018806509945750453, 'mean_word_length': 5.862282878411911, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 53, 'https://': 0, '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': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '479683', 'license': 'MIT', 'n_tokens_mistral': 1959, 'n_tokens_neox': 1836, 'n_words': 485, 'path': 'notebooks/quick_start.ipynb', 'repo_name': 'anguswilliams91/jbc-turing-rss-nowcasting', 'size': '5530'}