Прогнозиране на нови случаи на COVID19 в Португалия с помощта на процеси на Гаус

Практически уроци

Прогнозиране на нови случаи на COVID-19 в Португалия с помощта на процеси на Гаус

Използване на Python и Bayesian Statistics за прогнозиране на нови случаи за 30 дни

1. Въведение

Бих предпочел да направя този анализ на друга тема и в различен контекст. Днес Португалия е страната с най-голям абсолютен брой нови случаи на COVID-19 на един милион души [1]. Не искам случайно да превръщам тази статия в политическа статия, но има критична необходимост да се ограничи разпространението на вируса. Изисква се повече политическа смелост за предприемане на необходимите мерки и, от друга страна, по-стриктно поведение на населението, спазващо установените правила за локдаун.

Сега, след като премахнах тези мисли, нека се съсредоточим върху научната работа. Наборът от данни за COVID-19 вероятно е най-анализираният досега и са монтирани всякакви модели, за да се опитат да предскажат поведението му. Разпространението на пандемия има известна динамика, обикновено добре пригодена от отделни модели. Въпреки това в този случай имаме безпрецедентно глобално разпространение с много фактори, които го влияят. Например, държавите предприемат различни мерки за спиране на разпространението, някои от които с по-голямо въздействие върху броя на новите случаи, а други не изпитват много. Това е процес на учене чрез правене, като науката се използва колкото е възможно повече, за да помогне с политическите решения.

Обхватът на тази статия е върху прогнозата за броя на новите случаи в Португалия. Ще се опитаме да погледнем 30 дни в бъдещето и да видим какви различни сценарии могат да бъдат на масата за нас. За да преодолеем настоящия труден модел, ще използваме процеси на Гаус, за които е известно, че са изключително гъвкави.

Тази статия е научно упражнение; моля, не предприемайте никакви действия въз основа на неговите резултати.

Както обикновено, можете да намерите целия код в моя GitHub.

2. Процес на Гаус

Този раздел дава известна представа за теорията зад процесите на Гаус. Чувствайте се свободни да го пропуснете, ако се интересувате повече от приложението.

Процесът на Гаус е байесов непараметричен метод. Обикновено в байесовската рамка ние се интересуваме от извеждане на разпределение върху параметри на второ разпределение, което използваме, за да напаснем данните (вероятността). В този случай ние директно извеждаме разпределение върху функции.

По-прагматично погледнато, процесът на Гаус е случаен процес, при който присвояваме произволни променливи (или функции) f на краен брой дискретни точки за обучение, нашето X. Съвместното разпределение на тези променливи също е разпределение на Гаус:

където μ=(m(x_1),…,m(x_N)) са средните функции, K_ij=k(x_i,x_j) е ядрената функция, която дефинира ковариационната матрица и f=(f(x_1),…,f (x_N) са реализациите на стойностите на функциите.Например, ако ядрото счита x_i и x_j за подобни, може да се очаква стойностите на функциите в тези точки, f(x_i) и f(x_j) да имат подобни стойности.

2.1 Многомерно гаусово разпределение

Основата на процеса на Гаус е многомерното разпределение на Гаус, което е многомерното обобщение на разпределението на Гаус (вижте [2] към интерактивно обяснение). В многовариантния случай разпределението се определя от среден вектор μ и симетрична и положително определена ковариационна матрица Σ. Първият представлява очакваната стойност на всяка случайна променлива, докато последният описва два феномена: диагоналът изразява дисперсията на всяко измерение, а извън диагоналът ковариацията между всички случайни променливи. Той основно измерва как случайните променливи се променят заедно.

Многовариантното гаусово разпределение има следната съвместна плътност на вероятността:

където x е произволен вектор с размер d и |Σ| е детерминантата на dxd матрицата Σ.

За да разберем по-добре как тези параметри променят разпределението, нека начертаем два прости примера върху двумерно разпределение на Гаус. В първия пример променливите са независими, а във втория коварират (вижте стойностите извън диагонала).

Това може да се обобщи до:

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.gridspec as gridspec
import pandas as pd
import pymc3 as pm
import theano.tensor as tt
import arviz as az
from numpy.linalg import inv
from numpy.linalg import det
from scipy.optimize import minimize

pd.options.mode.chained_assignment = None
def multivariate_normal(x, d, mean, covariance):
    return (1. / (np.sqrt((2*np.pi)**d * np.linalg.det(covariance))) 
            * np.exp(-(np.linalg.solve(covariance, x-mean).T.dot(x-mean))/2))
def calculate_pdf(mean, covariance, d):
    x_ = 100
    xs = np.linspace(-5, 5, x_)
    ys = np.linspace(-5, 5, x_)
    x, y = np.meshgrid(xs, ys)
    pdf = np.zeros((x_, x_))
    for i in range(x_):
        for j in range(x_):
            pdf[i,j] = multivariate_normal(
                np.matrix([[x[i,j]], [y[i,j]]]), 
                d, mean, covariance)
    return x, y, pdf
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8,4), sharey=True)
d = 2 

mean_ind = np.matrix([[0.], [0.]])
covariance_ind = np.matrix([
    [1., 0.], 
    [0., 1.]])
x, y, p = calculate_pdf(
    mean_ind, covariance_ind, d)

ax[0].contourf(x, y, p, 100)
ax[0].set_xlabel('$X$', fontsize=13)
ax[0].set_ylabel('$Y$', fontsize=13)
ax[0].axis([-2.5, 2.5, -2.5, 2.5])

mean_corr = np.matrix([[0.], [0.]])
covariance_corr = np.matrix([
    [1., 0.8], 
    [0.8, 1.]])
x, y, p = calculate_pdf(
    mean_corr, covariance_corr, d)

con = ax[1].contourf(x, y, p, 100)
ax[1].set_xlabel('$X$', fontsize=13)
ax[1].set_ylabel('$Y$', fontsize=13)
ax[1].axis([-2.5, 2.5, -2.5, 2.5])

ax[0].set_title('Independent variables', fontsize=12)
ax[1].set_title('Correlated variables', fontsize=12)

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.02, 0.7])
cbar = fig.colorbar(con, cax=cbar_ax)
cbar.ax.set_ylabel('$p(x, y)$', fontsize=13)
plt.show()

Многовариантният Гаус е толкова полезен поради своите алгебрични свойства, като гарантира, че получаваме и гаусови разпределения при маргинализиране и кондициониране.

Маргинализацията означава интегриране на променливи от оригиналния набор от променливи. Резултатът е разпределението на подмножество от променливите без препратка към тези, които сме интегрирали. Така че, ако се интересуваме само от плътността на вероятността на X:

И както казах по-рано, това също води до разпределение на Гаус:

def univariate_normal(x, mean, variance):
    return ((1. / np.sqrt(2 * np.pi * variance)) * 
            np.exp(-(x - mean)**2 / (2 * variance)))
fig = plt.figure(figsize=(7, 7))
gs = gridspec.GridSpec(
    2, 1, height_ratios=[2, 1])
plt.suptitle('Marginal distributions', y=0.93)

ax1 = plt.subplot(gs[0])
x, y, p = calculate_pdf(mean_corr, covariance_corr, d)
con = ax1.contourf(x, y, p, 100)
ax1.set_xlabel('$x$', fontsize=13)
ax1.set_ylabel('$y$', fontsize=13)
ax1.axis([-2.5, 2.5, -2.5, 2.5])

ax2 = plt.subplot(gs[1])
x = np.linspace(-5, 5, 100)
px = univariate_normal(x, mean_corr[0,0], covariance_corr[0, 0])
ax2.plot(x, px, '--', label=f'$p(x)$')
ax2.legend(loc=0)
ax2.set_ylabel('density', fontsize=13)
ax2.yaxis.set_label_position('right')
ax2.set_xlim(-2.5, 2.5);

Същото се случва при кондициониране, което е вероятността една променлива да зависи от друга.

Нека изчислим p(X|y=1) за нашия пример:

y_cond = 1.
mean_x_cond_y = mean_corr[0,0] + (covariance_corr[0, 1]*(1/covariance_corr[1, 1])*(1 - mean_corr[1, 0]))
cov_x_cond_y = covariance_corr[0,0] - covariance_corr[0,1]*(1/covariance_corr[1,0]) * covariance_corr[0,1]
fig = plt.figure(figsize=(7, 7))
gs = gridspec.GridSpec(
    2, 1, height_ratios=[2, 1])
plt.suptitle('Conditional distributions', y=0.93)

ax1 = plt.subplot(gs[0])
x, y, p = calculate_pdf(mean_corr, covariance_corr, d)
con = ax1.contourf(x, y, p, 100)
ax1.set_xlabel('$x$', fontsize=13)
ax1.set_ylabel('$y$', fontsize=13)
ax1.axis([-2.5, 2.5, -2.5, 2.5])
ax1.plot([-2.5, 2.5], [y_cond, y_cond], color='lightblue', linestyle='dashed')

ax2 = plt.subplot(gs[1])
x = np.linspace(-5, 5, num=100)
px = univariate_normal(x, mean_x_cond_y, cov_x_cond_y)
ax2.plot(x, px, '--', label=f'$p(x|y=1)$')
ax2.legend(loc=0)
ax2.set_ylabel('density', fontsize=13)
ax2.yaxis.set_label_position('right')
ax2.set_xlim(-2.5, 2.5);

2.2 Гаусова регресия на процеса

Сега, когато разбираме многовариантните свойства на Гаус, е време да ги комбинираме, за да решим регресионни проблеми. Първо, трябва да въведем някои понятия.

Основната идея при извършване на регресия на процес на Гаус е, че нашият процес на Гаус ще присвои вероятност на потенциално безкрайните функции, които биха могли да паснат на нашите данни. Тази вероятност изразява несигурността на модела, което ни дава силна индикация за това доколко трябва да се доверим на точковата прогноза - т.е. средната стойност на полученото вероятностно разпределение.

2.3 Ядра

Ядрата дефинират типовете функции, които ще можем да извадим от нашето разпределение на функции. Въпреки че процесите на Гаус са непараметрични методи, ядрата имат параметри, които ни позволяват да контролираме формата на този специфичен тип функции. Нека разработим пример, за да разберем какво означава това.

Едно от най-използваните ядра е експоненциалното ядро на квадрат:

където дължината l контролира гладкостта, а η амплитудата на функцията.

az.style.use('arviz-darkgrid')
def sq_e_kernel(x1, x2, l=1.0, η=1.0):
    sqdist = np.sum(x1**2,1).reshape(-1,1) + np.sum(x2**2, 1) - 2 * np.dot(x1, x2.T)
    return η**2 * np.exp(-0.5 / l**2 * sqdist)
# Finite number of points
X = np.arange(0, 100).reshape(-1, 1)

cov = sq_e_kernel(X, X, l=5)
m = plt.imshow(cov, cmap="inferno", interpolation='none')
plt.colorbar(m);

Ако начертаем първия ред на ковариационната функция, можем да видим експоненциалното затихване въз основа на разстоянието между точките, което тази функция произвежда.

plt.plot(cov[:,0]);

2.4 Преди

Това, което дефинирахме по-горе, беше ковариационната матрица на нашето многовариантно разпределение. Това означава, че най-накрая можем да имаме достъп до реализации на нашите функции; трябва да вземем проби от нашето многомерно разпределение на Гаус. Можем да вземем проби, без реално да наблюдаваме каквато и да е тренировъчна точка. Това е изключително полезно при извършване на предварителни прогнозни проверки — ще използваме тази байесова функция, за да анализираме данните за COVID-19.

Нека дефинираме μ=0 за нашия процес на Гаус. Ще видим, че можем да моделираме всякакви видове поведение, използвайки само ковариационната матрица чрез различни ядра или комбинации от ядра. Стандартна практика е да оставим нашия μ=0.

samples = np.random.multivariate_normal(np.zeros(100), cov, size=3)
plt.plot(samples.T);

Можем да контролираме параметрите на ядрото или комбинация от ядра, за да променим формата на получените функции. По-долу можете да намерите 3 различни комбинации от параметри, използващи квадратно експоненциално ядро.

_, ax = plt.subplots(3, 2)
ax = ax.ravel()

l_p = [0.1, 2, 10]
η_p = [10, 2, 2]

for i in range(6):
    cov = sq_e_kernel(X, X, η=η_p[i//2], l=l_p[i//2])
    if not i%2:
        ax[i].imshow(cov, cmap="inferno", interpolation='none')
        samples = np.random.multivariate_normal(np.zeros(100), cov, size=3)
    else:
        ax[i].plot(samples.T)

2.5 Задна

За да завършим теоретичното въведение, трябва да включим нашите тестови точки, т.е. нови точки, където искаме да оценим нови функции. Това е мястото, където използваме условното свойство на многовариантния Гаус, което, както видяхме по-горе, също води до Гаусово разпределение.

В израза по-горе ние приемаме, че няма шум както в нашите данни за обучение, така и в прогнозите. Включването на шум в нашите данни за обучение се извършва чрез независимо добавяне на шум към всяко наблюдение,

За да вземем предвид шума в прогнозите, ние го добавяме към диагонала на K_**.

Ако в първия случай не се вземе предвид член на грешка, наборът от функции, произтичащи от кондиционирането, ще премине точно през всяка тренировъчна точка. В последния случай включването на шум може да се разглежда като добавяне на ядро с бял шум към нашето текущо, което помага да се моделират малки нередности в нашите прогнози.

def posterior(X_s, X_train, Y_train, l=1.0, η=1.0, noise_train=1e-8, noise_predict=1e-8):

    K = sq_e_kernel(X_train, X_train, l, η) + noise_train**2 * np.eye(len(X_train))
    K_s = sq_e_kernel(X_train, X_s, l, η)
    K_ss = sq_e_kernel(X_s, X_s, l, η) + noise_predict**2 * np.eye(len(X_s))
    K_inv = inv(K)
    
    mu_s = K_s.T.dot(K_inv).dot(Y_train)
    cov_s = K_ss - K_s.T.dot(K_inv).dot(K_s)
    
    return mu_s, cov_s

Нека използваме нашата последна функция, за да обусловим априора върху нови точки X_*. Първо, ние дефинираме нашите тренировъчни данни без никакъв шум и прогнози за прилягане също без да отчитаме никакъв шум.

X_train = np.arange(0,16).reshape(-1,1)
Y_train = np.sin(X_train)

X = np.linspace(0, 20, 100).reshape(-1, 1)
mu_new, cov_new = posterior(X, X_train, Y_train, noise_train=0, noise_predict=0)

samples = np.random.multivariate_normal(mu_new.ravel(), cov_new, 3)
uncertainty = 1.96 * np.sqrt(np.diag(cov_new))
plt.plot(X.ravel(),samples.T)
plt.fill_between(X.ravel(), mu_new.ravel() + uncertainty, mu_new.ravel() - uncertainty, alpha=0.1);

Когато добавим шум към нашите данни за обучение, но не добавим никакъв към нашия модел, можем да видим, че напълно преустройваме неправилното поведение на данните.

От друга страна, когато моделираме шума (на втория график), можем да видим, че можем да възстановим основния модел, присъстващ в данните, относително по-добре.

noise_train = 0.8
X_train = np.arange(0,16).reshape(-1,1)
Y_train = np.sin(X_train) + noise_train * np.random.randn(*Y_train.shape)

_, ax = plt.subplots(2, 1)
X = np.linspace(0, 20, 100).reshape(-1, 1)
mu_new, cov_new = posterior(X, X_train, Y_train, noise_train=0, noise_predict=0)

samples = np.random.multivariate_normal(mu_new.ravel(), cov_new, 3)
uncertainty = 1.96 * np.sqrt(np.diag(cov_new))
ax[0].plot(X.ravel(),samples.T)
ax[0].fill_between(X.ravel(), mu_new.ravel() + uncertainty, mu_new.ravel() - uncertainty, alpha=0.1);

mu_new, cov_new = posterior(X, X_train, Y_train, noise_train=0.2)
samples = np.random.multivariate_normal(mu_new.ravel(), cov_new, 3)
uncertainty = 1.96 * np.sqrt(np.diag(cov_new))
ax[1].plot(X.ravel(),samples.T)
ax[1].fill_between(X.ravel(), mu_new.ravel() + uncertainty, mu_new.ravel() - uncertainty, alpha=0.1);

2.6 Максимално a posteriori

Един от начините да оценим нашите параметри е чрез максимизиране на логаритмичната пределна вероятност (за да разберете как се извлича този израз, вижте [3]),

което е еквивалентно на минимизиране на отрицателната логаритмична пределна вероятност, която често се предпочита (вижте [4] за разширение на тези идеи).

def l_fn(X_train, Y_train):
    """Partially adapted from http://krasserm.github.io/2018/03/19/gaussian-processes/"""
    Y_train = Y_train.ravel()
    def negative_log_like(theta):
        K = sq_e_kernel(X_train, X_train, l=theta[0]) + theta[1]**2 * np.eye(len(X_train))
        return (0.5 * np.log(det(K))
                + 0.5 * Y_train.dot(inv(K).dot(Y_train))
                + 0.5 * len(X_train) * np.log(2*np.pi))
    return negative_log_like


res = minimize(l_fn(X_train, Y_train),
               [1.,1.],
               method='L-BFGS-B')

l_opt, noise_train_opt = res.x

print('l_opt:' + str(np.round(l_opt, 2)))
print('noise_train_opt:' + str(np.round(noise_train_opt, 2)))

mu_s, cov_s = posterior(X, X_train, Y_train, l=l_opt, noise_train=noise_train_opt)
plt.plot(X.ravel(), np.random.multivariate_normal(mu_s.ravel(), cov_s, 3).T)
plt.fill_between(X.ravel(), mu_s.ravel() + uncertainty, mu_s.ravel() - uncertainty, alpha=0.1);
l_opt:1.94
noise_train_opt:0.71

3. Модел

3.1 Данни

covid19 = pd.read_csv('./data/owid-covid-data20210120.csv')

covid = covid19[['location', 'continent', 'date', 'new_cases']] 
covid.columns = ['Country', 'Continent', 'Date', 'Count']
covid['Date'] = covid['Date'].astype('datetime64[ns]')
covid_pt = covid.loc[covid['Country']=='Portugal'].set_index('Date')
covid_pt = covid_pt.dropna()
dates = covid_pt.index
c_pt = covid_pt['Count'].fillna(0).reset_index().set_index('Date')['Count'].values
c_pt = np.abs(c_pt)
X = np.arange(c_pt.shape[0]).reshape(-1,1)

По-долу можем да видим развитието на новите случаи в Португалия. Както казах по-рано, разпространението на пандемията има известна динамика. Независимо от това, в този случай можем да видим, че са предприети действия с различни резултати във времето. В първата вълна на разпространение на вируса Португалия беше смятана за добър пример за ограничаването му. Това се дължи както на политически мерки, така и на поведението на населението като цяло. Същото не се случи във втората вълна, където липсата на планиране ни доведе до предизвикателна ситуация. В момента ние сме класирани като една от най-лошите страни в света, когато анализираме данните. Имаме работа с безпрецедентно нарастване на броя на случаите, което почти причинява колапс на нашата здравна система.

Данните включват стойности от 02–03–2020 г. до 19–01–2021 г. Не забравяйте, че стойностите винаги се отнасят за предходния ден, така че последният присъстващ ден в набора от данни всъщност е 18 януари.

n = c_pt.shape[0]

plt.plot(np.arange(n), c_pt, label='data')
plt.legend();

3.2 Изпълнение на модела

3.2.1 Ядра

Ядрата могат да бъдат класифицирани на стационарни и нестационарни (вижте другата ми статия, за да научите повече за тези концепции [5]). Стационарните ядра винаги ще се връщат към средната стойност на процеса на Гаус в области, които са по-далеч от наблюдаваните точки. Това се случва, защото ковариацията, която се моделира между точките, зависи само от тяхната относителна позиция, а не от тяхното абсолютно разстояние. Пример може да бъде степенно квадратично ядро или периодично ядро. От друга страна, линейното ядро е нестационарно ядро. Нашият модел е сбор от три GP за сигнала и един GP за шума:

Експоненциално квадратично ядро за моделиране на средносрочни нередности (различните вълни)
Периодичен термин за седмичната сезонност.
Линейно ядро за отчитане на тенденцията (експоненциална в този случай)
Шумът се моделира като ядро на бял шум.

Априорът на y като функция на времето е,

3.2.2 Средна функция

Ако не сме посочили средна функция в нашия модел, приемаме, че нашият GP има средна стойност нула. Ако използвахме само стационарни ядра, това би означавало, че функцията в крайна сметка ще се върне към нула, както прогнозираме в бъдещето. Това разумно предположение ли е? Вероятно случаят с тази пандемия не е такъв — ще трябва да живеем с нея в бъдеще. Предприетите мерки имат за основна цел да контролират скоростта на разпространение, а не точно да изчезнат от земята. Можем да дефинираме средна функция, различна от нула, и добавих кода, за да го направя по-долу. Въпреки това подходът, използван в нашия модел, ще използва вместо това линейно ядро.

from pymc3.gp.mean import Mean

class Linear(Mean):
    def __init__(self, b, a=0):
        Mean.__init__(self)
        self.a = a
        self.b = b

    def __call__(self, X):
        return tt.squeeze(tt.dot(X, self.b) + self.a)

3.3.3 Вероятност

В нашия случай имаме работа с данни за преброяване. Следователно, ние искаме да използваме един или повече GP като латентни процеси, които оценяват средната стойност на Поасон (и стандартното отклонение, тъй като те са еднакви).

Когато използваме разпределение като Поасон, трябва да сме сигурни, че параметърът, който оценяваме за разпределението, е ограничен до положителните реални числа. Тя може да бъде гарантирана чрез използване на конвенционалната функция за лог-линк — по този начин exp(f_i) винаги е положителен.

3.3.4 Използване на препараметризация за ускоряване

Добавям този малък раздел, за да ви дам видимост на параметризацията по подразбиране на PyMC3. Той използва нецентрирана параметризация на нашата многовариантна норма. Вече знаем, че извадките от Gaussian процес са реализации на многовариантна норма с параметризирана ковариационна матрица. Най-добрият начин да добавите многовариантна нормална плътност върху f към логаритмично-постериорната плътност е като използвате разлагането на Cholesky, за да получите коефициента L и го умножите по вектор от едномерни нормални случайни променливи.

Следователно предишната зависимост на плътността на f от неговите хиперпараметри (параметри на ядрото) се премахва.

3.3.5 Разпределения на хиперпараметри

Добра практика е винаги да се извършват някои предварителни прогнозни проверки и по този начин да се гарантира, че предишните разпределения на всички параметри са добре дефинирани. Поради гъвкавостта на процесите на Гаус, бих казал, че е още по-важно да го направим. Те могат лесно да надстроят данните и да ви дадат неочаквани резултати. Би било полезно, ако използвате информативни преди (вижте отлична статия по темата [6]).

X_ = np.linspace(0, 350,1000)[:, None]

with pm.Model() as model_prior:
    l_t = pm.InverseGamma('l_t', 4, c_pt.shape[0]/4)
    l_p = pm.InverseGamma('l_p', 4, c_pt.shape[0])
    c = pm.Normal('c', 0, 0.05)

    η_trend = pm.HalfNormal('η_trend',0.1)
    η_ts = pm.HalfNormal('η_ts', 0.01)
    η_per = pm.HalfNormal('η_per', 0.2)
    σ  = pm.HalfNormal("σ",  sigma=0.02)


    cov = (η_trend**2 * pm.gp.cov.ExpQuad(input_dim=1, ls=l_t)
           + η_ts**2 * pm.gp.cov.Linear(input_dim=1, c=c)
           + η_per**2 * pm.gp.cov.Periodic(1, period=7, ls=l_p)
            + pm.gp.cov.WhiteNoise(σ))
    
    
    gp = pm.gp.Latent(mean_func = pm.gp.mean.Zero(), cov_func = cov)
    
    f = gp.prior("f", X=X_)
    f_ = pm.Deterministic('f_', tt.exp(f))
    
    prior_checks = pm.sample_prior_predictive(samples=1000)
    
plt.plot(c_pt, color='darkorange')
plt.plot(X_, prior_checks['f_'].T, color='b', alpha=0.1)
plt.ylim(0, 15000)
plt.title('Draws from GP');

Можем да видим, че нашите предишни проби донякъде се съдържат в обхвата на нашите данни.

3.3.6 Всички заедно сега

По-долу можете да намерите модела, който използвахме, за да паснем на новите случаи на COVID-19 в Португалия, откакто пандемията достигна страната.

with pm.Model() as model:

    l_t = pm.InverseGamma('l_t', 4, c_pt.shape[0]/4)
    l_p = pm.InverseGamma('l_p', 4, c_pt.shape[0])
    c = pm.Normal('c', 0, 0.05)

    η_trend = pm.HalfNormal('η_trend',0.1)
    η_ts = pm.HalfNormal('η_ts', 0.01)
    η_per = pm.HalfNormal('η_per', 0.2)
    σ  = pm.HalfNormal("σ",  sigma=0.02)
    
    cov = (η_trend**2 * pm.gp.cov.ExpQuad(input_dim=1, ls=l_t)
           + η_ts**2 * pm.gp.cov.Linear(input_dim=1, c=c)
           + η_per**2 * pm.gp.cov.Periodic(1, period=7, ls=l_p)
            + pm.gp.cov.WhiteNoise(σ))
    
    
    gp = pm.gp.Latent(mean_func = pm.gp.mean.Zero(), cov_func = cov)
    f = gp.prior('f', X=X, reparameterize=True)


    y_pred = pm.Poisson('y_pred', mu=tt.exp(f), observed=c_pt)
    mp = pm.find_MAP(maxeval=20000, progressbar = True)

Сега е време да обусловим нашия модел с новите данни, т.е. 30-те нови времеви точки, които искаме да прогнозираме.

X_new = np.arange(c_pt.shape[0]+30).reshape(-1,1)
with pm.Model() as model:
    f_n = gp.conditional('f_n', Xnew=X_new)

    y_pred_new = pm.Poisson("y_pred_new", 
                            mu=tt.exp(f_n), 
                            shape=X_new.shape[0])

    pred_samples = pm.sample_posterior_predictive([mp], 
                                              vars=[y_pred_new], 
                                              samples=200,
                                              progressbar = False)

4. Резултати и дискусия

И накрая, нека начертаем нашите резултати.

from pymc3.gp.util import plot_gp_dist
fig = plt.figure(figsize=(12,5)); ax = fig.gca()

plot_gp_dist(ax, pred_samples['y_pred_new'], X_new, palette="Blues", fill_alpha=0.1, samples_alpha=0.1);
plt.plot(np.arange(c_pt.shape[0]),c_pt, label='data', color='darkorange')
plt.ylim(0,max(c_pt)*4)
plt.legend();

n_new = X_new.shape[0]

Увеличаването на предстоящите 30 дни, за които прогнозираме новините, не е добро. Моделът изчислява среден брой случаи, който достига до 24 047 на 18 февруари. Но по-важното е, че също така ни казва, че в случай на мерки като блокирането се окажат много ефективни, пикът от 24,047 може да бъде нисък от 8,753 (долната част на 50% интервал) или до 29,939 (горната част на 50% интервал). Можете също така да видите по-екстремните случаи по отношение на 90% надежден интервал.

Докато приключвах писането на тази статия, беше обявен броят на новите случаи за днес — добавих го като червен кръст. За съжаление можем да видим, че е в горния край на нашия интервал, между 50% и 90%. Може да е просто нередност (като резултат от повече тестове, регистър на данни, кумулативен ефект от предишни дни и т.н.) или да показва промяна в поведението. Ще трябва да анализираме данните от следващите няколко дни, за да го разберем.

mean_values = pd.DataFrame({'x':np.arange(n, n_new), 'y':np.mean(pred_samples['y_pred_new'], axis=0)[-30:]})
top_mean_values = pd.DataFrame(columns=['x', 'y'])
top_mean_values['x'] = [mean_values[(i)*6:(i+1)*6].sort_values(by='y', ascending=False)[:1].values[0][0] for i in range(5)]
top_mean_values['y'] = [mean_values[(i)*6:(i+1)*6].sort_values(by='y', ascending=False)[:1].values[0][1] for i in range(5)]

low_q25_values = pd.DataFrame({'x':np.arange(n, n_new), 'y':np.percentile(pred_samples['y_pred_new'], axis=0, q=[25]).ravel()[-30:]})
low_q25_values = low_q25_values.loc[low_q25_values.x.isin(top_mean_values.x)]
high_q25_values = pd.DataFrame({'x':np.arange(n, n_new), 'y':np.percentile(pred_samples['y_pred_new'], axis=0, q=[75]).ravel()[-30:]})
high_q25_values = high_q25_values.loc[high_q25_values.x.isin(top_mean_values.x)]

low_q10_values = pd.DataFrame({'x':np.arange(n, n_new), 'y':np.percentile(pred_samples['y_pred_new'], axis=0, q=[10]).ravel()[-30:]})
low_q10_values = low_q10_values.loc[low_q10_values.x.isin(top_mean_values.x)]
high_q10_values = pd.DataFrame({'x':np.arange(n, n_new), 'y':np.percentile(pred_samples['y_pred_new'], axis=0, q=[90]).ravel()[-30:]})
high_q10_values = high_q10_values.loc[high_q10_values.x.isin(top_mean_values.x)]
plt.fill_between(np.arange(n, n_new), np.percentile(pred_samples['y_pred_new'], axis=0, q=[10]).ravel()[-30:],
                 np.percentile(pred_samples['y_pred_new'], axis=0, q=[90]).ravel()[-30:], alpha = 0.1, color='b', label='90% CI');
plt.fill_between(np.arange(n, n_new), np.percentile(pred_samples['y_pred_new'], axis=0, q=[25]).ravel()[-30:],
                 np.percentile(pred_samples['y_pred_new'], axis=0, q=[75]).ravel()[-30:], alpha = 0.25, color='b', label='50% CI');
plt.plot(np.arange(n, n_new), np.mean(pred_samples['y_pred_new'], axis=0)[-30:], label='mean')
for i, j in zip(top_mean_values.x, top_mean_values.y):
    plt.annotate(str(int(j)), xy=(i,j))
for i, j in zip(low_q25_values.x, low_q25_values.y):
    plt.annotate(str(int(j)), xy=(i,j))
for i, j in zip(low_q10_values.x, low_q10_values.y):
    plt.annotate(str(int(j)), xy=(i,j))
for i, j in zip(high_q25_values.x, high_q25_values.y):
    plt.annotate(str(int(j)), xy=(i,j))
for i, j in zip(high_q10_values.x, high_q10_values.y):
    plt.annotate(str(int(j)), xy=(i,j))
    
plt.plot(324, 14647, 'rx', label='New cases Jan-19')
plt.legend();

5. Заключение

Тази статия има за цел да покаже силата и гъвкавостта на гаусовите процеси, когато се работи с модели, които са трудни за напасване. С такава мощност идва голяма отговорност, особено когато става въпрос за преоборудване. Когато се грижим да отделим време, за да изберем разумно нашите ядра и да дефинираме информативни приоритети за техните параметри, можем да съберем почти всякакъв вид данни с тях. Видяхме, че модел, толкова нередовен като португалските нови случаи на COVID-19 на ден, беше сравнително добре описан от нашия модел. Казвам сравнително добре, защото трябва да оценим нашите прогнози с показатели за грешки, да сравним с базовите модели и да приложим нашия модел към различни сценарии (други държави, други показатели като броя на смъртните случаи и т.н.). Друг важен аспект, който трябва да се има предвид при използване на гаусови процеси, е, че те не са известни като много мащабируеми (O(n³)), въпреки че има интересни подходи за увеличаване на тяхната мащабируемост (напр. [7]).

Тази статия е научно упражнение; моля, не предприемайте никакви действия въз основа на неговите резултати.

Бележка от редакторите на Towards Data Science: Въпреки че позволяваме на независими автори да публикуват статии в съответствие с нашите правила и насоки, ние не одобряваме приноса на всеки автор. Не бива да разчитате на творби на автор, без да потърсите професионален съвет. Вижте нашите Условия за четене за подробности.