TowardsDataScience 博客中文翻译 2020（六百四十）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

模拟指数增长

原文：https://towardsdatascience.com/modeling-exponential-growth-49a2b6f22e1f?source=collection_archive---------2-----------------------

疾控中心在 Unsplash 拍摄的照片

使用对数变换、指数增长和线性回归预测冠状病毒传播

随着当前冠状病毒的爆发，我们听到了很多关于指数增长的说法。在本文中，我展示了如何理解和分析指数增长。如果你想跟随，你可以使用那些示例数据和一个简短的 Python 笔记本。

为什么是指数增长？

指数增长是一个数学函数，可用于多种情况。这个公式告诉我们某一时刻的病例数，对于冠状病毒来说，这是被感染的人数。

在指数增长的其他用例中，这个数字可能是动物数量的大小，也可能是你银行账户上的金额(如果你足够幸运有好的利率)。

使用指数增长来模拟冠状病毒爆发的原因是，流行病学家已经研究了这些类型的爆发，众所周知，流行病的第一阶段遵循指数增长。

指数增长公式

指数增长的特征在于以下公式:

指数增长函数

其中:

• x(t) 是任意给定时间 t 的病例数
• x0 是案例开始时的数量，也叫初始值
• b 是每个病人感染的人数，即增长系数

指数增长的一个简单例子:基数 2

为了更清楚地说明这一点，我将举一个假设的例子:

• 我们从 1 个感染者的初始值开始，所以 x0 = 1
• 每个病人会感染另外两个人，所以增长率 b = 2
• 我们将考察从时间 0 到时间 14 的疫情发展

我们首先需要将 a 和 b 的值代入公式，以获得特定流行病的公式:

然后我们可以用这个公式来计算从 0 到 14 的每个 t 值的 y 值。当我们这样做时，我们在每个时间点获得了以下感染人数，如下表所示。这表明，从 1 个人开始，每个人的增长因子为 2，我们在 14 天后获得了超过 16000 个案例。

如果我们想用图表来表示，我们开始看到一个图表，看起来很像我们看到的关于冠状病毒的非常惊人的曲线:

生长因子为 2 的指数生长图

找到冠状病毒的确切分子式

现在，我们知道这个图或多或少有正确的形状，但我们需要做一个额外的步骤，使我们的分析有用。我们需要通过查看疫情传播的数据，找到电晕疫情的真正增长因素。

冠状病毒数据的摘录。来源:https://covid.ourworldindata.org/data/full_data.csv

寻找生长因子的线性回归

在查看数据时，我们只有每天的病例数，而没有增长因素。从每日的经验观察中找到生长因子的最好方法是使用一个叫做线性回归的统计模型。

给定 y 和 x 的经验观察值，线性回归允许我们在下面的公式中估计 a 和 b 的最佳值。在该公式中，y 是病例数，x 是时间。但是我们需要对指数增长函数进行一些重写，因为线性回归只能估计如下所示的公式:

线性回归所需的公式类型

重写线性回归的指数公式

首先，我们需要以线性回归的形式重写公式。我们需要的工具是对数。对数允许以正确的形式重写函数:

• 我们使用感染数量的记录而不是感染数量
• 我们使用生长因子的对数，而不是生长因子****

将线性回归应用于数据

**步骤 1—**Python 记事本的第一步是导入数据并应用日志转换:

步骤 2 — 然后我们使用 statsmodels 库来估计线性回归函数:

步骤 3 —根据表格制作预测函数

让我们回到线性回归的公式:

线性回归所需的公式类型

statsmodels 表给出了 coef (中间) : 下 a 和 b 的值

• 值 const 是我们线性回归中 a 的值: 0.4480
• 值时间是我们的线性回归中 b 的值:0.1128

因此，我们现在可以填写线性回归函数。

请记住:

• 我们说线性回归的一个 就是的对数的初始值
• 我们说线性回归的 b 是的对数增长因子

因此，我们知道:

• 初始值的对数等于 0.4480
• 生长因子的对数是 0.1128

为了找到实际的值，我们需要通过应用指数来“解开”它们。

现在，我们可以回到指数增长的原始公式，并填入这些值，以找出冠状病毒案例的实际值:

冠状病毒流行的实际公式！

对两周以后做预测！

现在我们已经估算了这个新公式，我们可以用它来预测任何我们想要的日期。值得注意的是，这里的预测只是一个例子，说明数学和统计学如何应用于流行病学。现实生活中的流行病学家将测试除指数增长之外的不同类型的模型，并在模型验证方面做大量工作，而目前的例子还没有这样做。

一旦找到最佳模型，就可以用它来进行预测。使用我们使用指数增长曲线估计的函数，如果我们想预测数据集最后一天(第 68 天)后的 2 周，我们只需将 t= 68 放入公式中，模型就预测当天有 3355 例感染。

进一步发展的可能性:

我已经展示了如何应用线性模型来预测指数增长过程。需要采取一些预防措施:

• 线性模型只是指数增长函数的最佳估计，它有一定的误差，我们可以在进一步的研究中检验
• 指数增长函数不一定是流行病的完美代表。我已经找到了最合适的指数增长函数，但是下一步要研究的可能是逻辑增长
• 指数级的增长只会在一开始符合流行病。在某种程度上，治愈的人将不再传播病毒，当(几乎)每个人都被感染或已经被感染时，增长将停止。

在这里，你可以找到一篇关于冠状病毒的逻辑增长的文章 ，它也考虑到了疫情的最后阶段。

感谢您阅读本文。我希望你喜欢它。不要犹豫，继续关注更多！

建模功能

原文：https://towardsdatascience.com/modeling-functions-78704936477a?source=collection_archive---------13-----------------------

从线性回归到逻辑回归

Clem Onojeghuo 在 Unsplash 上拍摄的照片

