Прогнозирование временных рядов с аппроксимацией параметрической кривой
Предскажите, когда остановится пандемия COVID-19 в вашей стране
Резюме
Люди сейчас находятся на карантине, гадая, когда закончится пандемия и жизнь вернется в нормальное русло. Основная цель этого руководства - понять, как получить данные о COVID-19 для вашей страны и спрогнозировать их распределение с помощью параметрической аппроксимации кривой.
В этой статье я шаг за шагом объясню как спрогнозировать распространение вируса в любой стране с помощью параметрической аппроксимации кривой. Я представлю полезный код Python, который можно легко использовать в других подобных случаях (просто скопировать, вставить, запустить) и пройтись по каждой строке кода с комментариями, чтобы вы могли легко воспроизвести этот пример (ссылка на полный код ниже).
Мы будем использовать наиболее часто используемый набор данных в дни карантина: набор данных CSSE COVID-19. Набор данных представляет собой временной ряд количества подтвержденных случаев заражения, сообщаемых каждой страной каждый день с момента начала пандемии. Этот набор данных находится в свободном доступе на github Университета Джона Хопкинса (ссылка ниже).
Полный код (Github):
Набор данных (Университет Джона Хопкинса):
Настраивать
Прежде всего, мы импортируем следующие библиотеки
## For data import pandas as pd import numpy as np ## For plotting import matplotlib.pyplot as plt ## For parametric fitting from scipy import optimize
Затем мы прочитаем данные в фрейм данных pandas.
dtf = pd.read_csv("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv", sep=",")
dtf.head()
В этом наборе данных каждая строка представляет собой временной ряд подтвержденных случаев в определенном географическом регионе. Я сгруппирую все на уровне страны и выберу только одну страну, я выберу Италия, поскольку она находится на пике заражения (вы можете выбрать любую другую страну или даже суммировать все на мировом уровне). Затем я создам новый столбец помимо столбца общее количество случаев: серию ежедневного увеличения общего количества, которое можно увидеть как количество новых случаев , рассчитывается как
новых случаев (t) = всего (t) - всего (t-1)
## group by country, sum and transpose
dtf = dtf.drop(['Province/State','Lat','Long'], axis=1).groupby("Country/Region").sum().T
## convert index to datetime
dtf.index = pd.to_datetime(dtf.index, infer_datetime_format=True)
## create total cases column
dtf = pd.DataFrame(index=dtf.index, data=dtf["Italy"].values, columns=["total"])
## create new cases column
dtf["new"] = dtf["total"] - dtf["total"].shift(1)
dtf["new"] = dtf["new"].fillna(method='bfill')
dtf.head()
dtf.tail()
Теперь у нас есть общее количество случаев в Италии и новые случаи за каждый день с 2020-01–22 до 2020–03–31 (сегодня), и они выглядят следующим образом:
А теперь приступим, ладно?
Дизайн модели
Модель является функцией независимой переменной и одного или нескольких коэффициентов (или параметров). Ошибка представляет собой случайные вариации данных, которые следуют определенному распределению вероятностей (обычно гауссовскому). Цель аппроксимации кривой - найти оптимальную комбинацию параметров, которая минимизирует ошибку. Здесь мы имеем дело с временными рядами, поэтому независимой переменной является время. С математической точки зрения:
y = f (время) + ошибка
Начнем с временного ряда общего числа случаев, нам нужно найти лучшую функцию для моделирования данных, а затем подобрать оптимальные параметры. Я попробую разные возможные модели, используя случайные коэффициенты, просто чтобы визуализировать кривые: линейная функция, экспоненциальная функция и логистическая функция.
'''
Linear function: f(x) = a + b*x
'''
def f(x):
return 10 + 1500*x
y_linear = f(x=np.arange(len(dtf)))
'''
Exponential function: f(x) = a + b^x
'''
def f(x):
return 10 + 1.18**x
y_exponential = f(x=np.arange(len(dtf)))
'''
Logistic function: f(x) = a / (1 + e^(-b*(x-c)))
'''
def f(x):
return 90000 / (1 + np.exp(-0.5*(x-20)))
y_logistic = f(x=np.arange(len(dtf)))
Теперь я построю общий временной ряд случая (черные точки) и 3 модели, определенные выше (цветные линии):
fig, ax = plt.subplots(figsize=(13,5)) ax.scatter(dtf["total"].index, dtf["total"].values, color="black") ax.plot(dtf["total"].index, y_linear, label="linear", color="red") ax.plot(dtf["total"].index, y_exponential, label="exponential", color="green") ax.plot(dtf["total"].index, y_logistic, label="logistic", color="blue") ax.legend() plt.show()
Казалось бы, экспоненциальная модель правильно соответствует данным… пока. Мы точно знаем, что у этого явления есть верхний предел, потому что вирус не может заразить больше, чем все население страны, поэтому рано или поздно рост остановится, и кривая будет плоской. Поэтому для этого больше подходит логистическая функция.
Давайте переключимся на новые временные ряды наблюдений: я попробую функцию Гаусса, используя некоторые случайные параметры только для визуализации.
'''
Gaussian function: f(x) = a * e^(-0.5 * ((x-μ)/σ)**2)
'''
def f(x):
return 6000 * np.exp(-0.5 * ((x-60)/6)**2)
y_gaussian = f(x=np.arange(len(dtf)))
Построим график фактических данных (черные полосы) и гауссовской модели, определенной выше (красная линия):
fig, ax = plt.subplots(figsize=(13,5)) ax.bar(dtf["new"].index, dtf["new"].values, color="black") ax.plot(dtf["new"].index, y_gaussian, color="red") plt.show()
Итак, мы решили продолжить:
- логистическая функция для моделирования общего временного ряда дел
- гауссовская функция для моделирования новых временных рядов дел.
Подобрать кривую
Пришло время выполнить подгонку, другими словами, мы собираемся найти оптимальные параметры (значения коэффициентов, которые минимизируют ошибку подгонки) для наших моделей. Мы будем использовать модуль optimize от scipy, который предоставляет функции для минимизации или максимизации целевых функций.
Давайте сначала подберем логистическую модель:
'''
Logistic function: f(x) = capacity / (1 + e^-k*(x - midpoint) )
'''
def logistic_f(X, c, k, m):
y = c / (1 + np.exp(-k*(X-m)))
return y
## optimize from scipy
logistic_model, cov = optimize.curve_fit(logistic_f,
X=np.arange(len(dtf["total"])),
y=dtf["total"].values,
maxfev=10000,
p0=[np.max(dtf["total"]), 1, 1])
## print the parameters
logistic_model
Затем мы можем сделать то же самое для гауссовской модели:
'''
Gaussian function: f(x) = a * e^(-0.5 * ((x-μ)/σ)**2)
'''
def gaussian_f(X, a, b, c):
y = a * np.exp(-0.5 * ((X-b)/c)**2)
return y
## optimize from scipy
gaussian_model, cov = optimize.curve_fit(gaussian_f,
X=np.arange(len(dtf["new"])),
y=dtf["new"].values,
maxfev=10000,
p0=[1, np.mean(dtf["new"]), 1])
## print the parameters
gaussian_model
Теперь, когда у нас есть подогнанные модели, мы, наконец, можем использовать их для прогнозирования.
Прогноз
До сих пор мы понимали, какие функции применять, и мы получили оптимальные параметры для ввода, другими словами, у нас есть 2 модели, одна для общих данных по случаям, а другая для данных ежедневного увеличения, и мы хотим предсказать будущее. . С этой целью мы применим эти две модели к новой независимой переменной: временным шагам от сегодняшнего дня до N. Для иллюстрации я сделаю прогноз на 30 дней вперед, начиная с сегодняшнего дня, поскольку наш набор данных уже содержит 69 временных шагов (строк), моя новая независимая переменная должна быть вектором в диапазоне от t = 70 до t = 100.
Чтобы получить хорошую визуализацию, я напишу полезную функцию для построения окончательных результатов:
'''
Plot parametric fitting.
'''
def utils_plot_parametric(dtf, zoom=30, figsize=(15,5)):
## interval
dtf["residuals"] = dtf["ts"] - dtf["model"]
dtf["conf_int_low"] = dtf["forecast"] - 1.96*dtf["residuals"].std()
dtf["conf_int_up"] = dtf["forecast"] + 1.96*dtf["residuals"].std()
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=figsize)
## entire series
dtf["ts"].plot(marker=".", linestyle='None', ax=ax[0], title="Parametric Fitting", color="black")
dtf["model"].plot(ax=ax[0], color="green")
dtf["forecast"].plot(ax=ax[0], grid=True, color="red")
ax[0].fill_between(x=dtf.index, y1=dtf['conf_int_low'], y2=dtf['conf_int_up'], color='b', alpha=0.3)
## focus on last
first_idx = dtf[pd.notnull(dtf["forecast"])].index[0]
first_loc = dtf.index.tolist().index(first_idx)
zoom_idx = dtf.index[first_loc-zoom]
dtf.loc[zoom_idx:]["ts"].plot(marker=".", linestyle='None', ax=ax[1], color="black",
title="Zoom on the last "+str(zoom)+" observations")
dtf.loc[zoom_idx:]["model"].plot(ax=ax[1], color="green")
dtf.loc[zoom_idx:]["forecast"].plot(ax=ax[1], grid=True, color="red")
ax[1].fill_between(x=dtf.loc[zoom_idx:].index, y1=dtf.loc[zoom_idx:]['conf_int_low'],
y2=dtf.loc[zoom_idx:]['conf_int_up'], color='b', alpha=0.3)
plt.show()
return dtf[["ts","model","residuals","conf_int_low","forecast","conf_int_up"]]
И последнее, но не менее важное: давайте напишем функцию для прогнозирования временного ряда:
'''
Forecast unknown future.
:parameter
:param ts: pandas series
:param f: function
:param model: list of optim params
:param pred_ahead: number of observations to forecast (ex. pred_ahead=30)
:param freq: None or str - 'B' business day, 'D' daily, 'W' weekly, 'M' monthly, 'A' annual, 'Q' quarterly
:param zoom: for plotting
'''
def forecast_curve(ts, f, model, pred_ahead=None, freq="D", zoom=30, figsize=(15,5)):
## fit
X = np.arange(len(ts))
fitted = f(X, model[0], model[1], model[2])
dtf = ts.to_frame(name="ts")
dtf["model"] = fitted
## index
index = pd.date_range(start=start,periods=pred_ahead,freq=freq)
index = index[1:]
## forecast
Xnew = np.arange(len(ts)+1, len(ts)+1+len(index))
preds = f(Xnew, model[0], model[1], model[2])
dtf = dtf.append(pd.DataFrame(data=preds, index=index, columns=["forecast"]))
## plot
utils_plot_parametric(dtf, zoom=zoom)
return dtf
Наконец, мы можем запустить его. Начнем, как обычно, с временных рядов общего числа обращений, а затем перейдем к временным рядам дневного прироста:
preds = forecast_curve(dtf["total"], logistic_f, logistic_model,
pred_ahead=30, freq="D", zoom=7)
preds = forecast_curve(dtf["total"], logistic_f, logistic_model,
pred_ahead=30, freq="D", zoom=7)
Согласно этим моделям, в Италии коронавирус уже замедляется, поскольку достигает максимальной способности заразиться, и в конце апреля общее количество случаев вырастет до 130 тысяч, а количество новых случаев упадет до нуля.
Заключение
Эта статья была руководством о том, как прогнозировать временной ряд с помощью параметрической подгонки кривой, в частности, мы взяли данные Covid-19 и сосредоточились на заражении Италии. Мы узнали, как обрабатывать данные для любой страны, как выбрать правильную модель, соответствующую данным, как найти оптимальные параметры и как использовать их для прогнозирования, когда пандемия COVID-19 остановится в выбранной стране.
Надеюсь, вам понравилось! Не стесняйтесь обращаться ко мне с вопросами и отзывами или просто для того, чтобы поделиться своими интересными проектами.
Эта статья является частью серии Прогнозирование временных рядов с помощью Python, см. Также:
