<jupyter_start><jupyter_text>Generalized Gaussian Process Using PyMC3's Latent Gaussian Process implementation- API: https://docs.pymc.io/api/gp/implementations.html- Example: https://docs.pymc.io/notebooks/GP-Latent.htmlExample data: `boston_medical_center_2020-04-29_to_2020-06-22.csv`- Column: `hospitalized_total_covid_patients_suspected_and_confirmed_including_icu`  GoalGiven a series of past daily counts (admissions, census, etc.) $$y_1, y_2, ..., y_T$$Assuming $T$ is today, want to predict counts $$y_{T+1}, y_{T+2}, ..., y_{T+F}$$ for $F$ days ahead.  ModelSuppose that $y$ is Poisson distributed over the exponential of a latent Gaussian Process, i.e.,$$y_t \sim \text{Poisson}( \exp(f_t) )$$where $f$ is modeled by a Gaussian Process$$ f_t \sim N(m(t), k(t,t'))$$with constant mean$$m(t) = c$$and squared exponential covariance$$k(t,t') = a^2 \exp\left(-\frac{(t-t')^2}{2l^2}\right)$$  Parameters and their PriorsGP mean: $$c \sim \text{TruncatedNormal}(4, 2, \text{low}=0)$$SqExp cov amplitude: $$a \sim \text{HalfNormal}(\sigma=2)$$SqExp cov time-scale: $$l \sim \text{TruncatedNormal}(10, 2, \text{low}=0)$$  TrainingLet the subscript "past" represent indices $1$ through $T$, and the subscript "future" represent indices $T+1$ through $T+F$.1. Specify a Latent GP with the mean and covariance functions defined above.1. Define a `prior` distribution over all $f$ (i.e., condition on both $t_\text{past}$ and $t_\text{future}$).1. Define the observed variable $$y_\text{past} \sim \text{Poisson}(\mu = \exp(f_\text{past}))$$ Set `mu` to be the exponential of the first $T$ $f$ values, i.e., only $f_\text{past}$.1. Use PyMC3's MCMC to draw samples from the posterior $$c^s, a^s, l^s, f^s \sim p(c, a, l, f_\text{past}, f_\text{future} | y_\text{past})$$<jupyter_code>import pymc3 as pm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import theano
import theano.tensor as tt
theano.config.gcc.cxxflags = "-Wno-c++11-narrowing"
import logging
logger = logging.getLogger('arviz')
logger.setLevel(logging.ERROR)
df = pd.read_csv('../mass_dot_gov_datasets/boston_medical_center_2020-04-29_to_2020-06-22.csv')
y = df['hospitalized_total_covid_patients_suspected_and_confirmed_including_icu'].astype(float)
T = len(y)
F = 7
t = np.arange(T+F)[:,None]
with pm.Model() as model:
    c = pm.TruncatedNormal('mean', mu=4, sigma=2, lower=0)
    mean_func = pm.gp.mean.Constant(c=c)
    
    a = pm.HalfNormal('amplitude', sigma=2)
    l = pm.TruncatedNormal('time-scale', mu=10, sigma=2, lower=0)
    cov_func = a**2 * pm.gp.cov.ExpQuad(input_dim=1, ls=l)
    
    gp = pm.gp.Latent(mean_func=mean_func, cov_func=cov_func)
    
    f = gp.prior('f', X=t)

    y_past = pm.Poisson('y_past', mu=tt.exp(f[:T]), observed=y)
    y_logp = pm.Deterministic('y_logp', y_past.logpt)
with model:
    trace = pm.sample(5000, tune=1000, target_accept=.99, random_seed=42, chains=1)
summary = pm.summary(trace)['mean'].to_dict()
for key in ['mean', 'amplitude', 'time-scale']:
    print(key, summary[key])

print('\nTraining score:')
print(np.log(np.mean(np.exp(trace.get_values('y_logp', chains=0)))) / T)
pm.traceplot(trace);<jupyter_output><empty_output><jupyter_text>Forecasting Procedure1. Define a predictive distribution for future $y$ values $$y_\text{future} \sim \text{Poisson}(\mu = \exp(f_\text{future}))$$1. Use PyMC3's `sample_posterior_predictive` and the posterior samples collected during training to produce forecasts.<jupyter_code>with model:
    y_future = pm.Poisson('y_future', mu=tt.exp(f[-F:]), shape=F)
    forecasts = pm.sample_posterior_predictive(trace, vars=[y_future], random_seed=42)

samples = forecasts['y_future']
low = np.zeros(F)
high = np.zeros(F)
mean = np.zeros(F)
median = np.zeros(F)

for i in range(F):
    low[i] = np.percentile(samples[:,i], 2.5)
    high[i] = np.percentile(samples[:,i], 97.5)
    median[i] = np.percentile(samples[:,i], 50)
    mean[i] = np.mean(samples[:,i])