**Table of Contents**[**Introduction**](#3846)1\. [Linear models](#e296)
2\. [Quadratic models](#005e)
3\. [Cubic models](#35bc)
4\. [Exponential models](#7b2b)
5\. [Logarithmic models](#ec20)
6\. [Sinusoidal models](#48d3)
7\. [Logistic models](#8108)[**Conclusion**](#9ef6)

介绍

我们将使用 Jupyter Notebook 绘制一个散点图，并建立一个从线性到逻辑的回归线模型。

线性模型

第一个是线性模型。线性模型被表示为𝑦=𝑚𝑥+𝑐.我们将使用 numpy.array 或 numpy.arange 来创建数据。如果你想了解更多关于线性关系的内容，请阅读线性关系的衡量标准。我们导入 Python 库 numpy 和 matplotlib。我们创建了一个年份和一个二氧化碳阵列。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlineyear=np.array([1980,1982,1984,1986,1988,1990,1992,1994,1996,1998,2000])
co2=np.array([338.7,341.1,344.4,347.2,351.5,354.2,356.4,358.9,362.6,366.6,369.4])

首先，我们使用 matplotlib 创建一个散点图。添加标题、标签、x 轴和 y 轴标签。你需要使用show()方法。您可以在没有它的情况下进行打印，但这将删除不必要的输出。

plt.scatter(year,co2,label='CO2')
plt.title("Year vs CO2")
plt.xlabel('Year')
plt.ylabel('CO2')
plt.legend()
plt.show()

x 轴上的整数

如上图所示，x 轴上有小数。在下面的代码中，我们使用前三行使它们成为整数。

from matplotlib.ticker import MaxNLocatorax = plt.figure().gca()
ax.xaxis.set_major_locator(MaxNLocator(integer=True))plt.scatter(year,co2,label='CO2')
plt.title("Year vs CO2")
plt.xlabel('Year')
plt.ylabel('CO2')
plt.legend()
plt.show()

x 轴上的整数

用numpy.polyfit和numpy.poly1d寻找线性模型

最简单的方法就是用numpy.polyfit。通过将 order 设置为 1，它将返回一个线性系数数组。在numpy.poly1d中使用它会返回一个使用系数的等式。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from sklearn.linear_model import LinearRegression
%matplotlib inlineax = plt.figure().gca()
ax.xaxis.set_major_locator(MaxNLocator(integer=True))year=np.array([1980,1982,1984,1986,1988,1990,1992,1994,1996,1998,2000])
co2=np.array([338.7,341.1,344.4,347.2,351.5,354.2,356.4,358.9,362.6,366.6,369.4])coef = np.polyfit(year, co2, 1)
equ = np.poly1d(coef)x_plot = np.linspace(1975,2005,100)
y_plot = equ(x_plot)
plt.plot(x_plot, y_plot, color='r')plt.scatter(year,co2,label='CO2')
plt.title("Year vs CO2")
plt.xlabel('Year')
plt.ylabel('CO2')
plt.legend()
plt.show()

散点图和线性回归线

使用 scikit-learn 查找线性模型

求回归斜率和截距的第二种方法是使用[sklearn.linear_model.LinearRegression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)。该类要求 x 值为一列。我们使用reshape(-1,1)修改年份数据。原始年份数据具有 1x 11 形状。您需要将年份数据调整为 11 乘 1。

year1=year.reshape((-1,1))
print(np.shape(year))
print(np.shape(year1))

我们导入[sklearn.linear_model.LinearRegression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)，重塑年份数据，使用LinearRegression().fit()拟合我们的数据。这将返回斜率coef_和 y 轴截距intercept_。coef_返回一个数组，所以我们用reg.coef_[0]取第一项。让我们打印出我们的回归线方程。

线性方程

from sklearn.linear_model import LinearRegressionyear1=year.reshape((-1,1))reg = LinearRegression().fit(year1,co2)slope=reg.coef_[0]
intercept=reg.intercept_print(f'The equation of regression line is y={slope:.3f}x+{intercept:.3f}.')

一起

我们一起画一个散点图和我们的线性回归线。我们使用从 1975 年到 2005 年的新 x 域，取 100 个样本作为回归线，np.linspace(1975,2005,100)。然后使用 x 域、斜率和 y 截距绘制回归线。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from sklearn.linear_model import LinearRegression
%matplotlib inlineax = plt.figure().gca()
ax.xaxis.set_major_locator(MaxNLocator(integer=True))year=np.array([1980,1982,1984,1986,1988,1990,1992,1994,1996,1998,2000])
co2=np.array([338.7,341.1,344.4,347.2,351.5,354.2,356.4,358.9,362.6,366.6,369.4])year1=year.reshape((-1,1))reg = LinearRegression().fit(year1,co2) 
slope=reg.coef_[0]
intercept=reg.intercept_plt.scatter(year,co2,label='CO2')
X_plot = np.linspace(1975,2005,100)
Y_plot = slope*X_plot+intercept
plt.plot(X_plot, Y_plot, color='r')
plt.title("Year vs CO2")
plt.xlabel('Year')
plt.ylabel('CO2')
plt.legend()
plt.show()print(f'The equation of regression line is y={slope:.3f}x+{intercept:.3f}.')

用 scipy 寻找线性模型

另一种寻找回归斜率和截距的方法是使用[scipy.stats.linregress](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html)。这将返回slope, intercept, rvalue, pvalue, stderr。

from scipy.stats import linregressslope, intercept, r_value, p_value, std_err = linregress(year,co2)
print(f'The equation of regression line is y={slope:.3f}x+{intercept:.3f}.')

线性回归图

为了画一条线，我们需要 x 个点。我们使用np.linspace，它是numpy.linspace，因为我们使用了import numpy as np。我们的数据是从 1975 年到 2000 年。所以我们用 1960 代表start，2005 代表stop，100 代表样本数。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress
%matplotlib inlineyear=np.array([1980,1982,1984,1986,1988,1990,1992,1994,1996,1998,2000])
co2=np.array([338.7,341.1,344.4,347.2,351.5,354.2,356.4,358.9,362.6,366.6,369.4])X_plot = np.linspace(1975,2005,100)
Y_plot = slope*X_plot+intercept
plt.plot(X_plot, Y_plot, color='r')
plt.show()

使用 scipy . Lin regression 的线性回归线

现在我们把散点图、回归线和回归方程放在一起。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from scipy.stats import linregress%matplotlib inlineax = plt.figure().gca()
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
slope, intercept, r_value, p_value, std_err = linregress(year,co2)
X_plot = np.linspace(1975,2005,100)
Y_plot = slope*X_plot+intercept
plt.plot(X_plot, Y_plot, color='r')
plt.scatter(year,co2,label='CO2')
plt.title("Year vs CO2")
plt.xlabel('Year')
plt.ylabel('CO2')
plt.legend()
plt.show()print(f'The equation of regression line is y={slope:.3f}x+{intercept:.3f}.')

散点图和线性回归线

练习 1

使用以下数据绘制散点图和回归线。求一个线性回归方程。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinetemp = np.array([55,60,65,70,75,80,85,90])
rate = np.array([45,80,92,114,141,174,202,226])

回答

你画了散点图和回归图吗？回归线应该是𝑦= 5.119𝑥—236.88。

[## 斯皮尔曼等级相关系数的微妙性

单调关系的未知部分

towardsdatascience.com](/the-subtlety-of-spearmans-rank-correlation-coefficient-29478653bbb9)

二次模型

我们使用 Numpy 的arange来创建从 0 到 9 的 10 个整数。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinetime = np.arange(10)
height = np.array([450,445,430,409,375,331,280,215,144,59])

我们来绘制上面的数据。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlineplt.scatter(time,height,label='Height of a ball')
plt.title("Time vs Height")
plt.xlabel('Time')
plt.ylabel('Height')
plt.legend()
plt.show()

二次回归图

numpy.polyfit拟合多项式。它需要 x，y 和拟合多项式的次数。二次为 2，三次为 3，以此类推。它返回一个数组，该数组的多项式系数为常数的高次幂。对于二次函数，它们是 a、b 和 c:

coef = np.polyfit(time, height, 2)
coef

让我们打印出二次回归线。

print(f'The equation of regression line is y=')
print(equ)

或者使用系数，回归线是:

print(f'The equation of regression line is y={coef[0]:.3f}x^2+{coef[1]:.3f}x+{coef[2]:.3f}.')

二次回归方程

我们再次使用 NumPy 的poly1d和polyfit。np.poly1d(coefficients)将使用我们的系数返回一个多项式方程。

equ = np.poly1d(coef)

我们可以找到任意 x 的值。例如，如果您想在 x=1 时找到 y 值:

equ(1)

x=1 时的 y 值

我们用这个来画回归线。我们使用numpy.linspace为 100 个样本定义从 0 到 10 的 x 值。并在equ中使用它作为 y 值。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinex_plot = np.linspace(0,10,100)
y_plot = equ(x_plot)
plt.plot(x_plot, y_plot, color='r')
plt.show()

我们把它们放在一起。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinetime = np.arange(10)
height = np.array([450,445,430,409,375,331,280,215,144,59])coef = np.polyfit(time, height, 2)
equ = np.poly1d(coef)x_plot = np.linspace(0,10,100)
y_plot = equ(x_plot)
plt.plot(x_plot, y_plot, color='r')plt.scatter(time,height,label='Height of a ball')
plt.title("Time vs Height")
plt.xlabel('Time')
plt.ylabel('Height')
plt.legend()
plt.show()print(f'The equation of regression line is y=')
print(equ)

练习 2

通过使用以下数据，在图表中绘制散点图和回归线。求二次回归方程。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlineangle = np.arange(20,80,10)
distance = np.array([371,465,511,498,439,325])

回答

你能画一条散点线和回归线吗？二次方程应该是:

二次回归练习答案

[## 卡方独立性检验简介

使用 Jupyter 笔记本的卡方初学者指南

towardsdatascience.com](/gentle-introduction-to-chi-square-test-for-independence-7182a7414a95)

立方模型

可以用和上面二次函数一样的方法。我们要用plyfit和poly1d。首先，我们准备数据。让我们画一个散点图。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlineengspeed = np.arange(9,23,2)
avespeed = np.array([6.45,7.44,8.88,9.66,10.98,12.56,15.44])plt.scatter(engspeed,avespeed,label='Speed of different boat engine')plt.title("Average speed of different boat engine")
plt.xlabel('Engine speed')
plt.ylabel('Boad speed')
plt.ylim(0,20)
plt.legend()
plt.show()

使用polyfit返回系数。对于三次函数，a、b、c 和 d 在:

立方函数

coef = np.polyfit(engspeed, avespeed, 3)
print(coef)

三次函数的系数

我们把所有的放在一起。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlineengspeed = np.arange(9,23,2)
avespeed = np.array([6.45,7.44,8.88,9.66,10.98,12.56,15.44])plt.scatter(engspeed,avespeed,label='Speed of different boat engine')coef = np.polyfit(engspeed, avespeed, 3)
equ = np.poly1d(coef)x_plot = np.linspace(8,25,100)
y_plot = equ(x_plot)
plt.plot(x_plot, y_plot, color='r')plt.title("Average speed of different boat engine")
plt.xlabel('Engine speed')
plt.ylabel('Boad speed')
plt.ylim(0,20)
plt.legend()
plt.show()a, b, c, d = coef
print(f'The equation of regression line is y={a:.3f}x^3+{b:.3f}x^2+{c:.3f}x+{d}.')

散点图和三次回归线

练习 3

使用以下数据绘制散点图和三次回归线。打印三次方程。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinex=np.arange(1,8)
y=np.array([0,0.012,0.06,0.162,0.336,0.6,0.972])

回答

你能画出散点图和回归线吗？回归线方程应为:

系数为[3.000000000 e-03，-1.16796094e-16，-9.000000000 e-03，6.00000000e-03]。这些意味着:

第二个实际上是 0。尝试以下方法，看看两者是否都是相同的 0.3。

print(300e-03)
print(300*10**(-3))

指数模型

我们将探索三种指数模型。

让我们设置数据。画一个散点图。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlineday = np.arange(0,8)
weight = np.array([251,209,157,129,103,81,66,49])plt.scatter(day,weight,label='Weight change')
plt.title("Day vs Weight")
plt.xlabel('Day')
plt.ylabel('Weight')
plt.legend()
plt.show()

我们要用[scipy.optimize.curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)。这需要一个函数，x 和 y 值，和初始值，p0以数组的形式。找到合适的p0需要一点反复试验。你必须测试不同的价值观。我们用p0=(1, 1e-6, 1)。它返回参数的最优值和 popt 的估计协方差。

𝑎⋅𝑒^−𝑏𝑥+𝑐

我们的第一个指数函数使用 a、b 和 c。我们将首先定义一个函数。这在curve_fit方法中使用。对于平滑曲线，我们用 100 个样本使用numpy.linspace从 0 到 7 设置 x 值。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inlinedef func(x, a, b, c):
    return a * np.exp(-b * x) + cpopt, pcov = curve_fit(func, day, weight, p0=[1, 1e-6, 1])x_plot=np.linspace(0,7,100)
plt.plot(x_plot, func(x_plot, *popt), 'r-')plt.scatter(day,weight,label='Day vs Weight')
plt.title("Day vs Weight a*e^-bx +c")
plt.xlabel('Day')
plt.ylabel('Weight')
plt.legend()
plt.show()# equation
a=popt[0].round(2)
b=popt[1].round(2)
c=popt[2].round(2)print(f'The equation of regression line is y={a}e^({b}x)+{c}')

散点图和指数回归线

𝑎⋅𝑒^−𝑏𝑥

第二个函数使用 a 和 b。我们相应地定义函数。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inlinedef func2(x, a, b):
    return a * np.exp(-b * x)popt, pcov = curve_fit(func2, day, weight, p0=[1, 1e-6])x_plot=np.linspace(0,7,100)
plt.plot(x_plot, func2(x_plot, *popt), 'r-')plt.scatter(day,weight,label='Day vs Weight')
plt.title("Day vs Weight a*e^-bx")
plt.xlabel('Day')
plt.ylabel('Weight')
plt.legend()
plt.show()# equation
a=popt[0].round(2)
b=popt[1].round(2)print(f'The equation of regression line is y={a}e^({b}x)')

指数回归线的第二个例子

𝑎⋅𝑏^𝑥

最后一个指数函数使用 a 和 b，我们相应地修改函数。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inlinedef func3(x, a, b):
    return a * b ** xpopt, pcov = curve_fit(func3, day, weight, p0=[1, 1e-6])
x_plot=np.linspace(0,7,100)
plt.plot(x_plot, func3(x_plot, *popt), 'r-')plt.scatter(day,weight,label='Day vs Weight')
plt.title("Day vs Weight a*b^x")
plt.xlabel('Day')
plt.ylabel('Weight')
plt.legend()
plt.show()# equation
a=popt[0].round(4)
b=popt[1].round(4)print(f'The equation of regression line is y={a}*{b}^x')

与 TI Nspire 结果比较

TI Nspire 的指数回归使用转换值 x 和𝑙𝑛的最小二乘拟合，将模型方程 𝑦 = 𝑎𝑏^𝑥 拟合到数据上*(𝑦)。它返回不同的值。*

实践

使用下列数据找出 ab^x.形式的指数函数，绘制散点图并绘制回归线。

import numpy as npweek = np.arange(1,21)
views = np.array([102365, 38716,21617,24305,9321,14148,2103,8285,5098,3777,831,1007,834,34,378,204,6,42,54,31])

回答

对数模型

半对数模型

通常，我们对指数函数使用半对数模型:

我们设置了模拟数据并绘制了散点图。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinetime = np.arange(0,30,4)
bacteria = np.array([20,150,453,920,1820,9765,15487,19450])plt.scatter(time,bacteria,label='Bacteria')
plt.title("Time vs Bacteria")
plt.xlabel('time')
plt.ylabel('bacteria')
plt.legend()
plt.show()

时间对细菌图

我们将使用[numpy.log](https://docs.scipy.org/doc/numpy/reference/generated/numpy.log.html)对细菌值进行自然记录。numpy.log是自然对数。这应该显示一个线性趋势。我们需要用ln(bacteria)修改标题和 y 标签。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinetime = np.arange(0,30,4)
bacteria = np.array([20,150,453,920,1820,9765,15487,19450])plt.scatter(time,np.log(bacteria),label='Bacteria')
plt.title("Time vs ln(Bacteria)")
plt.xlabel('time')
plt.ylabel('ln(bacteria)')
plt.legend()
plt.show()

我们使用在二次和三次函数中使用的numpy.polyfit。我们在numpy.polyfit()中使用1，这样它将返回一个线性回归。numpy.polyfit返回等式的所有系数。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinetime = np.arange(0,30,4)
bacteria = np.array([20,150,453,920,1820,9765,15487,19450])p = np.polyfit(time, np.log(bacteria), 1)
plt.plot(time, p[0] * time + p[1], 'g--', label='Semi-log graph')plt.scatter(time,np.log(bacteria),label='Bacteria')
plt.title("Time vs Bacteria")
plt.xlabel('time')
plt.ylabel('bacteria')
plt.legend()
plt.show()print(f'The equation of regression line is y={p[0]:.3f} * x + {p[1]:.3f}')

双对数模型

双对数模型用于幂函数。

让我们设置数据并绘制散点图。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinex=np.array([2,30,70,100,150])
y=np.array([4.24,16.4,25.1,30,36.7])plt.scatter(x,y,label='Log-log')
plt.title("Log-Log model")
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

我们使用[numpy.log](https://docs.scipy.org/doc/numpy/reference/generated/numpy.log.html)获取 x 和 y 值的自然对数。我们需要将 x 和 y 标签修改为 ln(x)和 ln(y)。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinex=np.array([2,30,70,100,150])
y=np.array([4.24,16.4,25.1,30,36.7])p = np.polyfit(np.log(x), np.log(y), 1)
plt.plot(np.log(x), p[0] * np.log(x) + p[1], 'r--', label='Regression line')plt.scatter(np.log(x),np.log(y),label='log-log')
plt.title("Log-log regression")
plt.xlabel('ln(x)')
plt.ylabel('ln(y)')
plt.legend()
plt.show()print(f'The equation of regression line is ln(y)={p[0]:.3f} * ln(x) + {p[1]:.3f}')

正弦模型

让我们试试正弦函数。我们设置数据并绘制散点图。既然要用scipy.optimize.curve_fit，那我们也导入一下吧。我们在指数模型中使用了它。我们建立了数据并绘制了散点图。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inlineyear=np.arange(0,24,2)
population=np.array([10.2,11.1,12,11.7,10.6,10,10.6,11.7,12,11.1,10.2,10.2])plt.scatter(year,population,label='Population')
plt.title("Year vs Population")
plt.xlabel('Year')
plt.ylabel('Population')
plt.legend()
plt.show()

我们定义一个叫做sinfunc的函数。这需要参数x, a, b, c, d。我们用[numpy.radians](https://docs.scipy.org/doc/numpy/reference/generated/numpy.radians.html)表示c。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inlineyear=np.arange(0,24,2)
population=np.array([10.2,11.1,12,11.7,10.6,10,10.6,11.7,12,11.1,10.2,10.2])def sinfunc(x, a, b, c, d):
    return a * np.sin(b * (x - np.radians(c)))+dpopt, pcov = curve_fit(sinfunc, year, population, p0=[1,0.4,1,5])
x_data = np.linspace(0, 25, num=100)plt.scatter(year,population,label='Population')
plt.plot(x_data, sinfunc(x_data, *popt), 'r-',label='Fitted function')
plt.title("Year vs Population")
plt.xlabel('Year')
plt.ylabel('Population')
plt.legend()
plt.show()a, b, c, d = poptprint(f'The equation of regression line is y={a:.3f} * sin({b:.3f}(x-{np.radians(c):.3f}))+{d:.3f}')

实践

使用下面的表格，画一个散点图并找到一个余弦回归函数。

回答

你可能有不同的系数。我用过

逻辑模型

我们收集数据并绘制散点图。我们使用plt.xlim和plt.ylim设置域从-10 到 10，范围从 0 到 250。我们将使用scipy.optimize.curve_fit进行逻辑回归。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inlinex=np.arange(0,10)
y=np.array([52,133,203,230,237,239.5,239.8,239.9,240,240])plt.scatter(x, y, label='Regression line')
plt.title("Logistic regression")
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-10,10)
plt.ylim(0,250)
plt.legend()
plt.show()

我们使用logifunc来定义我们的逻辑功能。我们用curve_fit找到popt中的函数参数。对于回归线，我们为函数设置一个新的域，x_data从-10 到 10。我们使用plt.plot绘制这条线。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inlinex=np.arange(0,10.0)
y=np.array([52,133,203,230,237,239.5,239.8,239.9,240,240])def logifunc(x,L,c,k):
    return L/ (1 + c*np.exp(-k*x))popt, pcov = curve_fit(logifunc, x, y, p0=[200,1,1])
x_data = np.linspace(-10, 10, num=100)plt.scatter(x,y,label='Logistic function')
plt.plot(x_data, logifunc(x_data, *popt), 'r-',label='Fitted function')
plt.title("Logistic")
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-10,10)
plt.ylim(0,250)
plt.legend()
plt.show()

当 y 数据点为负时

有时，您的数据在 y 坐标中可能有负值。

import pandas as pddf = pd.read_csv('[http://bit.ly/2tUIZjK'](http://bit.ly/2tUIZjK'))df.head()

数据的最小值必须为零。理想情况下，乙状结肠中点也为零。但是上面的数据集两者都不满足。使用等式(1–2)，并添加offset值适用于该数据集。

x=df.T.iloc[0]
y=df.T.iloc[1]def logifunc(x,L,x0,k,off):
    return L / (1 + np.exp(-k*(x-x0)))+offplt.scatter(x,y,label='Logistic function')popt, pcov = curve_fit(logifunc, x, y, p0=[50,185,0.1,-222])
plt.scatter(x,y,label='Logistic function')x_data = np.linspace(170,205,num=100)
plt.plot(x_data, logifunc(x_data, *popt), 'r-',label='Fitted function')
plt.legend()plt.show()
print(popt)

结论

scipy.optimize.curve_fit对许多功能都有用。唯一的问题是在p0中找到好的初始值。有时不同的p0值会返回不同的popt。可以试试 LMFIT 。

通过 成为 的会员，可以完全访问媒体上的每一个故事。

https://blog.codewithshin.com/subscribe

参考

• https://realpython.com/linear-regression-in-python/
• https://stack overflow . com/questions/41925157/logisticregression-unknown-label-type-continuous-using-sk learn-in-python

用图形数据库对医疗保健数据建模

原文：https://towardsdatascience.com/modeling-healthcare-data-with-graph-databases-3e3695bcae3c?source=collection_archive---------27-----------------------

使用 TigerGraph 和 Synthea 创建合成医疗保健系统

弗兰基·查马基在 Unsplash 上拍摄的照片

介绍

自从从纸质记录过渡到虚拟记录以来，医院一直在堆积数据。医疗保健系统的每个接触点、每个处方、手术和免疫都记录并存储在医院的电子健康记录中(EHR)。现在已经到了医院数据多到不知如何处理的地步。更糟糕的是，这种复杂数据的过饱和使得访问和分析数据的效率非常低。

那么，这场危机的解决方案是什么？

图形数据库！

图表非常适合存储和可视化医疗保健数据。它们旨在处理高度关联的信息，如病历。如果你对图形不熟悉，看看这篇的精彩文章，它介绍了一些图论的基础知识。

我将通过一个例子展示我们如何使用一个 TigerGraph 图形数据库来表示复杂的医疗保健数据。我们开始吧！

生成数据

在一个理想的世界里，我们可以用真实的病人数据来创建这个图表；然而，有许多规则和条例使得处理患者数据相当困难。而是可以用退而求其次的东西:合成数据。

使用 Synthea ，一个开源的合成病人生成器，我们可以创建一个完整的医疗保健生态系统，其中充满了病人、医院访问、保险提供商以及你能想到的一切。如果你以前从未遇到过 Synthea，看看我写的这篇解释其工作原理的简短文章。

Synthea 的输出数据分为几个 CSV 文件，如过敏、药物、遭遇、供应商等。

有些 CSV 文件的列标题不是列数据的明显指示器。查看 Synthea Wiki 有助于了解每个文件的内容。

样本输出数据格式

实际 CSV 文件的列与文档完全对应。

护理计划的 CSV 输出示例

对于这个医疗保健系统的例子，我生成了来自美国各地的 500 名患者的样本。

创建图表模式

完整的医疗保健模式

这个图表的模式非常复杂(正如您从图片中看到的)。我这样说并不是为了恐吓你，而是为了强调我已经说过的话。

医疗保健数据很复杂！

每个患者都与医疗保健系统有如此多的交互，一个简单的模式肯定无法捕获所有可用的数据和信息。我们的模式必须和我们的数据一样详细。

我们可以使用 TigerGraph 本地语言 GSQL 来创建我们的模式。整个脚本如下所示:

图形模式的 GSQL 脚本

光是这个脚本，我就可以写一整个博客。但是，为了使你免于阅读肯定会非常无聊的读物，我将只注意一些要点。

• 每个 CSV 文件主题都成为具有适当边的顶点。
• 所有的边都是无向的**,因为所有的关系都是双向的(即，一个病人有一种药物，但是该药物对应于该病人。)**
• 像性别**、种族和地址这样的属性可以是内部属性，但是我选择将它们分离出来以优化围绕这些属性的搜索。**
• 顶点 SnomedCode 存储使用的每一个医疗代码，这也有助于优化搜索

加载数据

我们可以编写 GSQL 加载脚本来加载我们的数据。让我们看看 CarePlans CSV 文件的加载脚本示例。

GSQL 加载脚本示例

同样，让我们简单地看一下这段代码的重要部分。

• 我们首先定义用来加载数据的文件。
• 然后，我们指定哪些列对应于模式中定义的顶点 id、顶点属性或边属性。
• 最后，我们声明我们的文件有一个头，分隔符是逗号

使用相同的格式，我们可以为其余的数据文件编写类似的加载作业。

对于我们 500 名患者的样本，加载我们所有的数据会产生大约 80 万个顶点和近****200 万条边**。**

示例查询

关于 GSQL 查询，我不会讲太多细节。我不想把重点放在实际编写查询上，而是放在查询运行的速度上——毕竟，这个博客旨在展示图形数据库的效率。如果你感兴趣的话，我有的另一个博客，里面有很多查询示例。和往常一样， TigerGraph 文档网站是寻找更多信息的好地方。

让我们运行一个简单的查询，获取与给定患者直接相关的所有顶点和边。

示例查询

该查询返回大量信息。它基本上调用了给定患者与医疗保健系统的每个接触点。通常，这对于数据库来说是一项艰巨的工作。但是有了我们的图形数据库，信息在几毫秒内就被检索出来了！真是快得不可思议！

这个速度也适用于大于 500 名患者的数据集。在一个拥有 1 亿患者(哇，这么多数据 ) 的样本系统中，收集相同信息的时间只有几秒钟**。**

使用图表进行查询非常高效，并且显示出比用于查询医疗保健数据的标准技术有了巨大的改进。

进一步探索

虽然图形数据库确实是保存数据的有效且高效的方式，但它们的好处不仅限于存储。我们已经看到了在大型数据集上执行查询有多快。但是，我们也可以利用图表进行可视化。例如，使用相同的数据库和稍微不同的查询，我们可以很容易地创建这个 3D 网络图，只需要几秒钟的渲染时间。

演示三维网络图

除了看起来很酷之外，这种 3D 视觉效果非常有用。虽然相同信息的 2D 表示会很混乱，无法阅读，但这个 3D 模型提供了一种开放而清晰的方法来查看我们大量相互关联的数据。此外，虽然美学部分是由 HTML 和 JavaScript 组成的，但数据是整个可视化的关键，它存在于图形数据库和查询中。

如果你想知道我是如何实现这个可视化的，看看我的另一篇博文 这里 ！

结论

与普通的关系数据库相比，图形是一个很好的选择，尤其是在表示高度关联的数据时。它们非常适合代表医疗保健网络，其中每个患者都与大量数据相关联。如果大规模实施，这项技术可以大大减轻 EHRs 的负担，并使数据的存储、分析和可视化更加有效。图表是医疗数据的未来！

我希望你喜欢这个博客，并学到一些新的东西。让我知道你的想法！

** [## 阿卡什·考尔-数据科学实习生-未来学家学院| LinkedIn

查看阿卡什·考尔在全球最大的职业社区 LinkedIn 上的个人资料。阿卡什有 5 个工作列在他们的…

www.linkedin.com](https://www.linkedin.com/in/akash-kaul-6a8063194/)**

几秒钟内建模:使用 PyCaret 作为数据科学快速决策的工具

原文：https://towardsdatascience.com/modeling-in-seconds-using-pycaret-as-a-tool-for-data-science-fast-decision-making-17c2c5642e17?source=collection_archive---------13-----------------------

弗兰基·查马基在 Unsplash 上的照片

我在浏览数据科学家的空闲时间时偶然发现了 Pycaret。这是一个多功能库，您可以在其中同时应用/评估/调整许多模型。

根据 PyCaret 文档:"py caret是 Python 中一个开源的低代码机器学习库，旨在减少从假设到洞察的周期时间。它使数据科学家和分析师能够高效地执行迭代的端到端数据科学实验，并且由于编码时间大大减少，他们可以更快地得出结论。”

这个库看起来非常类似于 R 中的 caret 库，但是是用 Python 实现的。

在从事数据科学项目时，通常我们需要花很长时间来理解数据(EDA 和功能工程)，那么如果我们可以将花费在项目建模部分的时间减少一半会怎么样呢？

我们已经有(许多)关于泰坦尼克号数据集的帖子，所以我将只展示该库的应用。不过，如果你很好奇想了解整个 EDA/数据操纵/特征工程，可以在 Kaggle 或者 GitHub 上查看我的内核。

1.建立图书馆

使用 PyCaret 的第一步是设置环境。这个设置的大部分是自动完成的，但是你可以手动设置一些参数。例如:

• 默认分流比是 70:30，但可以通过“列车大小”进行更改
• 默认情况下，k 倍交叉验证设置为 10
• “会话 id”是我们经典的“随机状态”

setup()函数的输出

setup()的输出是一个表格，其中包含您可以设置的所有可能的选项。如你所见，有很多选择。如果你想尝试不同的方法，这是非常有用的。例如，您可以使用其中一个特征选项来查看模型是否有所改进。

有趣的是，你可以用这个函数预处理你的数据。要查看所有选项，请查看库的文档。

2.比较模型

compare_models()的输出

compare_models()功能允许您一次比较多个型号。对我来说，这是使用 PyCaret 的一大优势。在短短的一行中，你就有了许多模型之间的表格比较。此外，它返回一些指标，如准确性、AUC 和 F1。另一件很酷的事情是库如何自动突出显示最好的结果。

一旦选择了模型，您就可以创建并调整它。在接下来的例子中，我将演示这些步骤。

3.创建和调整您的模型

让我们创建并调整我们的模型，看看库是如何处理它的。

使用 10 重交叉验证的 create_model()(左)和 tune_model()(右)的输出

我们创建了模型，并用 2 行代码对其进行了调整，如上所示。然而，我注意到的一件事是，模型的调整并没有提高我的模型的平均精度。我查了文档，上面写着:

"为了调整超参数，使用了 tune_model()函数。该功能在预定义的搜索空间中自动调整模型的超参数，并使用分层交叉验证对其进行评分。”

我在下面展示了为DecisionTreeClassifier调优模型的函数，因此我们可以对它的功能有一个概念。

from sklearn.tree import DecisionTreeClassifier param_grid = {"max_depth": np.random.randint(1, (len(X_train.columns)*.85),20), "max_features": np.random.randint(3, len(X_train.columns),20), "min_samples_leaf": [2,3,4,5,6],
              "criterion": ["gini", "entropy"],}

tune_model()功能将RandomizedSearchCV应用于考虑数据帧大小的模型。

4.模型性能可视化

PyCaret 的有趣之处在于它允许您轻松地评估您的模型性能。对我来说，创建这些图表有点乏味，所以有一个快速的方法来评估我的模型是完美的！有 15 种不同的地块可供选择。请参见下面的一些选项。

使用 PyCaret 可视化模型性能

与前面的步骤一样，生成这些图形的代码非常简单。我在下面的代码中展示了它:

5.混合模型

PyCaret 允许我们实现的另一种方法是混合模型。根据文献记载:blend_models()是一种集成方法，利用估计者之间的共识来生成最终预测。混合背后的想法是结合不同的机器学习算法，并在分类的情况下使用多数投票或平均预测概率来预测最终结果。

我们将尝试看看blend_models()是否会对我们的结果产生任何影响。blend_models()可以和一些预定义的模型一起使用，你可以用estimator_list通过。此外，在Classification的情况下，方法参数可用于定义“软”或“硬”，其中软使用预测的投票概率，硬使用预测的标签。为了简单起见，我就不通过车型列表了。

混合模型时的输出

我们的平均精确度从 0.8301 提高到 0.8316。没什么了不起的，但这是一点进步。

6.堆叠模型

PyCaret 还允许我们stack_models()。在文档中，它将堆叠模型定义为一种使用元学习的集成方法。叠加背后的想法是建立一个元模型，使用多个基本估计量的预测来生成最终预测。

与blend_models()不同的是，stack_models()要求通过estimator_list()，所以我选择了一些型号来尝试这种方法。

堆叠模型时的输出

如您所见，与其他方法相比，平均准确度没有提高。

7.总结

PyCaret 是一个多功能且易于使用的库。您可以用一行代码进行许多模型评估。

在比赛结束时，我使用 PyCaret 的最佳模型的平均准确度为 0.79 。然而，与任何新图书馆一样，仍有改进的余地。我将列出我在使用该库时发现的一些优点和缺点。

优点:

• 这使得项目的建模部分更加容易。您可以用一行代码创建许多不同分析。
• 在调优模型时，忘记传递参数列表。PyCaret 会自动为您完成。
• 您有许多不同的选项来评估模型，同样，只有一行代码。
• 由于它是建立在著名的机器学习库之上的，你可以很容易地将其与你的传统方法进行比较。

缺点:

• 这个库还是第一个版本，所以还不够成熟，容易出现 bug。
• 这是一个黑匣子，所以你真的看不到里面发生了什么。所以，我不会推荐新手使用。这会让学习过程变得有点肤浅。
• 在我的具体例子中，优化模型的函数并没有改进它。如果你们中的任何人对这件事有什么见解，请让我知道，这样我可以更新它。

暂时就这样了。我希望你喜欢这篇文章，如果你想查看我的一些作品，欢迎你访问我的 GitHub 或 LinkedIn 。

建模逻辑增长

原文：https://towardsdatascience.com/modeling-logistic-growth-1367dc971de2?source=collection_archive---------4-----------------------

冠状病毒的逻辑增长。疾控中心在 Unsplash 拍摄的照片

用 Python 中的 Scipy 对 Logistic 增长函数进行非线性最小二乘估计——以中国冠状病毒数据为例

在之前的一篇文章中，我解释了如何使用指数增长来模拟冠状病毒爆发的传播。它仅限于爆发的第一阶段，因为指数增长的最大限制是它永远不会停止增长。

在本文中，我将介绍这个模型的下一步:逻辑增长模型。如果你想跟着做，你可以在这里得到笔记本，在这里得到数据。

为什么是物流增长？

逻辑增长是一个数学函数，可用于多种情况。逻辑增长的特点是在开始阶段增长增加，但在后期，当你接近最大值时增长减少。例如，在冠状病毒的情况下，这个最大限度将是世界上的总人数，因为当每个人都生病时，增长必然会减少。

在逻辑增长的其他使用情况下，该数字可以是动物种群的大小，该动物种群呈指数增长，直到它们的环境不能为所有动物提供足够的食物，因此增长变得更慢，直到达到环境的最大容量。

使用逻辑增长来模拟冠状病毒爆发的原因是流行病学家已经研究了这些类型的爆发，并且众所周知，流行病的第一阶段遵循指数增长，并且总阶段可以用逻辑增长来模拟。

逻辑增长与指数增长

逻辑增长公式

逻辑增长的特征在于以下公式:

逻辑增长公式

其中:

• y(t) 是任意给定时间 t 的病例数
• c 是极限值，最大容量为 y
• b 必须大于 0

我还列出了这个公式的另外两个有趣的地方:

• 开始的情况数，也叫初始值是 : c / (1 + a)
• 最大增长率出现在 t = ln(a) / b 和 y(t) = c / 2

数学爱好者的逻辑增长公式的深度挖掘

对于那些想对公式有更好的数学理解的人:上面的函数实际上是从皮埃尔·弗朗索瓦·维赫斯特发现的下面的微分公式推导出来的。我会一步一步地解释:

• 绿色部分 dy / dt 表示这个公式不是针对人口规模的，而是针对人口增长的。
• 我们可以看到 y 和 c 在公式中，所以我们理解人口的增长取决于 y 的值(人口规模)和 c 的值(最大容量)
• 当 y 等于 c 时(即群体处于最大规模)，y / c 将为 1。因此，蓝色部分将为 0，因此增长将为 0。
• 当 y 比 c 小很多的时候(人口远离极限)蓝色部分将几乎是 1。因此，增长由橙色部分定义。橙色部分实际上是指数增长的公式。

逻辑增长的一个简单例子

为了更清楚地说明这一点，我将举一个假设的例子:

• 最大患病人数 c ，为 1000 人
• 我们从 1 个感染者的初始值开始，所以 c / (1 + a) = 1，给 1000 / (1 + a) = 1 ，给 a = 999
• 在感染开始时，每个病人会感染另外两个人，所以增长率b = 2
• 我们将考察从时间 0 到时间 10 的疫情发展

我们首先需要将 a 和 b 的值代入公式，以获得特定流行病的公式:

在 Python 中，这看起来如下:

然后我们可以用这个公式来计算从 0 到 10 的每个 t 值的 y 值。当我们这样做时，我们在每个时间点获得了以下感染人数，如下表所示。这表明，在开始时，感染人数增长很快，但随后就慢了很多，并在最大容量时结束。

如果我们想用图表来表示，我们开始看到一个图表，看起来很像我们看到的关于冠状病毒的非常惊人的曲线:

逻辑增长图

从冠状病毒数据到逻辑增长公式

现在，我们知道这个图或多或少有正确的形状，但我们需要做一个额外的步骤，使我们的分析有用。我们需要通过查看疫情传播的数据来找到电晕疫情的真实曲线。

在这种情况下，我将采用中国的数据:由于那里的增长已经大幅下降，逻辑曲线将是一个相当好的匹配。

中国冠状病毒数据摘录。来源:https://covid.ourworldindata.org/data/full_data.csv

中国物流增长的非线性最小二乘法

在查看数据时，我们只有每天的病例数。我们也有想要应用的公式，但是我们还没有公式中参数 a、b 和 c 的正确值。

不幸的是，不可能像指数模型那样将逻辑函数改写成线性回归。因此，我们需要一个更复杂的方法:非线性最小二乘估计。

第一步 —读入数据

第二步 —定义需要安装的物流功能

第三步 —随机初始化 a、b、c

在这种情况下，我使用 np.random.exponential，但是您可以使用您喜欢的任何东西。但是请注意:您在这里的选择可能会影响步骤 5 的性能。

第四步 —设置 a、b 和 c 的上下限

所有参数的下限都是 0。我把 b 的上界设为 3，因为在第一次尝试时我把它释放了，它变得太高了。a 和 c 的上限不会对曲线拟合产生负面影响，所以我让它们的上限相对较高。

步骤 5 —使用 Scipy 曲线拟合进行非线性最小二乘估计

在这一步中，Scipy 进行非线性最小二乘优化，并给出使模型的最小二乘误差最小化的 a、b 和 c 的值。

第 6 步 —绘制拟合函数与实际数据

正如你在下面的图表中看到的，逻辑模型与实际的中国冠状病毒数据相差不远。

第 7 步 —结束

我们发现了一个逻辑函数，它与从中国观察到的冠状病毒数据非常接近。为了在现实生活中真正使用这些信息，有必要进行大量的模型验证，比较不同模型的准确性和其他性能指标，并密切关注未来趋势是否符合所选模型。

这将远远超出本文的范围，本文只是试图说明如何拟合逻辑增长曲线，但是，使用该理论，我们可以陈述一些观察结果:

• 根据这个模型，c 是 81802，这将意味着中国感染人数的最大限度将是 81802。

我们还可以根据这个模型计算最大增长率出现的时间:

• 时间上的瞬间是:t =ln(a)/b = ln(61.5)/0.1929 =第 21 天
• 而那个时刻的感染人数是 y = c / 2 = 40500

有了一个完全有效的模型，这种类型的信息可以被决策者用来估计如何采取正确的措施。

我希望你已经清楚如何适应一个逻辑模型，以及如何在不同的用例中使用它。感谢您的阅读，敬请期待更多文章！

编者按: 走向数据科学 是一份以数据科学和机器学习研究为主的中型刊物。我们不是健康专家或流行病学家，本文的观点不应被解释为专业建议。想了解更多关于疫情冠状病毒的信息，可以点击 这里 。

用物理学模拟市场崩溃

原文：https://towardsdatascience.com/modeling-market-crashes-using-physics-485cebfddb8b?source=collection_archive---------12-----------------------

使用物理学概念识别市场崩溃的模式

格德·奥特曼德皮查拜的图片。

历史上有很多次市场崩盘(见这份名单)。重要的例子包括 1929 年的大崩盘，被称为“”(见下文):

图 1:1928 年至 1930 年间道琼斯工业平均指数的暴跌。

还有 1987 年的那个“黑色星期一”(见下):

图 2:1987 年 6 月 19 日至 1988 年 1 月 19 日的道琼斯工业平均指数(来源)。

最近，市场在冠状病毒的担忧中暴跌，这是自 2008 年金融危机以来最严重的崩盘:

图 3:最近，道琼斯工业平均指数在冠状病毒的恐惧中暴跌。

撞车是罕见事件的例子。正如这里将显示的，主要是在保罗和巴斯纳格尔之后，崩溃比通常的高斯或李维分布所能解释的更有可能发生。见下面两个分布的同步图。

图 4:高斯分布和李维分布的比较(来源)。

考虑到崩溃，市场的行为可以分为质量上不同的阶段:

• 价格波动相关性很弱的正常交易阶段
• 一个到崩溃的助跑阶段，其特征是价格波动之间的强相关性。

这种行为的变化将与多年来物理学中已知的一种现象相比较:相变。在相变中，有几个物理量发散的临界点。这些点将与市场崩溃进行比较。

崩溃的属性

正如在以前的一篇文章中详细讨论的那样(见下面的链接)，股票收益的分布更类似于李维分布，而不是高斯分布。前者有厚尾，在预测崩溃和其他类型的罕见事件时更有效。直观地假设，一次碰撞采样了尾部的非常远的边缘(用极端波动来标识)。

图 5:布朗运动和李维飞行的比较(来源)。

按照 Paul 和 Baschnagel 的说法，我们将把股票价格过程 S ( t )模拟为截断的无漂移李维飞行(参见上文李维飞行和布朗运动之间的比较)。

[## 股票价格的统计特性

模拟股票市场的动态

towardsdatascience.com](/statistical-properties-of-stock-prices-15145be752a2)

回想一下上面的文章，截断的李维分布是通过截断通常的李维分布的尾部得到的，使用的截断参数如下

等式 1:截断的李维分布。

其中 N 是一个归一化常数。下图显示了截断的李维分布和标准的李维分布。

图 6:截断的李维分布与标准李维分布的比较(来源)。

根据 Stanley 和 Mantegna 的说法，分布的中心部分将由 p ( l 表示，其中 l 是价格波动，由指数为 α = 1.4 的 l < 6 σ 近似，其中 σ 人们可以按照本文中的算法生成 Lévy 增量:

等式 2:无漂移的李维飞行。增量是从截断的 Lévy 分布中采样的。截止参数为δL = 6σ。

如果分布是对称的，指数 α 不等于 1，那么飞行变成:

等式 3:对称分布的 Lévy 飞行，其中 u 具有[-π/2，π/2]之间的均匀分布，v 具有平均值为 1 的指数分布。

情商。3 是对称分布的 Lévy 飞行，其中:

• u 是在[-π/2，π/2]之间均匀分布的随机变量
• v 具有平均值为 1 的指数分布。

注意截止的存在使得 p( l )的所有矩都是有限的。此外， S ( t→ ∞ )接近布朗运动。

除了最大偏移

让我们先定义某随机过程 x ( t )的最大值函数 M ( t ):

等式 4:给定间隔内样本路径 x(t)的最大值。

为了确定 M ( t )的分布，我们需要定义一个停止时间 T. 如果假设当 x ( t )第一次到达某个数字*a**，*时过程停止，我们可以定义:

等式 5:进程 x(t)的停止时间。

可以证明(参见 Paul 和 Baschnagel )标准布朗运动的最大值 M ( t )或最大偏移的概率密度为:

等式 6:标准布朗运动的最大值 M ( t )的概率密度。

在评估金融市场的风险时，这一结果至关重要。

我们对方程中参数的选择。如图 1 和 2 所示，极端事件的高斯行为是在大约 100 步的李维过程之后获得的。换句话说，价格波动的非高斯特征很快就消失了。

下降 人们也可以使用下降的分布来衡量极端事件。下降是价格时间序列中从相对最大值到下一个相对最小值的百分比下降。请注意，水位下降方法避免了在极端偏移分布分析中出现的固定时间范围问题(因为它对极端情况之间出现的不同时间延迟进行平均)。

图 7:提款 D 的相对数量 N(D)的线性对数图(来源)。

现在，在图 8 中，观察在 DJI 发现的两种状态。对于 D < 15% 保罗和巴斯纳格尔很好地符合一个指数。然而，两个图的比较表明，与模拟数据相比，DJI 在大约 30 %水平上的下降频率 D 大约高 10 倍。模拟数据预测提款的可能性要高得多。

图 8:DJIA 下降次数 N(D)与下降幅度 D 的对数线性图。这项分析每天进行，范围从 1900 年到 1994 年(来源)。

制度变化和相变

根据这些数据，有人推测(如引言中所述),这些罕见事件不仅仅是“正常”统计中的异常值，相反，它们表明了两种状态之间的转换:一种状态具有“正常”动态，另一种状态具有巨大的波动。

考虑图 9 中的以下两个崩溃。它们显示了 DJIA 在 1929 年和 1987 年崩盘前后的每日收盘报价。我们注意到这两个价格序列在崩盘前都是逐步上升的。这些阶跃是周期越来越小的振荡，在接近实际碰撞时消失(换句话说，随着碰撞时刻的临近，阶跃变得越来越短)。这种行为似乎是崩溃的前兆。

图 9:DJIA 接近 1929 年和 1987 年崩盘时的每日收盘报价。我们注意到，在崩溃之前，这两个系列都在逐步增加。此外，当我们接近崩溃时，步骤变得更短(来源)。

此外，价格下跌的幅度表明，它们不仅仅是“正常”价格序列的罕见统计波动。

合作性和关键性

当崩盘即将到来时，大量交易者自发地决定卖出他们的股票(有一种“合作性”)，动摇了市场的流动性。需求不匹配导致价格下降，市场进入非均衡状态。似乎发生的情况是，交易员之间形成了相关性，从而引发了随后的崩盘。

在热力学和统计力学中，当系统接近所谓的系统临界点时，会发生极其相似的现象。在这篇文章中，作者建立了一个简单的模型来提供市场行为和热力学之间的联系。

Cont-Bouchaud 模型 考虑市场中平均交易量相同的 N ( t )个交易者。交易者 i 可以:

• 如果 ϕᵢ = 1 则买入
• 卖出 if ϕᵢ = -1
• 如果 ϕᵢ = 0，请等待

供求差额由下式给出:

等式 7:供给和需求之间的差异。

提出了两个假设:

• 价格变化被假定为与供需不平衡成正比(这是小δS 所服从的)。
• 比例常数衡量市场对供求差异的敏感度(所谓的市场深度)。

通常，在市场中，交易者会相互“互动”。这种相互作用以两种不同的方式发生:

• 直接接触
• 间接的，通过跟踪价格的变化。

由于这种互动导致两个交易者采取共同的策略(他们要么买入，要么卖出，要么保持不动)，他们之间就产生了所谓的“纽带”。债券产生的概率是:

等式 8:债券创造概率，这对于大市场来说是个小数字。

考虑到 N(t) >>1，这是一个很小数量。参数 b 衡量交易者与同事互动的意愿。下一步是连接随机生成聚类的两个交易者，每个人对他们的活动都有不同的看法。

其中 ϕ ( c )=+1 或-1 包含活跃仓位的交易者，而ϕ(c)*=*0 包含非活跃仓位的交易者。参见图 10。

图 10:导致集群形成的步骤(改编自来源)。

等式中的和。7 可以重写为在 N ( c )簇上的和:

等式 9:等式中的和。7 以簇的形式书写。s 系数是聚类的大小。

其中 s ( c )是每个集群的大小。现在，我们将借用逾渗理论的一些结果，

图 11:逾渗(来源)。

根据该方程，无限程键逾渗模型的团簇尺寸概率分布为

方程式 10:无限范围键渗透模型的团簇尺寸概率分布。

这对 b 大约 1 有效。值 b =1 是逾渗阈值，有限簇具有 b < 1。当 b =1 时，一个单独的星系团出现并渗透整个系统。这是一个临界现象的例子，其中强时间相关性产生不同的动态机制。当所有的交易者都联系在一起时，就会发生渗透(平均而言)。如果b1，更多的交易者进入集群，开始控制整体行为。如果这个全跨度集群的成员决定出售，新的动力就会导致崩盘。我们必须有b1 来规避这样的场景。但是，我们必须记住，b 必须保持接近 1。否则，截断的 Lévy 行为的保存将被破坏。似乎发生的是，内在的市场动态使 b 接近 1。市场的本质应该是这样的，它将系统推向 b ≈ 1，并在“正常”行为期间将其保持在临界区域。现在，物理学充满了显示这种行为的系统。许多模型在没有任何外部控制(如温度)的情况下逐渐向动态吸引子演化，当它们接近该吸引子时，它们表现出临界行为，表现出幂律和厚尾 pdf。这里的类比是精确的，因为是系统本身的动态将它带向临界点。这叫做自组织临界。

图 12:森林火灾模型是显示自组织临界性的系统的一个例子(来源)。

因此，我们研究的模型认为股票市场是一个自组织的临界系统！

我的 Github 和个人网站 www.marcotavora.me 都有一些关于金融和其他主题的有趣资料，比如数学、数据科学和物理。看看他们！

用神经网络建模非线性动态系统

原文：https://towardsdatascience.com/modeling-non-linear-dynamic-systems-with-neural-networks-f3761bc92649?source=collection_archive---------21-----------------------

1.介绍

在 数学和科学中，非线性系统是一个系统，其中输出的变化与输入的变化不成比例。工程师、生物学家、物理学家、数学家以及许多其他科学家都对非线性问题感兴趣，因为大多数系统本质上都是非线性的。描述变量随时间变化的非线性动力系统，与简单得多的线性系统相比，可能显得混乱、不可预测或违反直觉。

我们将以工程系统为例，从统计/最大似然的角度来评估系统模型。工程系统的设计由原始设备制造商根据物理学/工程学的基本原理完成。数据科学家只获得输入和输出的集合。系统的建模——正式的系统识别——没有原始设计、领域知识和 OEM 的帮助；因此是不可能的。

那么，有没有可能根据 I/O 信息对系统建模；对于数据科学的非线性和动态系统来说也是如此吗？答案是肯定的，但即使在这种情况下，也需要一些领域专业知识来理解数据特性及其背后的数学基础。

仅基于观测数据的模型构建，没有关于其结构的任何详细信息；其特征仅在于其输入和输出，被称为黑盒模型。模型结构的元素没有物理意义。

当模型的结构和参数都完全已知时——完整的物理知识是可用的——我们就有了一个白盒模型。白盒模型可以根据先验信息构建，而不需要任何观测值。这也被称为实际系统的数字孪生或模拟。

2.什么是 NARX——非线性自回归外生建模？

图 NARX 模型的串并联架构。@版权

ARX 是自回归外生(ARX)的非线性推广，是线性黑箱系统辨识的标准技术。NARX 模型可用于建模各种各样的非线性动态系统。它是自动回归的，因为它的输出是以前值的回归。如果预测的映射函数是非线性的，则该模型称为 NARX。非线性映射可以是神经网络、高斯过程或其他。

图二。非线性映射函数。@版权

非线性映射函数被写成已知基函数的线性组合的情况将识别问题简化为线性回归。这是 NARX 模型非常受欢迎并成功应用于许多工业问题的一个原因。

对于静态系统，离散时间步长的输出仅取决于同一时刻的输入。静态网络可用于静态非线性系统建模。在动态系统中，给定时刻的输出不仅取决于当前的输入，还取决于系统先前的行为。动态系统是有记忆的系统——它们会随着时间而变化。

基于神经网络的 NARX 的架构如图 1 和图 3 所示。输入和输出都用抽头延迟线延迟。如果在输出信号路径中使用抽头延迟线，则可以构建反馈架构，其中前馈网络的输入或一些输入由网络的延迟输出组成。由此产生的网络是一个循环网络。该架构如图 3 所示。基本 NARX 网络用于多步预测。

假设目标的实际过去值是不可用的，并且预测本身被反馈到网络。因为在训练期间，我们可以访问实际的目标值，所以我们可以使用这些值来代替我们过去的预测。这有助于系统根据实际值而不是预测值进行训练。这是通过使用 NARX 网络的串并联版本实现的，如图 1 所示。该架构有两个优点。

首先是前馈网络的输入更加精确。

第二个是所得到的网络具有纯前馈结构，因此；静态反向传播可用于训练。

图 3。NARX 模型的并行体系结构。观察预测被用作延迟输入的闭环系统。@版权

3.发电厂示例

到为了让理解更具体，我们来建模一个电厂系统。数据来自 120MW 发电厂，参考 De Moor B.L.R .(版。)、 DaISy:用于系统识别的数据库，URL:http://homes.esat.kuleuven.be/~smc/daisy/,一座 120 MW 发电厂(法国桑布尔河畔)的数据。

数据采样时间为 1228.8 秒，大小为 200 个样本。输入有

1.气流

2.涡轮阀打开

3.过热器喷雾流量

4.气体阻尼器

5.气流

观察结果是

1.汽压

2.主阀杆温度

3.再热蒸汽温度

图 4 显示了电厂过程的输入传感器数据。

图 4。电厂工艺系统的输入

使用两个样本延迟、10 个隐藏神经元、Sigmoid 激活和 levengerg 损失函数来训练 NARX 网络。过去的外源输入和实际的目标观测值被用来预测当前观测值。如果没有记录实际目标观测值的过去值，网络将使用过去的预测值。但是，不建议用于工业先进过程控制的用例，因为预测会迅速偏离预期的估计。

图 5。输出预测和与观测值的比较。结果来自开环串并联架构。

结论

系统识别问题对于数据科学家来说是一个非常复杂的难题。

NARX 模型是能够以合理的精度对复杂动态系统进行建模的方法之一。

为了解决这些问题，需要用数据科学来增强领域专业知识、过程工程技能。

我们看到了一个串并联架构发电厂的建模示例。

下一步怎么样

还有哪些方法可以用来模拟非线性动态系统？将在下一篇文章中捕获它。

参考

1. Luigi Piroddi 博士,“非线性模型识别”课程讲义
2. 肖肯斯 1 和伦纳特·荣，非线性系统辨识，arXiv:1902.00683v1 [cs .2019 年 2 月 2 日
3. https://www.youtube.com/channel/UCAhRo1cCl1r_dFMVYaxTh5Q

使用可扩展和分层的变分推理对商店价格建模

原文：https://towardsdatascience.com/modeling-store-prices-using-hierarchical-variational-inference-6000f7701cb5?source=collection_archive---------35-----------------------

在本文中，我将使用来自 Kaggle 的 Mercari 价格建议数据,通过在 PyMC3 中实现的自动微分变分推理来预测商店价格。这个价格模型是由 Susan Li 在的这篇文章中使用 PyStan 完成的。我的一些作品受到了她的代码的启发。

我将更多地解释数据集、我执行的转换、我选择的模型以及我从中获得的准确性。请注意，本文的目的是探索不同的技术来更好地建模价格，并思考一些更好地构建数据的方法。我不会实现所有这些技术，但是我会在整篇文章中讨论它们。

价格建模的用例

每天，电子商务网站上都会发布数百个新的房源。因为电子商务网站竞争激烈，所以上市价格有很大的差异，如果上市价格不在一般客户的承受能力范围内，就非常容易失去一个客户。影响产品价格可行性的因素有很多。这包括产品的类型、状况、品牌以及免费送货等服务的可用性。在这个作业中，我们将把这些物品的价格作为这些因素的函数来建模。

预测列表的价格对电子商务公司非常有帮助，原因如下:

• 客户的价格决定会影响 Mercari 获得的在线流量，这与他们的声誉和广告利润等相关。因此，挂牌价格必须是最优的，而不是由客户随机生成的
• 知情的价格预测有助于卖家获得大量的客户关注，从而增加快速销售产品的机会，而无需等待更长的时间。
• 价格预测工具还改善了销售者和顾客的产品体验，因为价格是合理的，并且销售者不必花费大量精力来计算价格和降低价格以吸引顾客。
• 因此，价格预测工具将对卖家和客户的体验产生重大影响，同时提高 Mercari 的盈利能力和声誉。

数据

这些数据是从 Mercari 的 Kaggle Challenge 获得的，用于建立一个建议模型，根据客户提供的信息预测列表的价格。原始数据集包含大约 140 万个数据点。(梅尔卡利，2018 年)

对于每个数据点，我们都获得了以下指标:

• 条件 ID:指定产品的条件(从 1 到 5 不等)。
• 名称:列表的名称
• 类别名称:产品的一般类别。大约有 1500 个不同的类别。
• 品牌名称:品牌的名称，如果没有品牌，则为 Nan。
• 价格:客户列出的价格。
• 航运:如果航运可用
• 项目说明:列表的说明

数据清理、整合和转换

为了更好地分析数据，我们执行一些操作来合并文本数据:

• 数据中大约有 1500 个不同的类别。对所有这些类别建模是不现实的。因此，我没有使用 sklearn 的标签编码器，这是许多在线模型建议的选择。相反，我手动探索了这些数据，并将其分为 8 大类:男性、女性、电子产品、运动、美容、儿童、家居和其他。使用实际类别中的关键字，我将每个类别归类为八大类别代码之一。

**# Loading and processing the data here**
df = pd.read_csv(‘train.tsv’, sep=’\t’)
df = df[pd.notnull(df[‘category_name’])]category_codes = []**# Analysing the words in each category and classifying them as either of the 8 numerical categories**for i in df.index:
    category = df['category_name'][i]

    if 'Men' in category:
        category_codes.append(0)
    elif 'Women' in category:
        category_codes.append(1)
    elif 'Electronics' in category:
        category_codes.append(2)
    elif 'Sports' in category:
        category_codes.append(3)
    elif 'Beauty' in category:
        category_codes.append(4)
    elif 'Kids' in category:
        category_codes.append(5)
    elif 'Home' in category:
        category_codes.append(6)
    else:
        category_codes.append(7)df['category_code'] = category_codes

注意，上面的方法不是最佳的，但它是一种简单有效的情感分析方法。人们可以更深入地对每个产品的 category_name 和 item_description 进行情感分析，这也将加强处理后数据的质量。例如，可以分析诸如“良好状态”或“工作正常”等词的存在，并给它们打分/评级。不同种类的单词的存在将告知每个产品的总得分，这可以成为我们数据的一个特征。

• 条件 id 和运输代码未更改，因为它们是可用于分析的分类变量。
• 我不会在此分析中使用商品名称，但也可以使用情感分析来提取相关信息，如产品的尺寸或功能，以强化其价格。
• 对于品牌，可以建立一个简单的分类器，将较大的品牌分类为“2”，较小的品牌分类为“1”，没有品牌分类为“0”。为了简化模型，我不会实现这个策略，但是我鼓励您尝试一下。
• 我们把价格转换成对数标度。这使得价格分布更加平滑，从而便于建模。这也很重要，因为价格有一个乘数。你可以在这里阅读更多相关信息。但本质上，价格是通过乘以一个因子而变化的。例如，如果可以发货，价格将乘以 1.1。在贝叶斯统计中，对这些乘法因子进行建模更加困难，因为它们在采样期间更有可能发散。**因此，对数标度将乘法因子转换为加法因子，我们可以使用线性回归对其建模。**也使得价格成为实数线上的连续变量。

更好地理解数据

让我们看一下价格和对数价格分布，直观地了解这种转变:

价格与实际价格的对数分布。来源:自己的工作

在上图中，我们可以看到对数价格的分布类似于更容易建模的正态分布。

按航运可用性划分的原木价格分布。

上图挺有意思的。免运费产品的平均价格实际上小于免运费产品的平均价格，但我们可以预期相反的情况，因为运费会增加价格成本。这有两个问题:

1. 这些分布并不完全是正态分布，因此，执行 t 检验或比较平均值并不是比较这些分布的最佳方式。
2. 这里的比较分析也可能被混淆。可能的情况是，我们的无运费产品数据包含一般较便宜的产品，或者包含在分销中占主导地位的特定类别的产品。例如，大多数免运费的产品通常更便宜(也可以很容易地运送)，而更重、更大和更贵的产品通常没有免运费。因此，运费为 2 美元的 10 美元产品会被宣传为 12 美元，但更贵的产品是没有运费的产品，它们的价格(虽然不包括运费)会比包含运费的更便宜的产品高得多。

以上是数据中**结构偏差的一个经典例子，这种偏差源于我们作为人类所做的选择。**作为一名数据科学家，理解并解释这些源于非统计选择的偏差非常重要。没有简单的方法来纠正这种偏见。一种方法是执行个体匹配。基本上，我们会浏览数据集，挑选看起来完全相同或非常相似的记录。他们之间唯一的主要区别是免费送货的可用性。因此，我们将构建一个包含两个子集的新数据集，两个子集的产品分布非常相似，但一个没有发货，而另一个有更多发货。

您可以在文章末尾的笔记本中找到这个可视化的代码。欢迎您绘制自己的可视化图，并思考数据集中可能存在的偏见。

车型

我们先来考虑一个简单的线性回归。(虽然线性回归看起来非常简单，但它们功能强大且可扩展。因此，数据科学家更喜欢用线性回归建模的方式处理他们的数据。)

这里，c 是模型中的常数，而 a 和 b 是运输系数和条件系数。为了使它变得现实，我将使用变量池。这意味着我将为每一类产品创造不同的 a、b 和 c 值。因此，“电子产品”将有自己的一套参数，“家居装饰”也将有自己的一套参数。如果我没有对一个或多个参数进行池化，这将被称为部分池化。让我们更新等式:

这里,[8]表示参数反映了 8 个不同的参数，每个参数对应一个指定的类别。

a、b、c 的含义:

我们用 a、b 和 c 建立了一个模型，但是它们是什么意思呢？就线性回归而言，c 是截距，而 b 是条件系数和运输系数。在价格方面，“c”就像是每个类别的基础价格。此类别中的每个产品都有一个基准价格。该基线受产品状况的影响。这种影响的规模由‘b’控制。原木价格还受到运输的影响。这种影响的规模由“a”控制

建模基础设施

既然我们已经定义了模型的基本结构，我们需要构建基于这些参数的基础先验分布。我将在这里使用层次模型，因为层次提供了值的灵活性，并允许同一参数的不同值有它们的分布，而不是依赖于同一个值。

例如，每个“a”将从平均值和标准偏差的正态分布中取样，这些正态分布将从它们自己的先验中进一步取样。因此，每个“a”都有自己的正态分布，这对于满足联营假设很重要。否则，它们可能会相互关联，从而威胁到模型的有效性。

这里‘a’的超上级是 mu 和 sigma。μ是从正态分布采样的，而标准差是从 0 到 10 之间的均匀分布采样的。请注意，标准差必须为正，因此，我们不能使用正态先验，因为正态分布的支持是整个实数行。这些先验被选为未知先验，因为我们对这些系数或它们的影响知之甚少。我们将让模型找到并收敛于最佳可能值。

我们也为“b”和“c”复制相同的模型。我们将正态分布置于两者之上，并通过其超参数“mu”和“sigma”进一步参数化。最后，似然性被建模为另一个正态分布，该正态分布被置于具有噪声“sigma”的最终线性方程上。

执行推理

现在我们已经建立了我们的模型，我们想进行推理，以获得 a、b 和 c 的后验分布，它通知我们理解模型的能力。我们可以允许 MCMC 从中生成样本，或者我们可以执行变分推断，将这些分布放在我们的参数上，并优化它们的超参数值。让我们从这个分析中回顾一下我们的目标:对产品价格进行可伸缩的推断。

对 MCMC 及其缺点的思考；

Susan Li 在的这篇文章中使用 PyStan 进行了类似的分析。PyStan 是一个很棒的工具，它使用一个 NUTS 采样器来执行自动哈密尔顿蒙特卡罗采样。但是 HMC 是一个非常慢的算法。在这个数据集中，我们有大约 140 万个数据点。这使得算法的平均运行时间约为 2-3 小时。虽然有不同的方法来优化 PyStan 的效率，但它不是一个可伸缩的算法。你可以在这里了解更多关于 PyStan 的线性回归。

PyStan 的另一个重要限制是缺乏可伸缩的适应性。随着新数据的生成和处理，我们希望我们的算法能够实时更新它们的参数，以便模型每天都变得更加智能。我们也希望这一过程快速高效地进行。因此，PyStan 可能不是在线市场模式的最佳选择。你可以在这里阅读更多关于 PyStan 的信息。

自动微分变分推理

如果你不太了解 MCMC 或者变分推理(VI)，你应该读一下这篇文章。总而言之，VI 对给定的参数进行分布，然后优化这些分布。

VI 使用了一种叫做 KL-Divergence 的东西，它是一种成本函数，可以帮助我们理解我们的结果应该是什么样子，以及我们想要达到的目标之间的距离。但是 VI 不是全自动的，需要微调。因此，我们将使用一个称为自动微分变分推理(ADVI)的 VI 版本，它优化了一个不同的参数(类似于 KL-Divergence，但更容易优化)。该算法使用不同的库来区分 ELBO(我们这里的度量)并找到最大值来更新参数。它使用坐标下降(沿坐标轴优化)来实现参数更新方程。我知道这是很多信息，但是我推荐你阅读这篇文章来熟悉它。

**这里有一个重要的注意事项:**当我们在 VI 中执行推理时，我们对定义每个模型参数(a、b 和 c)的分布的变分参数(这里是超参数)执行推理。因此，模型变量成为一个过渡参数，其分布通过对变化参数(所有 mu 和所有 sigma)的推断来确定

在 PyMC3 中建立模型

对于这个实现，我们不会从头开始构建 ADVI 算法。相反，我们将使用 PyMC3 的内置特性来构建我们的模型，并指定超先验、先验、模型和推理策略。

为了使推断更快，我们将在 PyMC3 中使用迷你批处理。迷你批次减少了训练数据集的大小。在初始设置中，参数更新将要求参数查看整个训练数据集。不幸的是，我们有 140 万个数据点要看。对于大约 25000 次迭代(我们的目标)，这是一个很大的工作量。因此，迷你批次允许我们从训练数据集中取样。他们创建训练数据集的小批量(大小为“b”)，并且每次更新的训练都在小批量上进行。每次迭代后，批次都会更新。

经过大量的迭代，微型批次在优化 ELBO 方面变得非常有效，因为更新方程几乎可以查看每个数据点，但不是同时查看所有数据点，这使得该过程更具可扩展性。

下面是为我们构建模型的代码。(我的灵感来自于 PyMC3 网站上提供的这个分层部分池示例，其中一些代码是从他们的笔记本中获得的。)

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np
%matplotlib inline
%env THEANO_FLAGS=device=cpu, floatX=float32, warn_float64=ignore
import theano
import theano.tensor as tt
import pymc3 as pmcategory_code_new = df.category_code.values
shipping = df.shipping.values
condition_id = df.item_condition_id.values
log_price = df.log_price.values
n_categories = len(df.category_code.unique())**# Minibatches below. The batch size is 500**log_price_t = pm.Minibatch(log_price, 500)
shipping_t = pm.Minibatch(shipping, 500)
condition_t = pm.Minibatch(condition_id, 500)
category_code_t = pm.Minibatch(category_code_new, 500)**# Setting the hyperpriors (mu,sigma) for a,b,c**with pm.Model() as hierarchical_model:
    mu_a = pm.Normal('mu_alpha', mu=0., sigma=10**2)
    sigma_a = pm.Uniform('sigma_alpha', lower=0, upper=10) mu_b = pm.Normal('mu_beta', mu=0., sigma=10**2)
    sigma_b = pm.Uniform('sigma_beta', lower=0, upper=10) mu_c = pm.Normal('mu_c', mu=0., sigma=10**2)
    sigma_c = pm.Uniform('sigma_c', lower=0, upper=10)**# Setting the priors here. Note the shape variable.
# The shape must be equal the number of categories. (8 here)**with hierarchical_model:
    a = pm.Normal('alpha',mu=mu_a,sigma=sigma_a, shape=n_categories)
    b = pm.Normal('beta',mu=mu_b, sigma=sigma_b, shape=n_categories)
    c = pm.Normal('c', mu=mu_c, sigma=sigma_c, shape=n_categories)**# Setting the likelihood here as another gaussian**with hierarchical_model:
    **# Estimated log prices** log_price_est = c[category_code_t] + b[category_code_t] *
                     condition_t + a[category_code_t] * shipping_t **# Noise or standard deviation of the likelihood**
    noise = pm.Uniform('eps', lower=0, upper=100) **#Likelihood function. Note that we still need to provide the 
    # total size of the dataset for convergence check.** lp_like = pm.Normal('lp_like', mu=log_price_est, sigma=eps,
             observed=log_price_t, total_size=len(df))

优化时间！

让我们在 PyMC3 中将优化策略定义为 ADVI，并执行推理。

from pm.callbacks import CheckParametersConvergence**# Defining inference method and parameter checks**with hierarchical_model:
    inference = pm.ADVI()
    approx = pm.fit(100000, method=inference, callbacks= 
                    [CheckParametersConvergence()])plt.figure(figsize=(12,6))
plt.plot(-inference.hist, label='new ADVI', alpha=.8)
plt.plot(approx.hist, label='old ADVI', alpha=.8)
plt.legend(fontsize=15)
plt.ylabel('ELBO', fontsize=15)
plt.xlabel('iteration', fontsize=15);

ELBO 在多次迭代中的收敛性。旧 ADVI 和新 ADVI 提供了对每次迭代级别的优化的深入了解。来源:自己的工作

从图中，我们可以说 ELBO 收敛得很好，并在大约 25000 次迭代时达到最小值。我们观察到最小值出现在 490.22，而最大值出现在 4706.71。你可能会问——在 10000 次迭代左右，ELBO 中的峰值是什么？这些更新违背了我们的目标。发生这种情况可能有多种原因，最有可能的是微型批次中的巨大偏差，导致参数更新偏离其轨迹。但正是随机化的美妙之处(和数学原理)帮助我们最终得到正确的结果。

变分参数的后验分布

还记得我们没有执行 MCMC 吗？是的。所以我们得到的不是样本，而是变分参数的值，它们是模型参数的超先验，然后也是模型参数的超先验。我们通过拟合模型得到的变量“近似值”将为我们提供每个参数的平均值和标准偏差。我们将获得一个均值和标准差字典，可以用参数值进行索引:

means = approx.bij.rmap(approx.mean.eval())
sds = approx.bij.rmap(approx.std.eval())a_mean = means['alpha']
a_std = sds['alpha'] **# for each model and variational parameter**

我们对变分参数的分布不太感兴趣，因为它们的先验不是动态的。我们将快速浏览一下它们的均值和标准差:

1. mu_a ~ N(-0.00022505，0.7320357)
2. mu_b ~ N(6.097043e-05，0.7320404)
3. mu_c ~ N(0.00212902，0.7320338)

正如所料，它们以 0 为中心，因为这些管理部门的先验也以 0 为中心。标准差已经重新调整了(从最初的 10 到小得多的 0.732)。所有这些意味着分布彼此非常接近。因此，a、b 和 c 的实际值是从彼此非常接近的高斯样本中采样的。需要注意的重要一点是，虽然均值的后验概率是相似的，但基于所执行的推断，模型参数的实际后验概率可能会非常不同。

模型参数的后验分布

我们使用下面的代码来绘制后验分布的 PDF。(这里不能做直方图，因为我们没有一组样本。)

from scipy import stats**# List of categories**
category_names = ['Men','Women','Electronics','Sports','Beauty', 
                  'Kids', 'Home', 'Other']**# Model parameters names**
model_param = ['alpha','beta','c']fig, axs = plt.subplots(nrows=3, figsize=(12, 10))**# Plotting a set of PDFs for each Model parameter**
for var, ax in zip(model_param, axs):
    mu_arr = means[var]
    sigma_arr = sds[var]     ax.set_title(var) **# Title of the subplot** for i, (mu, sigma) in enumerate(zip(mu_arr.flatten(),
                                 sigma_arr.flatten())):
        sd3 = (-4*sigma + mu, 4*sigma + mu)
        x = np.linspace(sd3[0], sd3[1], 300)
        y = stats.norm(mu, sigma).pdf(x) **# PDF of the parameter**
        ax.plot(x, y, label='Coefficient of {}'. format(
                                           category_names[i]))
    ax.legend()
fig.tight_layout()

以下是混合模型参数的后验分布:

让我们分析模型参数:

1. 对于基准价格©，男装的基准最高，其次是女装、家居装饰和运动。其他产品有最便宜的基线。这可能是因为较小的杂项产品没有很多价值，因为它们不属于任何主要类别。
2. 对于 beta (b ),产品条件对某些类别的影响是积极的，而对其他类别则是消极的。更重要的是，它是对数线性模型中的一个系数。因此，这些系数代表实际价格的变化百分比。恶劣条件的负面影响最大的是女装，恶劣条件会导致价格下降 15%。这同样适用于男装和家居装饰。运动和童装下降了 5%。有一个电子价格的增长(约 12%)与一个坏的条件，这是违反直觉，但同样，这将需要验证匹配。
3. 对于阿尔法(a)，免费送货的影响很大程度上是负面的(正如我们在之前的初始分析中看到的)。男装的负面影响最小(约-18%)，而电子产品的负面影响最大(约-50%)。你应该更多地考虑数据的结构偏差，以及类别是否会加强偏差。**提示:**大多数电子产品需要更多的运输/处理费用，因此不提供免费运输。

从上述分析中，我们可以得出结论，最违反直觉的结果出现在电子产品和其他产品类别中，而其他产品没有出现极端结果。为了更好、更准确地建模，对该数据集执行匹配是值得的。

不确定因素

从后验分布可以看出，很多分布是非常确定的(比如女装的参数)。另一方面，运动产品的参数更加不确定(标准偏差为 0.025，约为实际价格的 2.5%变化)。因此，当我们根据这些结果建立预测模型时，我们需要小心这些不确定性，并将它们纳入我们的结果。

一个简单的预测模型

从技术上来说，我们可以从我们的线性回归结果中建立一个确定性预测模型:

means = approx.bij.rmap(approx.mean.eval())def predict(shipping, condition, category):
    c = means['c'][category]
    b = means['beta'][category]
    a = means['alpha'][category] return np.exp(c + b*condition + a*shipping)

但是 VI 没有提供预测区间，而 MCMC 提供了。要获得预测区间，我们只需构建 a、b 和 c 的置信区间。然后，我们使用这些区间来评估相同输入在 a、b 和 c 的下限/上限处的相同预测。

与 HMC 坚果取样器结果的比较

在同一个模型上，我们也可以执行 HMC 螺母抽样来获得后验分布。下面是实现这一点的代码:

with hierarchical_model: 
    step = pm.NUTS(scaling=approx.cov.eval(), is_cov=True)
    hierarchical_trace = pm.sample(100000, step=step, 
                         start=approx.sample()[0], progressbar=True, 
                                    tune=1000, cores=4)pm.summary(hierarchical_trace)

请注意，PyMC3 的 NUTS 采样器(主要是 summary 函数)有一些错误，因此，我没有广泛地使用它，但您仍然可以使用 MCMC 获得一些关于我们结果的收敛性的很好的见解。

**剧透警告:**锁链不会收敛。原因可能在于使用了微型批次，这导致链分叉，因为 MCMC 并不完全执行更新参数的迭代。它可以创建许多您想要生成的样本。因此，如果样本是基于不同的小批次生成的，这组样本可能会变得非常不同(就是这种情况)。在这种情况下，您将需要更多的迭代和链来完成收敛。

未来的步骤和改进

有许多方法可以改进目前的模式。我们讨论了一些更好地分析文本数据的方法。以下是 Mercari 在数据收集方面可以做得更好的一些事情。Mercari 还可以在收集更多准确且易于处理的数据方面做得更好。我们可以设置 5 个问题来帮助用户更好地对产品进行分类，而不是要求键入描述。这也将使我们能够收集易于分析的良好数据。例如，我们可以问这样的问题:“产品有多长时间了？”，‘产品有用吗？’、‘是翻新的还是有保修的？’，‘你买的时候花了多少钱？’。等等。

2.在构建价格预测器时，还应考虑客户的偏好。每个顾客在等待时间、预期价格、最初支付的价格、情感/社会价值(属于名人)等方面都有一定的偏好。因此，价格应该根据客户的偏好和情况进行调整。

测试笔记本

在这里随意克隆 Python 笔记本并测试不同的策略。

*注意:*如果你有任何问题或者你想联系，我的 LinkedIn 个人资料是这里的。

模拟摆动概率

原文：https://towardsdatascience.com/modeling-swing-probability-b02b5ab7fbb5?source=collection_archive---------33-----------------------

戴戟·梅本在 Unsplash 上的照片

多年来，棒球的一个方面让我着迷，那就是击球手和投手在击球时进行的心理博弈。每个玩家都在不断地试图进入对手的大脑，猜测他们下一步可能会做什么。击球手可能会利用他对投手的了解来预测他是否会尝试用快速球挑战击球手，或是引诱他将突破球追出区域。与此同时，投手正在利用关于击球手的信息来制定一系列投球，这些投球应该会让他因三振出局而回到替补席。

正是在这些智力游戏中，我发现统计学和 sabermetrics 可以得到最有效的应用。投手或击球手拥有的相关数据越多，他们的优势就越大。在我之前的一篇文章中，我量化了投手隐藏投球的程度，并讨论了击球手如何利用这些信息来识别投手可能会投出什么(不需要垃圾桶，太空人)。然而，在这篇文章中，我想调查投手如何利用打者的信息来给自己一个优势。特别是，我将试图测量某个击球手在给定投球时挥棒的概率。

型号:

为了构建这个模型，我将使用 Python 和它的机器学习库 LightGBM，我将把它用于它的梯度增强分类工具。如果你不知道什么是梯度增强分类算法，并想知道更多，请随时查看 StatQuest 的这个优秀视频进行深入了解。现在，关于我们的模型，你需要知道的是，它本质上是一个二元分类算法，其中两个分类要么是挥杆，要么是拿杆(注意:我们不关心挥杆的结果，只关心击球手是否挥杆)。

数据:

与大多数机器学习算法一样，如果我们希望我们的模型准确有效，我们需要大量的数据。幸运的是，有一个名为 pybaseball 的方便的 Python 库，使得获取这些数据变得轻而易举。我将使用这个库来收集 StatCast 逐节数据，以便训练模型。

特性:

特征选择是构建分类模型不可或缺的一部分。这些特征是算法将考虑的不同度量，以确定击球手是否会以特定的音高挥杆。因此，我们必须纳入对结果变量影响最大的指标。StatCast 跟踪大量不同的投球指标，如旋转速度、释放点等，但只有其中一些指标对决定击球手是否挥杆有用。我将在模型中使用的最重要的功能是:

1. 球场的位置
2. 投球的速度
3. 音高的垂直和水平变化
4. 投球的计数
5. 前一个投球的投球类型(除非它是击球手的第一个投球)

这些特征中的一些可能是显而易见的，但是有一些可能不是很直观。位置绝对是决定击球手挥杆与否的最重要因素；一个打者几乎不会在 3 英尺以上的区域挥棒，但经常会在中间挥出快速球。

速度和垂直/水平突破功能一举两得。首先，它们在决定击球手是否挥棒时非常重要(击球手经常挥棒打一个看起来像好球的滑球，但最终离开本垒板 6 英寸，快速球经常会冻结击球手期望的变速球)。其次，它们允许我隐式地告诉模型投掷的是什么类型的球。我不必明确地告诉模型具体的音高类型，而是算法可以根据音高的速度和移动来学习音高的行为。这在我们以后用模型做预测的时候会很有帮助。

投球的次数也是击球手决定挥杆与否的一个因素。例如，击球手通常对击球手的计数更有选择性，比如 3-0，除非完美，否则不会挥棒击球；但是当打者以 0-2 落后时，他通常会挥棒攻击任何接近好球带的地方。

最后，前一个投球的投球类型在决定击球手是否会挥棒方面起着更微妙的作用。一个刚看到曲球的击球手可能不太可能挥出随后的曲球，因为他刚看到一个曲球。

现在我们有了所有这些数据，我们终于可以开始训练模型了。

训练模型:

为了训练这个模型，我选择了红队的一垒手乔伊·沃托作为我们的击球手。我选择沃托是因为他在挥杆时极其挑剔，我想看看我们的模型是否会反映这一点。在用乔伊·沃托过去 5 个赛季的逐节数据训练模型后，我们得到了 0.795 的准确率。这意味着我们的模型能够在 79.5%的时间里正确预测沃托是否在投球。虽然并不完美，但这似乎很有希望。期望接近 100%的准确性是不现实的，因为人为误差和棒球的随机特性使得不可能绝对确定地预测。经过一点点的超参数调整，我能够获得高达 80.5%的精度分数，这对于我们的模型来说已经足够好了。

使用模型:

既然我们的模型已经训练好了，我们终于可以开始用它来做预测了。虽然知道一个特定的音高是否会被挥杆是有帮助的，但输入特定的音高来看它们是否会产生挥杆是非常低效的。即使我们这样做了，我们的模型也只有 80%的准确性，所以会有相当大的误差。相反，我们可以绕过这两个问题，在好球区周围绘制一个特定的球距图，并使用我们的模型来计算沃托在该区域每个位置的挥杆概率。

*从捕手的角度看好球带

这个图表给了我们模型的预测概率，如果乔伊·沃托在击球手的第一次投球时投出泰勒·格拉斯诺的快速球，他将挥棒击球。以这种方式使用模型的预测概率更有帮助，提供了更多的洞察力，而不是让模型简单地指示投球是否会被挥杆。该模型将预测挥杆概率为 51%的投球和挥杆概率为 98%的投球，因为两者都被挥杆，但它们显然是两种非常不同的投球。通过可视化数据，我们可以看到 Votto 以某一音高摆动的可能性，而不是简单地得到是或否的答案。

看前面的图表，我们可以看到，在球板中心上方的球路上，沃托挥出快速球的概率约为 80%,但当你到达好球带的边缘时，这种概率会迅速降低。这似乎与我们对乔伊·沃托的了解一致。他是一个非常挑剔的击球手，有着精英的眼光，所以当你远离本垒板的中心时，他的挥杆概率就会下降。这一点在击球的第一次投球时得到了进一步的强调，在这种情况下，没有必要挥动任何不容易击中的东西。基于第一个结果，我们的模型看起来是准确的，但是让我们更深入一点来确定。我们可以测试模型准确性的一种方法是通过观察 Votto 在 0-0 计数中的摆动率，看它是否与我们的模型给出的相似。

看这个我们可以看到我们的模型的情节看起来惊人的相似。这是一个好迹象，表明我们的模型是准确的。看起来确实有些细微的差异，但那很可能是因为挥杆率图表考虑了所有的球种，而不只是快速球。为了进一步测试这个模型，让我们来看看沃托在另一场比赛中的挥杆概率。

该图显示了在 0-2 计数中，沃托的挥杆概率与克莱顿·克肖的曲球。这个模型似乎又一次符合我们的预期。在 0-2 的计数中，Votto 最有可能在任何看起来可能是好球的投球上挥杆，这是由他在低于区域的投球上相对较高的挥杆概率所表明的。此外，由于落在区域顶端的曲球很少被挥击，所以当你在好球带中走得更高时，沃托的挥击概率会降低是有道理的。

应用模型:

既然我们已经看到我们的模型是相对准确的，我们可以考虑如何在现实世界中应用它

第一个也是最明显的应用是增强投手的打者侦察报告。虽然挥杆率数据是可用的，并且可以在 Baseball Savant 等网站上很好地可视化，但这些数据通常非常有限或过于笼统。如果一个打者在他的职业生涯中只面对过几次投手，那么样本量太小，不足以揭示任何有意义的见解。类似地，如果你试图对一种特定的投球类型进行归纳和建模，你就失去了特异性。这样的模型将克莱顿·克肖的缓慢循环曲球与泰勒·格拉斯诺的快速突破曲球归为一类，这损害了模型在应用于个案时的准确性。通过使用我们在这里制作的模型，您可以获得两个世界的最佳效果。即使一个打者没有面对过某个投手，他也很可能面对过几个不同的投手，他们有相似的投球。通过使用单个投手的平均投球动作/速度，该模型可以相对准确地预测击球手的行为，即使关于他们具体比赛的数据很少。这有助于制作适合特定投手的球探报告，并提供对击球手行为的洞察，否则这些行为会被忽略。

其次，这个模型也可以用于模拟，试图预测击球手或投手在比赛中的表现。使用像这样的模型来预测挥杆概率，同时使用其他算法来预测挥杆结果，可能有助于预测特定投手的击球结果。

如果你做到了这一步，感谢你的阅读，并随时给我留下任何提示和/或批评！如果你想查看这个模型的代码，这里有一个到 GitHub 库的链接。

模拟纽约市的交通事故风险

原文：https://towardsdatascience.com/modeling-the-risk-of-traffic-accidents-in-new-york-city-9292c5ae4ab4?source=collection_archive---------28-----------------------

使用 convLSTM2D 深度学习模型对交通事故风险进行分类

来源:帕特米克凯恩，经由皮克斯拜 (CC0)

问题陈述

仅在 2017 年，美国就有超过 3.7 万人死于道路交通事故。2010 年与交通事故相关的成本(来自 2015 年发表的一份报告)为 2420 亿美元，占当年国内生产总值的 1.6%。以运输为生的公司将受益于关于哪里更可能发生事故的信息，因此他们可以避开这些地区，降低事故风险。优步和 Lyft 等拼车公司依赖于快速安全地从一个地方到达另一个地方，以获得客户满意度。优步为其司机提供保险，承担事故成本，因此如果发生事故，面临的风险不仅仅是客户满意度下降。这个项目将集中于整个纽约市事故风险的时空预测。拼车公司在生成路线以避开事故风险增加的区域时，可以使用这样的信息作为额外的数据层。

我们将使用的数据是由纽约市提供的 NYPD 机动车碰撞数据。

• https://data . cityofnewyork . us/Public-Safety/NYPD-机动车-碰撞/h9gi-nx95

数据以 CSV 格式提供，包括事故时间、地理坐标和其他信息，如事故原因、涉及的车辆、死亡人数等。数据将被清理，以确保每个条目包含有效的坐标和时间戳。在这一阶段，我们计划将数据表示为栅格，以确定一天中特定时间事故风险增加的区域。随着时间的推移而可视化的 2D 矩阵基本上是一个图像识别问题，因此我们计划开发一个卷积 LSTM 深度学习模型。像这样的方法是有先例的。

数据争论

我们决定将此作为一个分类问题，提出这样一个问题:在这一天的这个时候，这个地区发生事故的风险是否会增加。数据集包含 1.58 MM 的事故行，包括日期、时间、经度和纬度，以及涉及的车辆数量、促成因素、受伤人数等。在之前的项目中，我们对纽约市范围内发生的所有事故都很感兴趣，因此不关心是否有缺失的位置值。对于本项目，事故的位置和时间非常重要，因为我们试图建立一个模型来预测在选定的时间和日期发生事故的风险增加的位置。

通过上面的链接下载的. csv 文件中的数据被导入到熊猫数据框架中。列标题被转换为小写，标题中的空格被替换为下划线以允许链接。创建了一个新的 datetime 列，并将其设置为索引，以便使用 pandas 进行索引。loc 访问器。数据集中的每条线代表向 NYPD 报告的事故的时间和位置。最初检查数据时，日期时间数据的行数与经度和纬度的行数不匹配。更仔细的检查揭示了超过 185，000 个 NaN 实例，以及超过 900 个经度和纬度的 0 值。0 经纬度在赤道上，加纳以南。由于事故的日期和位置对这项工作至关重要，所有的 NaN 值首先被转换为 0，并与现有的 0 值一起删除。这些条目的删除导致了与经度和纬度相同数量的日期时间值。

所有事故都是用它们的纬度和经度值绘制的。这显示了多个外围点，很明显不在纽约市范围内。对经度和纬度的限制越来越严格，直到通过目视检查获得了一幅完整的纽约市地图。接下来，我们将经纬度转换为世界大地测量系统，这是将地球曲面投影到平面上的标准，就像地图一样。这一点很重要，因为地面距离测量会更精确。

并列显示坐标(绿色)和 WGS 投影(蓝色)

这种差异很小，但在地面测量中会造成几十米的误差。美国东海岸位于通用横轴墨卡托坐标系的 18T 区域。更多信息见此处。应用此功能可将纬度和经度转换为米东距和米北距，这是一种更精确的地面距离测量方法。这些值将用于这项工作的剩余部分。

电子设计自动化(Electronic Design Automation)

为了感受数据的一致性，我们首先按小时对事故进行分组，并绘制了两周的直方图。这表明在一天中似乎有一个双峰分布，大约在上午 8 点和下午 4 点达到峰值，这与正常的高峰时间交通一致。

从周一开始，绘制了两周的直方图，显示每小时的事故数量。

双峰分布对许多人来说是明显的，并且总是在下午发生最多的事故。还有一段时间是凌晨，事故比较少，大概是因为路上车少。然后，我们绘制了整个数据集内每月的事故总数。

数据集期间每月的交通事故总数。

事故数量似乎是季节性的，年初时最低，夏季达到高峰，接近年底时逐渐减少。然而，我们不只是对预测事故何时发生感兴趣。我们感兴趣的是事故将在何时何地发生。为此，我们将按小时和地点汇总事故。如前所述，我们将把它放在 2D 数组中。最初，我们选择将数据离散成 64 x 64 的网格。这返回了一个更加像素化的数据视图，但是城市的许多特征仍然是可辨别的。

所有的事故数据都被收集到一个 64 x 64 的数组中

像这样分组为 64×64 个箱使得每个箱的宽度大约为 712×765 米。虽然这不是绕过事故的最实际距离，但它将允许我们展示我们的模型预测交通事故的能力。为了在合理的时间内进行计算，我们决定查看最近 6 个月的事故数据，该数据集包括 2018 年 12 月 16 日至 2019 年 4 月 16 日的数据。数据被分成 70:20:10 的训练:测试:验证集。使用训练集和测试集改进的模型，验证集是一个保留集，仅应用于训练模型。

模型结构

我们选择用 2D 事故阵列来表示这个模型。这两个维度由规则的距离间隔定义，每个箱中的事故数量相加。这给了我们一个纽约市事故的像素化视图，但它适合看看这种方法是否有效。像这样排列的数据可以被认为是一幅图片，或者当考虑多个时间片时，是一部电影的帧。有鉴于此，我们可以把它当作一个“预测电影的下一帧”的问题。这是一个有趣的问题，因为这是一个时间序列问题，也是一个图像问题。时间序列问题通常受益于长短期记忆(LSTM)模型。这些是递归神经网络，具有存储先前信息以应用于当前帧的能力。这对于时间序列问题非常有用，因为在时间序列中可能会出现规则的模式。因为我们的问题也是一个图像问题，我们可以从卷积神经网络(CNN)中受益。CNN 使用核函数从图像中提取重要特征。我们使用 Keras 进行了深度学习，Keras 是一种运行在 TensorFlow 之上的高级 API。该包包括 convLSTM2D 层，它类似于 LSTM 层，但输入层和变换都是卷积的。

我们在这项工作中使用的深度学习 convLSTM2D 模型的图。

第一次卷积输入后，我们应用了 4 个 convLSTM2D 层，每层之间都进行了 relu 激活和归一化。最终的 conv2D 层使用 sigmoid 激活和二进制交叉熵损失函数将张量折叠回 2D 图像，这提供了我们对事故发生概率的 64 x 64 2D 阵列的预测。

该层所需的数据的输入形状是具有形状(样本、时间、行、列、通道)的 5D 张量，在我们的情况下是(样本、时间，64，64，1)。因为我们有一个 64×64 的矩阵，并且它是一个二进制分类问题，所以 1 个通道。在这种情况下，时间是时间步长，我们必须选择要使用多少时间点来预测感兴趣的小时。对于这项工作，我们选择使用之前的 6 个小时来预测 7 日。我们选择了 64 个具有 5×5 内核的过滤器。

特征工程

由于所涉及的因素极其复杂，直接预测交通事故是一项极其困难的工作。几乎不可能包括导致事故的所有因素，包括但不限于司机分心、天气状况、交通状况和道路状况。我们将预测事故风险，而不是直接预测事故。目前，我们的事故数据是按小时在张量中时空排列的。每个网格代表一天中该时间该网格内发生的事故数量。要将事故转换为事故风险，我们需要找到过去 x 天内每个单元的平均事故数。对于这项工作，我们选择寻找 6 天内事故的平均数。例如，如果最近 6 天上午 8 点的一个格网单元包含以下事故数:3、0、6、2、4、1，则关注日上午 8 点的事故风险为(3 + 0 + 6 + 2 + 4 + 1)/6 = 16/6 = 2.67 事故/小时。这种转换应用于所有数据，因此建模中使用的数据是 64 x 64 的事故率数组，基于同一时间的事故和过去 6 天的网格单元。

一小时内发生的事故。

n 的事故率是用前 6 天每天的同一小时计算的

同一小时的事故率是使用过去 6 天中一天的同一小时计算的(右图)。对于构建我们的模型来说，事故率比单独的事故数字提供了更多的信息。最终模型总共包含 2，876，737 个可训练参数，当在使用 GPU 的 Google Colab 笔记本中运行 50 个时期时，总训练时间约为 6 小时。

该模型似乎在第 50 个纪元时收敛，因此我们没有将训练时间延长到这之后。

评估模型

虽然 Keras 提供了许多内置分析，但为了简单起见，我们计算了自己的混淆矩阵。如前所述，该模型输出一个 64 x 64 的事故风险概率数组。

纽约市事故发生地点概率的基本事实(左)和预测(右)。

因为我们想要二进制分类而不是概率，所以我们需要将结果分层为 0 或 1。我们使用 predict_classes 函数来实现这一点，该函数为我们进行了分层。任何 0.5 或以上的概率变成 1，其余的为 0。

纽约市交通事故地点的地面实况(左)和预测(右)。

每个预测和地面实况 64 x 64 阵列被转换成 pandas 数据帧，该数据帧允许逐元素的乘法和减法。当两者相乘时，因为它们都是二进制矩阵，所以只有正确识别的值或真正值会保留下来。当我们从预测中减去地面真实矩阵时，真正值变成零(1 -1)，错误预测的事故或假正值变成 1(1-0)，错过的预测或假负值变成-1(0-1)。然后，可以通过从矩阵的大小(64 x 64 = 4096)中减去真正值、假正值和假负值的总和来计算真负值。

用于构建混淆矩阵的基础事实(右上)、预测(右上)、乘积图(左下)和差异图(右下)。

上面是四幅图像，分别显示了事故风险增加的基本事实(左上)、预测矩阵(右上)、显示正确预测像元的乘积矩阵(左下)以及差值图(右下)，如果像元被错误预测，则差值图返回-1(白色)，如果像元未被预测，但存在于基本事实(右)中，则差值图返回+1(黑色)。由于这两个阶层之间的不平衡(负值远远大于正值，这是一件好事！)，报告模型的准确性不是一个合适的指标。相反，我们可以使用马修相关系数，或 phi 系数作为二进制预测质量的衡量标准。

MCC = TP * TN—FP * FN/sqrt((TP+FP)(TP+FN)(TN+FP)*(TN+FN))

其中 TP，TN，FP，FN 分别=真阳性，真阴性，假阳性，假阴性。应用此公式将返回一个介于 1 和-1 之间的值，其中 1 表示基本事实和预测完全一致，0 表示随机机会，而-1 表示模型和基本事实完全不一致。

跨训练集(左)、测试集(中)和验证集(右)计算的迈克尔相关系数。

计算我们三个集合的迈克尔相关系数，我们看到总体上非常好的分数，训练集、测试集和验证集的平均值分别为 0.82、0.82 和 0.83。令人高兴的是，由于这些数据根本不包括在模型训练中，所以坚持组也与训练和测试组得分相似。然而，我们可以从上面图中较低值的重复模式和下面直方图中的负偏斜看出，我们的预测中存在一些异质性。

每个机组的 MCC 值分布。

理想情况下，MCC 直方图将是平均值周围的正态分布。这表明我们的模型在预测一天中特定时间的事故风险方面较弱。如果它在预测某个区域方面很弱，我们不会在 MCC 得分的急剧下降中看到这一点。

首先，我们在混淆数据帧中添加了一个日期时间索引。然后，我们将阈值设置为 0.7，并过滤出发生这些情况的小时的索引值。绘制这些值的直方图，我们看到绝大多数值出现在凌晨 1 点到 3 点之间。

训练组(左)、测试组(中)和验证组(右)的 MCC 值低于 0.7 的时间分布。

如果我们为我们的模型重新计算平均 MCC 值，如果我们排除这些数据点，我们的训练集、测试集和验证集的平均值都略微增加到 0.84。

去除低于 0.7 的值后，MCC 在每组中的分布。

结论

我们生成了一个深度学习模型，用于预测纽约市的交通事故风险。该模型基于 Keras convLSTM2D 架构，该架构接受事故率的 2D 阵列，该阵列是过去 6 天每小时每个网格单元的平均值。

我们证明了该模型对未来未知数据的预测性，其迈克尔相关系数为 0.84。

用原始事故数而不是事故率来训练同一模型的尝试未能收敛，导致对事故数量和位置的严重高估。

要成为真正有用的模型，每个网格的分辨率最好是城市街区的分辨率。将时间分辨率增加到比 1 小时更频繁的周期也是令人感兴趣的。这当然会大大增加计算时间。

编辑

如果你感兴趣，这里有我在这项工作中使用的笔记本的链接。第一个用于数据处理，第二个用于convlsm 2d模型。

模拟疫情的传播！

原文：https://towardsdatascience.com/modeling-the-spread-of-a-pandemic-92c64a6c6d8d?source=collection_archive---------50-----------------------

叶夫根尼·切尔卡斯基在 Unsplash 上拍摄的照片

变更数据

科学家如何模拟流行病？政府是如何制定封锁计划的？我们怎么知道曲线是否已经变平了？

编者按: 走向数据科学 是一份以数据科学和机器学习研究为主的中型刊物。我们不是健康专家或流行病学家，本文的观点不应被解释为专业建议。想了解更多关于疫情冠状病毒的信息，可以点击 这里 。

我们正生活在我们有生之年从未想象过的最前所未有的时代。我们的生活早已偏离了过去的常态，不幸的是，这个世界再也不会像从前一样了。与其强调已经发生的变化，不如接受新常态，继续前进，重新开始，这将是一个明智的想法。

19 世纪见证了永远改变世界的三件大事——工业化、互联网和传染病！

好吧，让我们来解决房间里的大象——源自中国的冠状病毒感染了全球近 2200 万人，并造成 77 万人死亡。这里是截至 2020 年 8 月 17 日的世界现状。

截图，摄于 2020 年 8 月 17 日 https://www.worldometers.info/coronavirus/

这是这篇文章的议程。您可以选择跳到以下您喜欢的任何部分:

• 新冠肺炎数据可视化
• 指数增长，逻辑增长？他们到底在说什么？
• 了解疫情的基本 SIR 模型和构建模块
• 各国政府如何提出不同的策略？

本文中使用的所有代码都可以从 GitHub 库下载

新冠肺炎数据可视化

数据可视化#1:按国家列出的案例总数

下面的数据可视化贯穿了从 2020 年 2 月到 2020 年 7 月底的时间线，并显示了报告病例总数排名前 20 位的国家是如何移动位置的。

数据可视化#2:按国家划分的死亡总数

另一个数据可视化显示了前 20 个国家在报告的新冠肺炎死亡总数方面的位次变化。

数据可视化#3:显示第一例报告病例前天数的世界热图

下面的数据可视化显示了它传播到全世界的速度。

距首次报告病例还有几天。互动版点击这里

从上面的交互式数据可视化中，我们可以看出:

• 泰国是第一个受到影响的国家(前 13 天)，然后是日本(前 15 天)和韩国(前 20 天)。
• 一个月内，病毒传播到印度、东南亚国家、澳大利亚、美国和加拿大。
• 除了葡萄牙和爱尔兰，大多数西欧国家在头一个月内感染了这种病毒，而东欧国家在大约两个月内都是安全的。
• 俄罗斯在第一个月左右就被感染了。
• 在斯堪的纳维亚国家中，挪威是最后一个被感染的国家，而瑞典和芬兰是在前一个月内被感染的。
• 相对来说，传播到南美和非洲大陆需要更长的时间。
• 朝鲜——没有数据，老规矩，别问我为什么？😛

数据可视化#4:显示前 1000 个案例的天数的世界热图

现在，让我们看看在每个国家达到第一个 1000 个病例需要多长时间:

前 1000 例，互动版，点击此处

从上面的交互式数据可视化中，我们可以看到，要达到第一千例:

• 虽然泰国和日本是第一个被感染的国家，但是他们花了将近两个月的时间才达到第一个千例的里程碑。
• 土耳其、伊朗和塔吉克斯坦只用了两周时间。
• 韩国用了一个多月。
• 大半个欧洲用了近一个月。
• 美国、印度和澳大利亚用了将近 2 个月。
• 巴西和大多数南美国家用了不到一个月的时间。
• 南非用了不到一个月。

尽管这种病毒在一些国家传播相对较慢，但在其他一些国家却传播很快。不同地区的增长率可以用指数增长来表示。

在下一节中，我们将讨论用于病毒建模的不同增长模型，我们还将讨论政府如何提出各种缓解策略。

指数增长，逻辑增长？他们到底在说什么？

在过去的几个月里，我打赌你可能已经在新闻频道上无数次听到这些术语了。他们用这些术语来解释病毒是如何传播的，等等。但是你可能想知道他们到底在说什么？

好吧，为了理解指数增长到底是怎么回事，让我们来处理一个假设的情况。比方说，我让你在我拥有的 Airbnb 网站上住一整月，这是一个惊人的交易——开始的第一天你必须支付 1 美分，但每隔一天，直到一个月，你必须向我支付前一天支付的三倍。这是一笔难以抗拒的交易吗？我看起来像个白痴吗？没关系，你可以有你无声的笑的时刻，但如果你回来报名，我打赌你非常需要一堂数学课。

我用的是所谓的 指数增长 骗你进去。随着指数增长，我的每日利率从第一天的 0.01 美元上升到一个月底的 2 万亿美元以上。如果我们修改我们的交易，所以你每天只付我两倍的钱，你最后一天的最终账单会减少到 1000 万美元，如果我们再次修改交易，所以你每天只增加 50%的钱，你的最终账单会是更合理的 1917.51 美元。

下表总结了这些结果:

指数增长的示例

下面是我们用来计算上表中数字的简单指数方程:

指数增长公式

在上面的公式中，x(t)是你在第 t 天支付给我的，x₀是初始值(你在第一天支付给我的)，r 是增加的百分比。

明白了吗？

科学家和流行病学家利用新冠肺炎的这种指数增长来了解病毒是如何传播的。换句话说，他们试图将每天记录的病例数据点拟合到这个模型中，并试图估计增长因子，从而预测未来。更具体地说，我们可以找出病毒在一个国家相对于另一个国家传播的速度。

但是指数增长有个问题！如果你还没有注意到，增长并没有停止。只要你输入一些参数 t，就会有一些输出值。但这并不是以前所有病毒爆发的情况，因为我们知道它们都在某个时候停止了。

我们的下一位嘉宾来了， 物流成长 ！

与指数增长不同，逻辑增长不会无限增长。

它遵循三个阶段——在第一阶段非常缓慢地增长，在中间阶段非常迅速，在第三阶段再次开始缓慢增长，最后变平。这正是任何病毒爆发都会遵循的模式。

这是解释逻辑增长的简单公式。

逻辑增长公式

科学家们还研究了一种叫做 冈珀兹增长 的东西，这是一种特殊的逻辑增长，也类似于一个 sigmoid 函数。与在第一阶段和最后阶段完全对称的逻辑增长不同，Gompertz 函数增长缓慢，并缓慢变平。下面的公式表示 Gompertz 增长

Gompertz 增长公式

在 logistic 和 Gompertz 增长公式中，参数 c 代表 logistic 函数的渐近线。换句话说，c 的值是疫情结束前我们将出现的最大病例数。

下面的图表给你一个在所有三种模型中增长率如何变化的感觉。

说明指数、逻辑、Gompertz 的增长率

科学家们使用指数函数和逻辑函数来研究我们目前处于流行病的哪个阶段。流行病的第一阶段可以用指数函数更好地表示，然而，第二阶段可以用逻辑函数更好地表示。

利用约翰·霍普斯金大学编制的新冠肺炎报告数据，我们可以对指数和逻辑斯蒂模型进行曲线拟合。

以下是我们将要做的事情——我们将把除 lat 1 day(2020 年 7 月 29 日)之外的所有数据拟合到模型中，并尝试预测这一天的情况。我们将使用[scipy.optimize.curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)一个 Scipy 函数，该函数使用非线性最小二乘回归来拟合一个选定国家的数据。我们可以使用 R2 分数，这是一个测试模型拟合度的常用指标，衡量模型预测曲线新点的能力。值非常接近 1，意味着模型可以很好地预测数据，而值接近 0 甚至负值，意味着模型非常不适合。

下面是相同的代码:

曲线拟合代码

以下是一些选定国家的结果:

中国

从下面的结果来看，中国的指数模型 R2 得分为负，逻辑模型 R2 得分为正。这标志着中国已经越过指数阶段，进入逻辑阶段。R2 分数非常接近 1 表明非常适合逻辑模型。事实上，曲线的平坦化发生在 2020 年 3 月 13 日。

**** Prediction for China-CHN as of 2020-07-29 ****
Initial no. of cases on 2020-01-19 -> 216.0
No. of cases on 2020-07-28 -> 86783.0
~~ Exponential Model ~~
['p=0.035', 'sigma_p=0.001']
MSLE: 45.617, Exp of RMSLE: 9.670, R2 score: -9.730
~~ Logistic Model ~~
['a=10.000', 'b=0.113', 'c=84166.683'] ['sigma_a=1.115', 'sigma_b=0.005', 'sigma_c=390.112']
MSLE: 0.418, Exp of RMSLE: 1.909, R2 score: 0.963
Day of flattening of the infection curve ->  13, Mar 2020

韩国

韩国也是如此。曲线的平坦化发生在 2020 年 5 月 1 日。

**** Prediction for South Korea-KOR as of 2020-07-29 ****
Initial no. of cases on 2020-02-21 -> 155.0
No. of cases on 2020-07-28 -> 14203.0
~~ Exponential Model ~~
['p=0.032', 'sigma_p=0.001']
MSLE: 31.421, Exp of RMSLE: 9.670, R2 score: -11.591
~~ Logistic Model ~~
['a=3.858', 'b=0.074', 'c=12089.302'] ['sigma_a=0.468', 'sigma_b=0.005', 'sigma_c=113.041']
MSLE: 0.132, Exp of RMSLE: 1.439, R2 score: 0.881
Day of flattening of the infection curve ->  01, May 2020

法国和意大利

对法国和意大利来说，逻辑模型显示了拟合优度。这两个国家的曲线在 6 月变平。

**** Prediction for France-FRA as of 2020-07-29 ****
Initial no. of cases on 2020-03-01 -> 100.0
No. of cases on 2020-07-28 -> 183079.0
~~ Exponential Model ~~
['p=0.056', 'sigma_p=0.001']
MSLE: 37.960, Exp of RMSLE: 9.670, R2 score: -3.601
~~ Logistic Model ~~
['a=10.000', 'b=0.055', 'c=166804.289'] ['sigma_a=1.059', 'sigma_b=0.003', 'sigma_c=1840.738']
MSLE: 1.355, Exp of RMSLE: 3.203, R2 score: 0.961
Day of flattening of the infection curve ->  21, Jun 2020

**** Prediction for Italy-ITA as of 2020-07-29 ****
Initial no. of cases on 2020-02-24 -> 132.0
No. of cases on 2020-07-28 -> 246286.0
~~ Exponential Model ~~
['p=0.054', 'sigma_p=0.001']
MSLE: 43.168, Exp of RMSLE: 9.670, R2 score: -3.856
~~ Logistic Model ~~
['a=10.000', 'b=0.056', 'c=243030.882'] ['sigma_a=0.895', 'sigma_b=0.002', 'sigma_c=2068.458']
MSLE: 1.156, Exp of RMSLE: 2.931, R2 score: 0.973
Day of flattening of the infection curve ->  13, Jun 2020

新西兰

新西兰政府的反应被认为是最好的反应之一。对新西兰来说，2020 年 4 月 18 日开始变平

**** Prediction for New Zealand-NZL as of 2020-07-29 ****
Initial no. of cases on 2020-03-23 -> 102.0
No. of cases on 2020-07-28 -> 1207.0
~~ Exponential Model ~~
['p=0.024', 'sigma_p=0.001']
MSLE: 18.770, Exp of RMSLE: 9.670, R2 score: -22.147
~~ Logistic Model ~~
['a=5.567', 'b=0.213', 'c=1159.461'] ['sigma_a=0.323', 'sigma_b=0.006', 'sigma_c=2.736']
MSLE: 0.005, Exp of RMSLE: 1.072, R2 score: 0.985
Day of flattening of the infection curve ->  18, Apr 2020

美国、印度和巴西

当世界其他国家开始看到隧道尽头的一些曙光时，美国、印度和巴西仍然深陷疫情。新病例持续上升。

虽然美国的逻辑模型表明非常适合(R2 接近 1)，但印度和巴西的模型有些不太适合，但仍被认为处于逻辑阶段，因为 R2 > 0。下面的预测表明，这三个国家有望在 2020 年 11 月底之前使曲线变平。

**** Prediction for Brazil-BRA as of 2020-07-29 ****
Initial no. of cases on 2020-03-15 -> 121.0
No. of cases on 2020-07-28 -> 2442375.0
~~ Exponential Model ~~
['p=0.080', 'sigma_p=0.000']
MSLE: 32.941, Exp of RMSLE: 9.670, R2 score: -0.713
~~ Logistic Model ~~
['a=9.998', 'b=0.025', 'c=2159468.826'] ['sigma_a=2.364', 'sigma_b=0.006', 'sigma_c=637699.122']
MSLE: 5.820, Exp of RMSLE: 11.161, R2 score: 0.735
Day of flattening of the infection curve ->  19, Nov 2020

**** Prediction for India-IND as of 2020-07-29 ****
Initial no. of cases on 2020-03-17 -> 125.0
No. of cases on 2020-07-28 -> 1483156.0
~~ Exponential Model ~~
['p=0.076', 'sigma_p=0.000']
MSLE: 27.423, Exp of RMSLE: 9.670, R2 score: -0.560
~~ Logistic Model ~~
['a=9.994', 'b=0.024', 'c=1066388.221'] ['sigma_a=3.159', 'sigma_b=0.008', 'sigma_c=441540.841']
MSLE: 6.032, Exp of RMSLE: 11.657, R2 score: 0.650
Day of flattening of the infection curve ->  27, Nov 2020

**** Prediction for United States-USA as of 2020-07-29 ****
Initial no. of cases on 2020-03-03 -> 103.0
No. of cases on 2020-07-28 -> 4290263.0
~~ Exponential Model ~~
['p=0.079', 'sigma_p=0.001']
MSLE: 45.859, Exp of RMSLE: 9.670, R2 score: -1.501
~~ Logistic Model ~~
['a=10.000', 'b=0.024', 'c=4262579.365'] ['sigma_a=1.038', 'sigma_b=0.003', 'sigma_c=503190.213']
MSLE: 5.149, Exp of RMSLE: 9.670, R2 score: 0.903
Day of flattening of the infection curve ->  17, Nov 2020

理解疫情的基本 SIR 模型和构建模块

虽然指数模型和逻辑模型有助于我们了解我们处于流行病的哪个阶段，但是也有一定的局限性:

• 他们不考虑人口规模，除非我们为这些函数创建一些界限。
• 他们只模拟感染，但不考虑康复和死亡。
• 他们没有考虑感染率或恢复率。
• 他们没有考虑改善治疗对更快康复的影响。
• 他们不考虑疫苗接种的影响，如果在一段时间内有疫苗接种的话。

为了分析如此复杂的场景和模拟，科学家和流行病学家使用了另一种模型，称为 SIR(易感-感染-康复)模型。

这个模型由克马克和麦肯德里克于 1927 年首次发表。在这个模型中，一个地区的全部人口被分成三个部分，称为易感(S)、感染(I)、恢复/移除®。

易感房室包括易受感染的所有人群。在任一时间点，这是总人口减去感染、康复和死亡的人口。

已感染® 区间包括目前已感染病毒且有传播疾病能力的人群。从技术上来说，我们需要至少有一个被感染的个体才能引发流行病。

恢复/移除® 区间包括所有不再被感染且不能感染任何人的人。这包括已经康复、产生免疫力或已经死亡的人。因此，这个群体有时被称为删除。

一般来说，当一个人康复后，他们的身体会产生抗体，从而获得对疾病的免疫力，因此他们不会再被感染。除非病毒能够变异，这可能适用于大多数传染病。

接下来，我们需要了解病毒是如何传播的。

了解病毒如何传播的第一个基础是知道病毒的 繁殖数(R0) 。

生殖数量是一个感染者在患病过程中会感染的平均个体数量。

如果 R0=2.0，这意味着平均每个感染者会将疾病传播给另外两个人。对于任何要成为成功疫情的病毒爆发，它需要具有非常高的 R0 值。

R0 始终作为一个范围发布，因为它根据各种因素不断变化。以下是一些常见病毒爆发的 R0，

一些常见疾病爆发的 R0 值。来源:https://SPH . umich . edu/pursuit/2020 posts/how-scientists-quantify-outbreaks . html

从上表可以看出，科学家估计新冠肺炎的 R0 在 1.5-3.5 之间。它不像一些致命的流行病如麻疹、天花等那样糟糕，但规模如此巨大的原因可能是传播的性质。

任何流行病要减缓传播并逐渐消失，都需要维持 R0 ≤1

下一个构件是病毒的 恢复时间 。

恢复时间是个人从病毒中恢复的预期天数。

缩短恢复时间也会减缓病毒的传播。

使用这两个构建模块，我们可以计算以下参数:

• Gamma (γ): 这代表日恢复数。假设一个人需要 14 天恢复，那么γ = 1/14
• β/感染率(β): 这是病毒的每日繁殖数。假设如果 R0=4，那么β=4/14

下面是一个基本的 SIR 模型的图解，带有来自人口的每个状态的转移率。

随着时间的推移，人们在易感到感染到移除之间过渡。每天都会发生以下情况:

• 易感人口因新增感染者数量而减少。
• 感染人口随着新感染人数的增加而增加，随着新恢复人数的减少而减少。
• 恢复/移除人口随着新恢复的数量而增加。

这些变化率可以用下面的微分方程进行数学表达:

SIR 模型微分方程

使用下面的代码，我们可以为 1000 人的小群体模拟一个简单的 SIR 模型，有 10 个初始感染和 0 个恢复。

SIR 模型

这是简单模拟的结果。

SIR 模型的示例模拟。

可以看出，随着时间的推移，会发生以下情况:

• 在流行病开始的时候，人群中易感个体的比例非常高。这将增加感染率，因为更多来自感染者的相互作用将与易感个体发生。此时，恢复的个体数量非常少**，因为没有多少被暴露。**
• 在疫情中期，易感个体在人群中的比例开始下降。很快新恢复的个体每天变得大于新感染的个体，并且感染**的数量开始减少。****
• ****在一场疫情的结束时，在 附近【R】***数量开始趋于平稳【N】***随着大家都有了病毒，又恢复了 **。感染人数 和 易感人数(S)随着不再有更多的人易感，在 0 附近开始趋于平稳。

各个政府如何提出不同的策略？

在上一节中，我们已经看到了针对小规模人群的基本 SIR 模型的简单模拟，现在让我们来看看大规模的基本 SIR 模型。

让我们考虑一个假设的情况，一个拥有 800 万人口的大都市被 R0=5 的疫情困住，当地政府面临控制局面的压力。

这是这种疾病的简要介绍:

• 再现数，R0=5.0
• 被感染的正常恢复时间是 14 天。
• 初始感染人数为 1000 人
• 初始恢复人数为 0

政府想看到的第一件事是这件事有多大。该镇最优秀的科学家通过一些快速模拟(见下图)显示，大约 375 万人(该市 47.8%的人口)将在疫情爆发的 37 天内被感染。换句话说，什么都不做的代价是非常非常巨大的！

说明疫情如何在大都市蔓延。

因此，政府对即将发生的事情有所察觉。所以他们接下来想的是我们如何控制局面？

科学家们通过一些额外的模拟来显示 R0 的影响，他们强调 R0 需要降至 1 或更低，以减缓病毒的传播，并最终阻止疫情。

毕竟这不是火箭科学。通过查看下面的图表，政府将会明白，降低 R0 会产生巨大的影响——R0 每降低 1 个单位，感染的峰值就会变小，并推迟几周，因此政府可以争取一些时间来改善医疗设施或实施策略或开发疫苗(如果可行的话)。

R0 对疫情传播的影响

政府不希望经济关闭，所以他们想看看某种洗手或公共行为政策限制是否会有所帮助？

科学家们进行了一些模拟，并做出了以下评论:

• 这种病毒的本质是通过感染者和易感人群之间成功的相互作用来传播。
• 洗手和公共接触政策的改变只是通过限制成功互动的次数来降低病毒传播的可能性。因此，R0 只是减少到某个数值，比如说 2.5。
• 即使在这种情况下，至少有 23%的人口在 88 天后被感染。所以这一次，在感染达到高峰之前有一些延迟。

因此，政府明白这是行不通的，所以拒绝了这个想法。

说明洗手指南的效果。

政府能想到的下一件事是封锁，但他们需要知道以下几点:

• 他们需要多长时间实施锁定—1 周、2 周还是 3 周？
• 他们应该封锁多久——30 天还是 60 天？
• 封锁的严格程度应该是怎样的——适度还是严格？

除了寻找上述问题的答案，政府还需要优化以下内容:

• 最大限度地减少感染高峰——确保总体感染比例非常低。
• 我们需要选择一个最佳的锁定选项，以免对经济造成太大伤害，因为一个真正受到伤害的经济稍后会卷土重来，对我们造成同样严重的打击。
• 尽可能争取一些时间，以便扩大医疗设施，或许也给科学家们一些时间来开发疫苗或其他东西。

同样，科学家通过对各种场景进行模拟来帮助政府，并得出以下结果:

不同锁定选项的图示

你应该如何阅读上面的图表？

• 每一行代表锁定开始的日期，例如，第一行的所有图表对应于在爆发后第 7 天开始的锁定，意味着在 1 周内，第二行是关于在第 2 周开始的锁定，依此类推。
• 在每一行中，有 30 天与 60 天、适度与严格的模拟选项。
• 每个图表显示 R0、感染高峰日的时间以及感染高峰日的感染人口比例。
• 当然，我们希望感染比例最低，但与此同时，我们也需要了解感染高峰期还有多长时间。

我们可以从上面的图表中观察到以下情况:

• 在每个锁定选项中，感染的高峰只是被延迟了，但不会完全消失。
• 与相应的适度锁定选项相比，所有严格锁定选项都为我们赢得了更多时间。例如，如果我们考虑从第 14 天开始的 30 天封锁，严格封锁会在 75 天内导致 47%的感染，而中度封锁会在 60 天内导致 45%的感染(相对更快)——赢得一些时间至关重要。

不管潮水有多高，我们宁愿晚一点被打，因为这样我们可以做一些安全准备或制定逃生路线。

• 在所有中度封锁中，封锁解除后，感染高峰就会上升。
• 等到 35 天后才实施封锁看起来太晚了，因为感染高峰会上升，所以完全排除了这种选择。

通过考虑前面讨论的所有优化因素，政府提出了一个从第 14 天开始实施的 30 天严格封锁计划。

我们现在正处于从第 14 天开始的为期 30 天的严格封锁之中。利用赢得的时间，政府在两周内迅速增加了医疗设施，并从第 28 天开始部署改进的治疗方案。

以下是正常锁定方案与改进处理方案的对比:

改进了 30 天锁定中的恢复时间图示

从上面的图表中我们可以观察到以下情况

• 改进的治疗方案大大减少了 50%的恢复时间，这是非常好的。
• 感染高峰减少了 50%，这意味着现在只有 22%的人被感染，这是一个好消息。
• 感染高峰推迟了，这意味着我们可以再争取几天时间——这太好了。

随着改进的治疗方案的到位，政府想知道如果我们现在将封锁期延长到 60 天，会对感染高峰产生任何影响吗？

下图显示了相同的改进治疗方案，但是是在 60 天严格封锁的情况下。

改进了 60 天锁定中的恢复时间图示

30 天的禁闭与 60 天的禁闭有什么不同？

尽管在一天结束时感染的人口比例相同，但是通过 60 天的封锁，我们赢得了更多的时间(135 天)，这对制定下一步措施至关重要。

因此，政府将封锁期延长至 60 天。

政府现在意识到永远推迟这一趋势是没有意义的，因此它与一组科学家合作，为疫苗开发提供一些资金。

如果有疫苗呢？下一个问题是有效控制感染的每日容量应该是多少？

下面显示了 60 天锁定环境中的另一个模拟:

疫苗容量说明

这些数字看起来很有希望。

对于我们正在应对的大都市(800 万人口)的规模，我们似乎需要每天接种 3 万到 5 万次疫苗，以保持较低的感染率。

因此，现在政府可以致力于提供必要的基础设施，希望疫苗能尽快问世。

在本文中，为了简单起见，我们只讨论了基本的 SIR 模型。但是，该模型没有考虑以下情况:

• 生命动力学——每天都有新生儿出生和死亡。
• 移民每天都在发生。
• 大多数病毒都有潜伏期——这可以用 SEIR 模型来模拟，它是 SIR 的变体，引入了另一个称为“暴露”的隔室。
• 有些人是渐近的，不处于感染状态等。
• 住院、重症监护室、呼吸机等。

在真实世界场景中，科学家采用最准确的数据，例如新记录的出生或死亡数量、因其他原因发生的死亡数量、ICU 中的人数等，从而做出准确的预测。

因此，在本文中，我们看到了如何将新冠肺炎数据拟合到指数和逻辑模型中，从而做出预测，并推断出某个国家何时曲线变平，或者如果没有变平，预计何时会变平。

此外，我们还看到科学家如何使用 SIR 模型模拟各种场景，以帮助政府制定不同的病毒预防策略。

本文使用的所有代码都可以从GitHub 资源库下载