<jupyter_start><jupyter_text>--- BACKUPCode snippets (partially) considered for the main story, most intended to be run after already acquiring the data,  calculating thte second dose delay, and defining utility functions.  weekly quantiles<jupyter_code>weekly_avg = df["newTot"].resample("W", label="left").mean()
# label="left" uses week start dates to label our averages, rather
# than week end dates (which is the default)
# only keep the date (not time) to tidy up axis labels
weekly_avg.index = weekly_avg.index.date

plt.figure(figsize=figsize)
weekly_avg.plot.barh()
plt.xlabel("Mean Daily Doses [millions]")
plt.ylabel("Week Starting")
quantiles = [0.1, 0.25, 0.75, 0.9]
weekly_avg.quantile(q=quantiles)<jupyter_output><empty_output><jupyter_text>The quantile values above show that 25% of the weeks to date had daily averages of 0.36 million doses administered per day or lower, and 25% had more than 0.43 million doses per day (or equivalently 75% had 0.43 million or lower).<jupyter_code>forecast_q = {}
daily_doses_q = {}

for q in quantiles:
    daily_doses_q[q] = weekly_avg.quantile(q=q)
    forecast_q[q] = forecast_const(
        df,
        avg_second_delay,
        daily_doses_q[q],
        uk_pop=priority_totals["All Adults"], 
    )

plot_cumulative_doses(
    df_forecast, forecast_date=last_data, figsize=(15, 8), 
)

q_fill = [
    # start quantile, end quantile, fill opacity (alpha)
    (0.1, 0.25, 0.1),
    (0.25, 0.75, 0.25),
    (0.75, 0.9, 0.1)
]

for q_start, q_end, alpha in q_fill:
    for column in ["cumFirst", "cumSecond"]:
        plt.fill_between(
            forecast_q[q_start].index,
            forecast_q[q_start][column],
            forecast_q[q_end][column],
            color=col_format[column]["color"],
            alpha=alpha,
        )<jupyter_output><empty_output><jupyter_text>Forecast with doses today = same as doses N days agoIf the previous week repeated, when will each group be vaccinated?<jupyter_code>def forecast_lookback(
    df,
    avg_second_delay,
    days_lookback=7,
    uk_pop=priority_totals["All Adults"], 
    end_date=datetime(2021, 12, 1),
):
    """
    Forecast vaccine data assuming the number of vaccines given
    today is the same as 'days_lookback' days ago.
    """
    def lookback(df, date, days=days_lookback):
        """
        Return total number of doses given a number of days before
        the input date.
        """
        if date - timedelta(days=days) < df.index.min():
            return df.iloc[0]["newTot"]
        else:
            return df.loc[date - timedelta(days=days), "newTot"]

    df_forecast = forecast_vaccines(
        df,
        avg_second_delay,
        doses_fn=lookback,
        uk_pop=uk_pop, 
        end_date=end_date,
    )

    return df_forecast

df_forecast = forecast_lookback(df, avg_second_delay)
plt.figure(figsize=figsize)

total = df_forecast["newTot"]
last_data = df.index.max()
total[total.index <= last_data].plot(color="k", label="Actual Total", linewidth=3)
total[total.index >= last_data].plot(color="k", linestyle="--", label="Forecast Total")

plt.legend(loc="upper left")
plt.ylim([0, 1.1 * total.max()])
plt.ylabel("Total Doses [millions]")<jupyter_output><empty_output><jupyter_text>Forecast - Random sample doses<jupyter_code>import numpy as np
from scipy.stats import lognorm
import matplotlib.pyplot as plt

data = df["newTot"].dropna()

s = lognorm.fit(data)
print(s)
plt.figure()
sns.histplot(data, stat="density", bins=20)
x = np.linspace(0, 1, 100)
y = lognorm.pdf(x, *s)
plt.plot(x, y)
plt.show()


n_forecasts = 100
rng = np.random.default_rng(seed=123)

def rnd_doses(df, date):
    return lognorm.rvs(*s, random_state=rng)


rnd_forecasts = [
    forecast_vaccines(
        df,
        avg_second_delay,
        doses_fn=rnd_doses,
        uk_pop=priority_totals["All Adults"], 
        end_date=datetime(2021, 12, 31),
        min_second_delay=28,
    )
    for _ in range(n_forecasts)
]

rnd_cumFirst = pd.concat(
    [rnd_forecasts[i]["cumFirst"].rename(f"forecast{i}") for i in range(n_forecasts)],
    axis=1, names=[0, 1]
)
rnd_cumSecond = pd.concat(
    [rnd_forecasts[i]["cumSecond"].rename(f"forecast{i}") for i in range(n_forecasts)],
    axis=1, names=[0, 1]
)

quantiles = [0.025, 0.25, 0.5, 0.75, 0.975]

q_1st = rnd_cumFirst.quantile(quantiles, axis=1)
q_2nd = rnd_cumSecond.quantile(quantiles, axis=1)

plot_cumulative_doses(
    pd.DataFrame({"cumFirst": q_1st.loc[0.5], "cumSecond": q_2nd.loc[0.5]}),
    forecast_date=last_data, figsize=(15, 8), 
)

q_fill = [
    # start quantile, end quantile, fill opacity (alpha)
    (0.025, 0.25, 0.1),
    (0.25, 0.75, 0.25),
    (0.75, 0.975, 0.1)
]

for q_start, q_end, alpha in q_fill:
    plt.fill_between(
        q_1st.loc[q_start].index,
        q_1st.loc[q_start],
        q_1st.loc[q_end],
        color=col_format["cumFirst"]["color"],
        alpha=alpha,
    )
    plt.fill_between(
        q_2nd.loc[q_start].index,
        q_2nd.loc[q_start],
        q_2nd.loc[q_end],
        color=col_format["cumSecond"]["color"],
        alpha=alpha,
    )
rnd_newTot = pd.concat(
    [rnd_forecasts[i]["newTot"].rename(f"forecast{i}") for i in range(n_forecasts)],
    axis=1, names=[0, 1]
)

q_tot = rnd_newTot.quantile(quantiles, axis=1)

plt.plot(q_tot.loc[0.975])
plt.plot(q_tot.loc[0.025])
dates = [
    d[d["cumSecond"] >= priority_totals["All Adults"] - 1e-9].index.min()
    for d in rnd_forecasts
]
print(np.argmin(dates), min(dates))
print(np.argmax(dates), max(dates))

plt.figure(figsize=figsize)
plt.plot(rnd_forecasts[np.argmin(dates)]["newTot"].rolling(window=7).mean())
plt.plot(rnd_forecasts[np.argmax(dates)]["newTot"].rolling(window=7).mean())
rnd_forecasts[np.argmax(dates)].loc[rnd_forecasts[np.argmax(dates)]["newTot"] > 0, "newTot"].mean()
#0.37 -> 0.44<jupyter_output><empty_output><jupyter_text>Stability of Forecasts<jupyter_code>def run_forecasts(
    df,
    run_start,
    days_lookback=7,
    uk_pop=priority_totals["All Adults"], 
    end_date=datetime(2022, 3, 1)
):
    """
    Calculate the completion date of the UK vaccination programme
    with different forecast start dates. A forecast is run as if
    it was every day starting from the date run_start, up to the
    date of the latest available data.
    """
    forecast_start_dates = pd.date_range(    
        start=run_start, end=df.index.max()
    )
    # fill missing second dose delay values with nearest available value
    df["delaySecond"] = df["delaySecond"].bfill()

    completion_dates = pd.Series(index=forecast_start_dates, dtype=float)
    for forecast_start in forecast_start_dates:
        second_delay = df.loc[forecast_start - timedelta(days=1), "delaySecond"]
        df_forecast = forecast_lookback(
            df[df.index <= forecast_start],
            second_delay,
            days_lookback=days_lookback,
            uk_pop=uk_pop, 
            end_date=end_date
        )
        if df_forecast["cumSecond"].max() >= (uk_pop - 1e-7):
            finish_date = df_forecast[df_forecast["cumSecond"] >= (uk_pop - 1e-7)].iloc[0].name
            completion_dates.loc[forecast_start] = finish_date
        else:
            print(forecast_start.date(), ": Population not vaccinated by ", end_date.date())
            completion_dates.loc[forecast_start] = end_date
    
    return completion_dates
