武汉加油，中国加油

基本再生数

**基本再生数（basic reproduction number，R0）**是指在一个全是易感染态个体构成的群体中，一个感染态的个体在恢复之前平均能感染的人数。在流行病学中，R0>1 表示疾病将爆发，<1 则表示疾病走向消亡，故 R0 是判断流行病是否爆发的重要条件之一。但在真实疫情的传播过程中，由于政府干预政策实施、个体行为改变（戴口罩、减少出行等）、易感人群数量减少（因患病人数增多或使用疫苗等）等外在因素影响，不再满足基本再生数定义的理想模型条件。为刻画流行病随时间的演化过程，我们将传播过程中某一时刻下（t）平均一个感染态个体感染的人数定义为有效再生数，记为 Rt。在实际疫情的控制过程中，当 Rt<1 时，即平均一个感染态个体感染的人数小于 1 时，认为疾病已被控制同时疾病也将走向消亡。

SIR模型

传染病模型有着悠久的历史，一般认为始于1760年Daniel Bernoulli在他的一篇论文中对接种预防天花的研究。真正的确定性传染病数学模型研究的前进步伐早在20世纪初就开始了，Hamer、Ross等人在建立传染病数学模型的研究中做出了大量的工作，直到1927年Kermack与McKendrick在研究流行于伦敦的黑死病时提出了的SIR仓室模型，并于1932年继而建立了SIS模型，在对这些模型的研究基础上提出了传染病动力学中的阈值理论。Kermack与McKendrick的SIR模型是传染病模型中最经典、最基本的模型，为传染病动力学的研究做出了奠基性的贡献。
模型中把传染病流行范围内的人群分成三类：S类，易感者（Susceptible），指未得病者，但缺乏免疫能力，与感病者接触后容易受到感染；I类，感病者（Infective），指染上传染病的人，它可以传播给S类成员；R类，移出者（Removal），指被隔离，或因病愈而具有免疫力的人。
来自百度百科

数据来源

https://voice.baidu.com/act/newpneumonia/newpneumonia/?from=osari_pc_3
以及每日深圳公布的详细感染病例

代码部分

