预测疫情新增确诊数，我是这么做的

平常我的工作和机器学习相关，我看到数据的第一反应就是构建模型来预测，于是我今天花了几个小时做了下这个工作。

下面是两周的全国层面的疫情数据，我们能不能，基于这 14 条数据，预测一下明天的新增确诊人数呢？做到技术与时事的结合。

其中各个字段名称含义如下：

• new：每天新增确诊人数，昨天新增为 2829

• all：当前总共确诊人数

• watch：当前医学观察人数

先说结论吧：

预测明天公布的新增确诊人数为 2997 人，有浮动哈，具体浮动多少我也不敢说，怕脸疼！

分析

说说我是怎么分析的吧，所有确诊人数和新增确诊人数我们看一个就行了，这里我选的是每日新增确诊人数，来看看每日新增确诊人数的折线图（目的是观察相等时间间隔下数据的趋势）

可以看到一个随时间稳步攀升的趋势，这是我们期望预测的指标，也就是我们的因变量。

接下来看看因变量，也就是我们的当前医学观察人数：

再给这两个变量画个图：

  

   

    

     

    
    

     

      plt.figure()
     
    
   

    

     

    
    

     

      plt.plot(df[
      "watch"], df[
      "new"])
     
    
   

    

     

    
    

     

      plt.show()
     
    
  

可以去掉 27 号的数据：

  

   

    

     

    
    

     

      clean_df = df.drop(df.index[
      7])
     
    
   

    

     

    
    

     

      plt.figure()
     
    
   

    

     

    
    

     

      plt.plot(clean_df[
      "watch"], clean_df[
      "new"])
     
    
   

    

     

    
    

     

      plt.show()
     
    
   

    

     

    
    

     

      # 绘图并计算相关程度 高达
      0.9935173395339865
     
    
   

    

     

    
    

     

      clean_df.watch.corr(clean_df.
      new)
     
    
  

这里我们知道的是当天的医学观察数和当天的新增确诊数有较强的线性关系，现在有一个问题是我们并不能通过当天的医学观察数来预测当天新增确诊数，因为医学观察数要到预测的那天才知道。

我是这样解决：

• 观察今日的医学观察数和明日的新增确诊数有什么关系，这样我们只要知道今日的医学观察数就可以进行预测

• 先构建当日的医学观察数和当日的新增确诊数模型，再构建一个模型基于历史的医学观察数预测后面一天的医学观察数，这样就得到了最新的医学观察数，同样可以根据此时的医学观察数预测对应的新增确诊数

先看第一个方案，数据重新排列一下，如下：

绘制当前数据的图像：

  

   

    

     

    
    

     

      plt.figure()
     
    
   

    

     

    
    

     

      plt.plot(new_df[
      "watch"], new_df[
      "new"])
     
    
   

    

     

    
    

     

      plt.show()
     
    
  

目测线性关系没有很强，算一下看看：

  

   

    

     

    
    

     

      # 
      0.9840477846537675 还可以 没上面的好
     
    
   

    

     

    
    

     

      new_df.watch.corr(new_df.
      new)
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      # 对数据进行对数化
     
    
   

    

     

    
    

     

      log_df = np.log(new_df)
     
    
   

    

     

    
    

     

      # 
      0.9931578091568266
     
    
   

    

     

    
    

     

      log_df.watch.corr(log_df.
      new)
     
    
  

不错，可以看到这个方案可行，此方案实现就可以直接预测出今日医学观察数对应的明日新增确诊数。

再看第二个方案，我们已经有了当日医学观察数和当日新增确诊数的模型，现在缺的是当日医学观察数，我们可以根据时间序列的ARIMA模型来根据历史数据预测当日医学观察数。

  

   

    

     

    
    

     

      from statsmodels.tsa.stattools 
      import adfuller as ADF
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      # 并不平稳 p值太大，不平稳
     
    
   

    

     

    
    

     

      ADF(df[
      "watch"])
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      # 对数化 平稳了
     
    
   

    

     

    
    

     

      log_watch_df = np.log(clean_watch_df)
     
    
   

    

     

    
    

     

      ADF(np.log(df[
      "watch"])), ADF(log_watch_df)
     
    
  

就这样，我们确定了预测方案：

• 第一个方案：输入今日医学观察数，得出明日新增确诊数

• 第二个方案：通过时间序列模型预测命题医学观察数，再输入明日医学观察数，输出新增确诊数

开干吧！开始构建模型。

模型

首先构建两个线性回归模型：

• 输入今日医学观察数，预测明日新增确诊数