completion_dates = run_forecasts(df, datetime(2021, 2, 1))
plt.figure(figsize=figsize)
completion_dates.plot()
plt.ylabel("Date Adult Population Fully Vaccinated")
plt.xlabel("Forecast Date")<jupyter_output><empty_output><jupyter_text>Effect of Easter Bank HolidayStabilised mid-September?There are several forecasts between late February and early March where the vaccine programme is not completed until January 2022, 2-3 months later than most forecasts at that time that have end dates in October. This looks a bit strange (why such big jumps?), but if we plot one of the forecasts with a late completion date we can see what's happening:<jupyter_code>df_feb24 = forecast_lookback(
    df[df.index <= datetime(2021, 2, 24)],
    df.loc[datetime(2021, 2, 24), "delaySecond"],
    end_date=datetime(2022, 2, 1),
)
plot_cumulative_doses(
    df_feb24,
    forecast_date=datetime(2021, 2, 24),
    figsize=(15, 8), 
    title=f"UK Vaccination Forecast (using data up to {datetime(2021, 2, 24).date()})"
)<jupyter_output><empty_output><jupyter_text>In the plot above, we see that in this forecast almost all adults are fully vaccinated by 20th August. However, before the remaining first doses can be given, a new batch of people requiring their second dose appear and must be given priority (without spare capacity for new first doses). 73 days later, on 1st November, the last few adults can start to be given their first dose. You'll notice that this 73 day delay in giving the last first doses is closely linked to the 72 day gap we assumed earlier between giving first doses and second doses. We can think of vaccines being given in alternating 72 day cycles of 1st doses and 2nd doses (at least in our simple forecast), and in the example above we have had to go through one additional cycle to complete the vaccination programme (compared to most forecasts which have completion datea in autumn 2021).BUT- assuming whole pop gets vaccinated etc.- delay 2nd doses few days more in this case?  prophet<jupyter_code>from prophet import Prophet

df_prophet = pd.DataFrame(
    {"ds": df.index, "y": df["newTot"], "cap": df["cumTot"]}
)

uk_pop = priority_totals["All Adults"]
#df_prophet["cap"] = (2 * uk_pop) - df_prophet["cap"]
df_prophet["cap"] = df["newTot"].max()
df_prophet["floor"] = 0

df_prophet

m = Prophet(interval_width=0.95, changepoint_prior_scale=0.1, growth="logistic")
m.add_country_holidays(country_name='UK')
m.fit(df_prophet)
#m.fit(df_prophet[df_prophet.index > datetime(2021, 2, 1)])
#df_prophet_w = df_prophet.resample("W-MON").mean()
#df_prophet_w["ds"] = df_prophet_w.index
#m.fit(df_prophet_w)

future = m.make_future_dataframe(periods=150)
future["cap"] = df["newTot"].max()
future["floor"] = 0

forecast = m.predict(future)
forecast
m.plot(forecast);
from prophet.plot import add_changepoints_to_plot
fig = m.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), m, forecast)
m.plot_components(forecast);
samples = m.predictive_samples(future)
samples["trend"]
df["newTot"].plot(marker="o", linestyle="None")
m.predict(df_prophet).set_index("ds")["yhat"].plot()
forecast[["ds", "yhat_lower", "yhat", "yhat_upper"]]
from prophet.diagnostics import cross_validation
df_cv = cross_validation(m, horizon = '14 days')
df_cv
from prophet.diagnostics import performance_metrics
df_p = performance_metrics(df_cv)
df_p
from prophet.plot import plot_cross_validation_metric
fig = plot_cross_validation_metric(df_cv, metric='mape')

df.resample("W-MON").mean()["newTot"].plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12586318018598655, 'fraction_numerical': 0.025872387441303747, 'mean_word_length': 4.182251908396947, 'pattern_counts': {'":': 6, '<': 33, '<?xml version=': 0, '>': 36, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '657956', 'license': 'MIT', 'n_tokens_mistral': 4066, 'n_tokens_neox': 3773, 'n_words': 926, 'path': 'backup.ipynb', 'repo_name': 'jack89roberts/covid-vax', 'size': '10861'}