# from login.mysql_login import conn as CONN
import pandas as pd,numpy as np
from scipy.integrate import odeint
import math,datetime
import redef getCureTime(cure_time:str):cure_time = cure_time.split('。')[-2]cure_time = cure_time.replace('已','')cure_time = cure_time.replace('痊愈出院','')cure_time = cure_time.replace('治愈出院','')cure_time = re.sub('（.+?）','',cure_time)if cure_time == '目前情况稳定':return Noneif cure_time == '23日':cure_time = '1月' + cure_timecure_time = '2020年' + cure_timecure_time = datetime.datetime.strptime(cure_time, "%Y年%m月%d日")return cure_timeif __name__ == '__main__':# sp_ncov_status_all = pd.read_sql('''# select  occur_period_f,#         admin_area_name,#         confirm as '确诊病例(人）',#         cure as '治愈人数（人）',#         death as '死亡人数（人）',#         suspect as '疑似病例（人）'# from qhdata_support_spider.sp_ncov_status# WHERE admin_area_code = 0086# ORDER BY occur_period_f# ''',CONN)## sp_ncov_status_sz = pd.read_sql('''# select occur_period_f,#        admin_area_name,#        confirm as '确诊病例(人）',#        cure as '治愈人数（人）',#        death as '死亡人数（人）',#        suspect as '疑似病例（人）'# from qhdata_support_spider.sp_ncov_status# WHERE admin_area_code = 440300000000#     ''', CONN)sp_ncov_status_all = pd.read_excel('source/sp_ncov_status_all.xlsx')sp_ncov_status_sz = pd.read_excel('source/sp_ncov_status_sz.xlsx')detail_case = pd.read_excel('source/深圳市冠状病毒病例个案数据（汇总）-20200211.xlsx')# ##############################数据预处理##############################################for i in ['确诊病例(人）','治愈人数（人）','死亡人数（人）']:#     sp_ncov_status_all['当日' + i] = sp_ncov_status_all[i] - sp_ncov_status_all[i].shift()sp_ncov_status_sz['当日' + i] = sp_ncov_status_sz[i] - sp_ncov_status_sz[i].shift()detail_case['疑似到确诊时间间隔'] = detail_case[['发病时间','确诊时间']].\apply(lambda x:None if x[0] == '-' elseNone if x[1] == '-' else(x[1] - x[0]).days,axis=1)detail_case['疑似到确诊时间间隔'] = detail_case['疑似到确诊时间间隔'].apply(lambda x:None if x < 0 else x)detail_case['年龄'] = detail_case['年龄'].apply(lambda x:x.replace('。','') if type(x) == str else x)detail_case['年龄'] = detail_case['年龄'].astype(int)detail_case['治愈时间'] = detail_case[['目前是否治愈','完整病例介绍原文']]\.apply(lambda x:getCureTime(x[1]) if x[0] == '是' else None,axis=1)detail_case['恢复天数'] = detail_case[['确诊时间','治愈时间']].\apply(lambda x:None if x[0] == '-' elseNone if x[1] == '-' else(x[1] - x[0]).days,axis=1)# 以第一个不明原因肺炎发现者的报道时间2019年12月8日为感染的第一天first_date = datetime.datetime.strptime('2019-12-08', "%Y-%m-%d")sp_ncov_status_all['t'] = sp_ncov_status_all['occur_period_f'].apply(lambda x:(datetime.datetime.strptime(x, "%Y-%m-%d") - first_date).days + 1)# ###################################################################################### ##############################数据统计################################################case_num = detail_case.shape[0]# 性别比例sex_info = pd.pivot_table(detail_case,index=['性别'],values=['病例序号'],aggfunc={'病例序号': [len,lambda x:len(x)/case_num]}).reset_index()sex_info.columns = ['性别','占比','人数']# 年龄分层temp_df = detail_case[['年龄','病例序号']]temp_df['年龄分组'] = temp_df.年龄.apply(lambda x:'1.青少幼年(20以下)' if x < 20 else'2.青壮年(20-50)' if ((x >= 20) and (x < 50)) else'3.中年(50-70)' if ((x >= 50) and (x < 70)) else'4.老年(70以上)')age_info =temp_df[['年龄分组','病例序号']].\groupby('年龄分组').\agg({'病例序号':[len,lambda x:len(x)/case_num]}).reset_index()age_info.columns = ['年龄分组','人数','占比']# 疑似到确诊时间间隔diagnostic_interval = pd.DataFrame(detail_case.describe()['疑似到确诊时间间隔']).Tdiagnostic_interval_1 = pd.pivot_table(detail_case,index=['疑似到确诊时间间隔'],values=['病例序号'],aggfunc={'病例序号': [len,lambda x:len(x)/case_num]}).reset_index()diagnostic_interval_1.columns = ['疑似到确诊时间间隔','占比','人数']# 恢复天数restore_day = pd.DataFrame(detail_case.describe()['恢复天数']).Trestore_day_1 = pd.pivot_table(detail_case,index=['恢复天数'],values=['病例序号'],aggfunc={'病例序号': [len,lambda x:len(x)/case_num]}).reset_index()restore_day_1.columns = ['恢复天数','占比','人数']# ###################################################################################### ##############################R0计算逻辑##############################################avg_day = diagnostic_interval['mean'][0].round(2) # 疑似到确诊时间间隔引用深圳市的平均值(天)# 选择拥有新增最新数据的2.4-2.6日的确认数据与对应的2.1-2.3疑似数据(相隔4天左右),计算转化率pp = round((3887 + 3694 + 3151) / (4562 + 5173 + 5072), 2) # 这里计算出来的p大约0.72左右# SARS传播的p的取值在0.5-0.8之间;初期报导提到59个疑似病例有41位最终确诊，因此p的一个参考值为41/59=0.695;# 以第一个不明原因肺炎发现者的报道时间2019年12月8日为t=1# 2020年2月6日实际上被感染人数为Y(t),t=61Y_t = 3151 + round(4833*p,0)t = 61T_g  = 8.4# 生成时间（generation time）T_i = (5+6+4+3)/4 # 潜伏时间(incubation time)lamda = math.log(Y_t,math.e)/tRho = T_i/T_g # 潜伏期占生成时间的比例Rho# 论文里面说平均潜伏期5.1天，不过我没有找到。钟南山说3天。R_0 = round(1 + lamda*T_g + Rho*(1 - Rho) * math.pow((lamda*T_g),2),2)# ###################################################################################### ##############################R0数据计算##############################################def basicReproductionNumber(Y_t,t,T_i,T_g  = 8.4):lamda = math.log(Y_t,math.e)/tRho = T_i/T_g  # 潜伏期占生成时间的比例RhoR_0 = 1 + lamda * T_g + Rho * (1 - Rho) * math.pow((lamda * T_g), 2)return R_0sp_ncov_status_all_temp = sp_ncov_status_all.dropna()sp_ncov_status_all_temp['Y_t'] = sp_ncov_status_all_temp['新增确诊（人）'] + sp_ncov_status_all_temp['新增疑似（人）'] * psp_ncov_status_all_temp['Y_t'] = sp_ncov_status_all_temp['Y_t'].round(0)sp_ncov_status_all_temp['R0'] = sp_ncov_status_all_temp[['Y_t','t']].\apply(lambda x:basicReproductionNumber(x[0],x[1],T_i,8.4),axis=1)# ###################################################################################### ##############################SIR计算逻辑################################################# https://towardsdatascience.com/modelling-the-coronavirus-epidemic-spreading-in-a-city-with-python-babd14d82fa2# SIR模型中的总人口数, N.N = 11081000 # 这里采用2019年武汉市常住人口数(百度的)# 感染者I与移除者RI0, R0 = 44, 0 # I0选择最初的武汉华南海鲜市场的44个感染者作为最初的感染者# 易感人群SS0 = N - I0 - R0# 恢复率(移除率)gamma# 这里是选取深圳感染案例的平均恢复天数，实际上不同地区的医疗水平不同会导致不同的恢复天数。gamma = 1. / restore_day['mean'][0]# 接触率beta# 这里的R_0是全国的R_0，但是实际上应该将R_0分成湖北省内与全国除湖北省的两种情况。# 现在省内R_0明显比全国的要高的多，而省外明显要比全国的低的多。在此我们额外添加数个R_0来观测[0.9,2,3,4]# R_0 = 4beta = R_0 * gamma # 基本再生数 R0, 等于接触率除以移出率# 最长时间(天)day = 250t = np.linspace(0, day, day)def deriv(y, t, N, beta, gamma):S, I, R = ydSdt = -beta * S * I / NdIdt = beta * S * I / N - gamma * IdRdt = gamma * Ireturn dSdt, dIdt, dRdty0 = S0, I0, R0ret = odeint(deriv, y0, t, args=(N, beta, gamma))S, I, R = ret.T# ###################################################################################### ##############################画图###################################################import matplotlib.pyplot as pltimport seaborn as snsplt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号fig = plt.figure(figsize=[16, 9])  # 设定图片大小temp_df = age_info[['年龄分组','人数']]sns.barplot(x='年龄分组', y='人数',data=temp_df[['年龄分组','人数']])plt.savefig('pic/年龄分层.png')plt.close()fig = plt.figure(figsize=[16, 9])  # 设定图片大小ax1 = fig.add_subplot(111)  # 添加第一副图ax2 = ax1.twinx()  # 共享x轴，这句很关键！sns.barplot(x='性别', y='人数', data=sex_info[['性别','人数']], ax=ax1, alpha=0.8)  # 画柱状图,注意ax是选择画布sns.pointplot(x='性别', y='占比', data=sex_info[['性别','占比']], ax=ax2)  # 画在ax2画布上plt.savefig('pic/性别比例.png')# plt.show()plt.close(fig)fig = plt.figure(figsize=[16, 9])  # 设定图片大小temp_df = diagnostic_interval_1[['疑似到确诊时间间隔','人数']]sns.barplot(x='疑似到确诊时间间隔', y='人数',data=temp_df[['疑似到确诊时间间隔','人数']])plt.savefig('pic/疑似到确诊时间间隔.png')plt.close()fig = plt.figure(figsize=[16, 9])  # 设定图片大小temp_df = restore_day_1[['恢复天数','人数']]sns.barplot(x='恢复天数', y='人数',data=temp_df[['恢复天数','人数']])plt.savefig('pic/恢复天数.png')plt.close()fig = plt.figure(figsize=[16, 9])  # 设定图片大小temp_df = sp_ncov_status_all_temp[['occur_period_f', 'R0']]sns.pointplot(x='occur_period_f', y='R0',data=temp_df)plt.savefig('pic/R0变化趋势.png')plt.close()# SIR模型fig = plt.figure(facecolor='w',figsize=[16, 9],dpi=144)plt.style.available # style样式 https://neusncp.com/user/blog?id=203plt.style.use('seaborn-darkgrid')plt.suptitle('R_0: %s' % R_0, fontsize=16, fontweight='bold')plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号ax = fig.add_subplot(111, axisbelow=True)ax.plot(t, S / 10000, 'b', alpha=0.5, lw=2, label='易感人群')ax.plot(t, I / 10000, 'r', alpha=0.5, lw=2, label='感染者')ax.plot(t, R / 10000, 'g', alpha=0.5, lw=2, label='移除者')ax.set_xlabel('时间(天)')ax.set_ylabel('人数(万人)')ax.xaxis.set_major_locator(plt.MultipleLocator(10))ax.set_xlim(0, day)ax.yaxis.set_tick_params(length=0)ax.xaxis.set_tick_params(length=0)ax.grid(b=True, which='major', c='w', lw=2, ls='-')legend = ax.legend()legend.get_frame().set_alpha(0.5)for spine in ('top', 'right', 'bottom', 'left'):ax.spines[spine].set_visible(False)plt.savefig(f'pic/SIR变化趋势(R_0={R_0}).png')plt.close()# #####################################################################################