plt.figure(figsize=(8,6))
x_future = np.arange(1,F+1)
plt.errorbar(x_future, median,
             yerr=[median-low, high-median],
             capsize=2, fmt='.', linewidth=1,
             label='2.5, 50, 97.5 percentiles');
plt.plot(x_future, mean, 'x', label='mean');
x_past = np.arange(-4,1)
plt.plot(x_past, y[-5:], 's', label='observed')
plt.legend();
plt.title('Forecasts');
plt.xlabel('Days ahead');<jupyter_output><empty_output><jupyter_text>Heldout Scoring Procedure1. Partition the data to treat the first 80% as "past" and the last 20% as "future." We'll use $y_\text{past}$ as the training set, and $y_\text{future}$ as the validation set.1. Train the model using $y_\text{past}$.1. Define the predictive distribution $$y_\text{future} \sim \text{Poisson}(\mu = \exp(f_\text{future}))$$ Set `observed` to be the observed $y_\text{future}$ values.1. Define a `Deterministic` distribution that computes the logp of the observed variable $y_\text{future}$.1. Use `sample_posterior_predictive` to compute the log probability of $y_\text{future}$ conditioned on each posterior sample {$c^s, a^s, l^s, f^s$}.1. Use Monte Carlo integration to estimate the log probability of the heldout set: $$\log p(y_\text{future} | y_\text{past}) = \log \frac{1}{S} \sum_{s=1}^S p(y_\text{future} | c^s, a^s, l^s, f^s, y_\text{past})$$<jupyter_code>df = pd.read_csv("../mass_dot_gov_datasets/boston_medical_center_2020-04-29_to_2020-06-22.csv")
y = df['hospitalized_total_covid_patients_suspected_and_confirmed_including_icu'].astype(float)

T = int(.8 * len(y))
y_tr = y[:T]
y_va = y[T:]
F = len(y_va)

t = np.arange(T+F)[:,None]
with pm.Model() as model:
    c = pm.TruncatedNormal('mean', mu=4, sigma=2, lower=0)
    mean_func = pm.gp.mean.Constant(c=c)
    
    a = pm.HalfNormal('amplitude', sigma=2)
    l = pm.TruncatedNormal('time-scale', mu=20, sigma=5, lower=0)
    cov_func = a**2 * pm.gp.cov.ExpQuad(input_dim=1, ls=l)
    
    gp = pm.gp.Latent(mean_func=mean_func, cov_func=cov_func)
    
    f = gp.prior('f', X=t)

    y_past = pm.Poisson('y_past', mu=tt.exp(f[:T]), observed=y_tr)
    y_logp = pm.Deterministic('y_logp', y_past.logpt)
with model:
    trace = pm.sample(5000, tune=1000, chains=1, target_accept=.99, random_seed=42)
summary = pm.summary(trace)['mean'].to_dict()
for key in ['mean', 'amplitude', 'time-scale']:
    print(key, summary[key])

print('\nTraining score:')
print(np.log(np.mean(np.exp(trace.get_values('y_logp', chains=0)))) / T)
pm.traceplot(trace);
with model:
    y_future = pm.Poisson('y_future', mu=tt.exp(f[-F:]), observed=y_va)
    lik = pm.Deterministic('lik', y_future.logpt)
    logp_list = pm.sample_posterior_predictive(trace, vars=[lik], keep_size=True)
print('Heldout score:')
print(np.log(np.mean(np.exp(logp_list['lik'][0]))) / F)
with model:
    y_pred = pm.Poisson('y_pred', mu=tt.exp(f[T:]), shape=F)
    forecasts = pm.sample_posterior_predictive(trace, vars=[y_pred], random_seed=42)

samples = forecasts['y_pred']
low = np.zeros(F)
high = np.zeros(F)
mean = np.zeros(F)
median = np.zeros(F)

for i in range(F):
    low[i] = np.percentile(samples[:,i], 2.5)
    high[i] = np.percentile(samples[:,i], 97.5)
    median[i] = np.percentile(samples[:,i], 50)
    mean[i] = np.mean(samples[:,i])

xticks = np.arange(F)
plt.figure(figsize=(8,6))
plt.errorbar(xticks, median,
             yerr=[median-low, high-median],
             capsize=2, fmt='.', linewidth=1,
             label='2.5, 50, 97.5 percentiles');
plt.plot(xticks, mean, 'x', label='mean');
plt.plot(xticks, y_va, 's', label='observed');
plt.legend();
plt.title('Forecasts');
plt.xlabel('Day');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17841119028767485, 'fraction_numerical': 0.02454473475851148, 'mean_word_length': 4.963021243115657, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 2, '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': '265114', 'license': 'MIT', 'n_tokens_mistral': 3087, 'n_tokens_neox': 2880, 'n_words': 636, 'path': 'notebooks/GP_model.ipynb', 'repo_name': 'tufts-ml/single-hospital-count-forecasting', 'size': '7578'}