• 输入明日医学观察数，输出明日新增确诊数

  

   

    

     

    
    

     

      from sklearn.linear_model 
      import LinearRegression
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      # 第一个线性回归模型
     
    
   

    

     

    
    

     

      X_train = clean_df[
      "watch"].values.reshape(
      -1,
      1)
     
    
   

    

     

    
    

     

      y_train = clean_df[
      "new"].to_list()
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      lin_reg01 = LinearRegression()
     
    
   

    

     

    
    

     

      lin_reg01.fit(X_train, y_train)
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      # 第二个线性回归模型 new_df指的是上述重新排列后的数据集
     
    
   

    

     

    
    

     

      X_log_train = np.log(new_df[
      "watch"].values.reshape(
      -1,
      1))
     
    
   

    

     

    
    

     

      y_log_train = np.log(new_df[
      "new"].to_list())
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      lin_reg02 = LinearRegression()
     
    
   

    

     

    
    

     

      lin_reg02.fit(X_log_train, y_log_train)
     
    
  

其中第二个模型可以直接进行预测了：

  

   

    

     

    
    

     

      # 
      0202的医学观察数为
      152700，输入
     
    
   

    

     

    
    

     

      res = lin_reg02.predict(np.log([[
      152700]]))[
      0]
     
    
   

    

     

    
    

     

      # 得出结果为
      2931
     
    
   

    

     

    
    

     

      print(np.exp(res))
     
    
  

这里可以看到我们已经得出第一个结果了2931，我还是想试试第二个方案的结果，2931可以作为一个参考，后面可以将两个结果平均一下。

第二个方案，线性回归模型我们已经构建好了，接下来弄个时间序列模型来预测后面一天的医学观察数。

  

   

    

     

    
    

     

      from statsmodels.tsa.arima_model 
      import ARIMA
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      order = sm.tsa.arma_order_select_ic(log_watch_df, max_ar=
      7, max_ma=
      7, ic=[
      'aic', 
      'bic', 
      'hqic'])
     
    
   

    

     

    
    

     

      # 给出最优p q值
     
    
   

    

     

    
    

     

      p, q =order.bic_min_order
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      train_X = log_watch_df[:]
     
    
   

    

     

    
    

     

      arima_model = ARIMA(train_X, (p, 
      0, q)).fit()
     
    
   

    

     

    
    

     

      # 预测未来两天数据
     
    
   

    

     

    
    

     

      # 
      2020
      -02
      -03    
      11.974123
     
    
   

    

     

    
    

     

      # 
      2020
      -02
      -04    
      11.945058
     
    
   

    

     

    
    

     

      predict_data_02 = arima_model.predict(start=
      len(train_X), end=
      len(train_X) + 
      1, dynamic = False)
     
    
   

    

     

    
    

     

      # 预测历史数据
     
    
   

    

     

    
    

     

      # 
      2020
      -01
      -28    
      10.853324
     
    
   

    

     

    
    

     

      # 
      2020
      -01
      -29    
      10.997501
     
    
   

    

     

    
    

     

      # 
      2020
      -01
      -30    
      11.591271
     
    
   

    

     

    
    

     

      # 
      2020
      -01
      -31    
      11.714462
     
    
   

    

     

    
    

     

      # 
      2020
      -02
      -01    
      11.775772
     
    
   

    

     

    
    

     

      # 
      2020
      -02
      -02    
      11.920444
     
    
   

    

     

    
    

     

      predict_data = arima_model.predict(dynamic = False)
     
    
  

看看预测值和实际值：

  

   

    

     

    
    

     

      comp = pd.DataFrame()
     
    
   

    

     

    
    

     

      comp[
      'original'] = train_X[:]
     
    
   

    

     

    
    

     

      comp[
      'predict'] = predict_data[:]
     
    
   

    

     

    
    

     

      comp.plot(figsize=(
      12,
      6))
     
    
  

此时我们已经知道2020-02-03 11.974123是模型预测的明天的值，我们转换一下再将前面的误差综合一下，2020-02-03的预测医学观察数是：163739，根据第一个模型，输入明日医学观察数，得出明日新增确诊人数：

  

   

    

     

    
    

     

      lin_reg01.predict([[
      163739]])
     
    
   

    

     

    
    

     

      # 输出
     
    
   

    

     

    
    

     

      array([
      3063.42234138])
     
    
  

目前可以得出以下结果：

• 对于第一个方案，预测2020-02-03新增确诊人数是：3063

• 对于第二个方案：预测2020-02-03新增确诊人数是：2931

我来取个中间值吧：2997。

于是，我预测2020-02-03新增确诊人数是：2997，此结论仅仅作为参考，也肯定会在一定范围内波动，不要当真哈。

欢迎关注公众号【老胡的储物柜】

历史推文：

• 据非官方研究：学习Python可抑制冠状病毒

• 【2019 nCoV】在微信群自动监控疫情最新进展

• 程序员如何在工作中进步