部分说明

1.计算公式的来由：
中国循证医学杂志的《武汉新型冠状病毒感染肺炎基本再生数的初步预测》
柳叶刀的《Nowcasting and forecasting the potential domestic and
international spread of the 2019-nCoV outbreak originating
in Wuhan, China: a modelling study》

2.Tg的来由：
柳叶刀的《A familial cluster of pneumonia associated with the 2019
novel coronavirus indicating person-to-person transmission:
a study of a family cluster》

红框是没有去武汉的深圳一家人中的一人，从中可以看出，该人在其他家庭成员回来后大约经过8.4天后发病。出现初步症状是在8号左右。

3.Tg的来由：
柳叶刀的《A familial cluster of pneumonia associated with the 2019
novel coronavirus indicating person-to-person transmission:
a study of a family cluster》

不过在论文《武汉新型冠状病毒感染肺炎基本再生数的初步预测》(http://www.cjebm.com/article/10.7507/1672-2531.202001118)中写的平均潜伏时间是5.1天，具体我在源论文中没有找到。我只找到这个表格有个类似的数据。

结果展示

年龄分组

性别分组

疑似到确诊时间间隔分布

恢复天数

全国R0变化趋势

简易SIR变化趋势

这里的t是以第一个不明原因肺炎发现者的报道时间2019年12月8日为感染的第一天。
SIR模型中的移除者可以使治愈，可以是隔离，也可以是死亡，请注意。
从由于R_0采用的是全国数据来计算的，但是实际上应该将R_0分成湖北省与全国除湖北省两种情况来看待。从现在看全国(除湖北)根本没有爆发大规模感染，所以我假设全国(除湖北)的R_0小于1。而湖北到今日(2月14日，与2019年12月8日间隔69天)其确诊、疑似病患数量仍然上升，可以看出其R_0比全国R_0(2.58)要高，所以这里使用3,4来测试。
从图与实际结合推测，全国(除湖北)的SIR比较接近于R_0=0.9情况，而湖北则从现在仍然是爆发期来看，比较接近R_0=3-4之间.