TowardsDataScience 博客中文翻译 2019（三百九十一）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

利用向量误差修正模型预测方股票价格

原文：https://towardsdatascience.com/predicting-fang-stock-prices-ecac4ddd27c1?source=collection_archive---------10-----------------------

Source: All images in this article were generated by the author in R using dygraphs

基于向量误差修正模型的协整时间序列数据建模

金融数据是时间序列数据最著名的资源。一个多世纪以来，股票价格的记录被一丝不苟地记录下来，提供了丰富的数据集。所谓的“方”公司，即脸书、亚马逊、网飞和谷歌/Alphabet 的股票因其令人印象深刻的回报率和这些高科技公司的“酷”因素而受到了很多关注。但是这些股票真的相互关联吗？或者这是华尔街营销一堆不相关商品的一个例子？还有一种可能是，投资者认为他们之间存在关联，这足以导致他们的股票一起波动。即使我们不知道它们为什么一起移动，如果它们一起移动，我们至少可以确定你是否可以预测一只股票和其他股票。

我从雅虎财经获得了方股票的收盘价数据。脸书的首次公开募股是在 2012 年，所以这是样本期的限制因素。我们的样本期从 2012 年 5 月 18 日到 2019 年 8 月 23 日。我截掉了 2019 年的所有交易日作为测试数据集，把剩下的观测值留给训练数据集。这意味着数据集缺乏这些公司的股票在衰退中如何表现的数据。不管怎样，我继续我的分析。我测试了整合度、格兰杰因果关系和约翰森整合。基于这些结果，我继续进行向量误差修正模型(VECM)，这是一种同时估计至少有一个协整关系的多个时间序列的技术。

平稳性测试

平稳性是指时间序列的均值方差不依赖于时间。当用图表表示时，它看起来像白噪声。由于对非平稳时间序列建模所引起的几个问题，主要是自相关，你需要在建模之前检查一个时间序列是否是平稳的。如果它不是静止的，你需要通过差分得到它。

上图显示了 4 个绝对看起来不稳定的时间序列。但是，在继续之前，我们需要对它进行正式测试。

无差异股票的自相关函数(ACF)图表明，股票肯定不是平稳的。差异变量的 ACF 看起来可能是稳定的。这表明每只股票都是综合 I(1)。增强的 Dickey-Fuller (ADF)测试可以证实这一点。

脸书的自相关函数图

亚马逊的自相关函数图

网飞的自相关函数图

谷歌的自相关函数图

无差异和有差异的 ADF 测试的输出强烈支持股票是整合的 I(1)的假设。在建模之前，需要对这些股票价格进行一次区分。许多用于集成时间序列模型的时间序列包会自动完成这项工作。

针对脸书的增强迪基-富勒测验

Non-Differenced## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 0.405   0.760
## [2,]   1 0.441   0.771
## [3,]   2 0.483   0.783
## [4,]   3 0.508   0.790
## [5,]   4 0.542   0.800
## [6,]   5 0.608   0.819
## [7,]   6 0.639   0.828
## [8,]   7 0.692   0.843
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -1.11   0.661
## [2,]   1 -1.16   0.644
## [3,]   2 -1.19   0.634
## [4,]   3 -1.17   0.642
## [5,]   4 -1.14   0.651
## [6,]   5 -1.14   0.651
## [7,]   6 -1.19   0.635
## [8,]   7 -1.19   0.635
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -2.42   0.399
## [2,]   1 -2.25   0.469
## [3,]   2 -2.08   0.541
## [4,]   3 -2.01   0.572
## [5,]   4 -1.88   0.628
## [6,]   5 -1.60   0.745
## [7,]   6 -1.47   0.803
## [8,]   7 -1.28   0.881
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01Differenced## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -41.5    0.01
## [2,]   1 -29.8    0.01
## [3,]   2 -24.6    0.01
## [4,]   3 -21.7    0.01
## [5,]   4 -20.1    0.01
## [6,]   5 -18.3    0.01
## [7,]   6 -17.3    0.01
## [8,]   7 -16.4    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -41.5    0.01
## [2,]   1 -29.8    0.01
## [3,]   2 -24.6    0.01
## [4,]   3 -21.7    0.01
## [5,]   4 -20.2    0.01
## [6,]   5 -18.4    0.01
## [7,]   6 -17.4    0.01
## [8,]   7 -16.5    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -41.5    0.01
## [2,]   1 -29.8    0.01
## [3,]   2 -24.6    0.01
## [4,]   3 -21.7    0.01
## [5,]   4 -20.2    0.01
## [6,]   5 -18.4    0.01
## [7,]   6 -17.4    0.01
## [8,]   7 -16.5    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

亚马逊增强的迪基-富勒测试

Non-Differenced## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag  ADF p.value
## [1,]   0 1.54   0.969
## [2,]   1 1.60   0.974
## [3,]   2 1.60   0.974
## [4,]   3 1.67   0.977
## [5,]   4 1.76   0.980
## [6,]   5 1.62   0.975
## [7,]   6 1.63   0.976
## [8,]   7 1.75   0.980
## Type 2: with drift no trend 
##      lag      ADF p.value
## [1,]   0 -0.03873   0.953
## [2,]   1  0.00509   0.956
## [3,]   2  0.00206   0.956
## [4,]   3  0.04246   0.959
## [5,]   4  0.08628   0.963
## [6,]   5  0.01241   0.957
## [7,]   6  0.02187   0.957
## [8,]   7  0.07776   0.962
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -1.93   0.605
## [2,]   1 -1.90   0.619
## [3,]   2 -1.89   0.622
## [4,]   3 -1.86   0.637
## [5,]   4 -1.80   0.663
## [6,]   5 -1.87   0.631
## [7,]   6 -1.87   0.634
## [8,]   7 -1.80   0.664
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01Differenced## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -41.9    0.01
## [2,]   1 -29.1    0.01
## [3,]   2 -24.3    0.01
## [4,]   3 -21.7    0.01
## [5,]   4 -17.8    0.01
## [6,]   5 -16.4    0.01
## [7,]   6 -15.9    0.01
## [8,]   7 -16.9    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -41.9    0.01
## [2,]   1 -29.2    0.01
## [3,]   2 -24.4    0.01
## [4,]   3 -21.8    0.01
## [5,]   4 -17.9    0.01
## [6,]   5 -16.5    0.01
## [7,]   6 -16.0    0.01
## [8,]   7 -17.1    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -41.9    0.01
## [2,]   1 -29.2    0.01
## [3,]   2 -24.5    0.01
## [4,]   3 -21.8    0.01
## [5,]   4 -18.0    0.01
## [6,]   5 -16.5    0.01
## [7,]   6 -16.1    0.01
## [8,]   7 -17.1    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

为网飞增加了迪基-富勒测验

Non-Differenced## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 0.837   0.885
## [2,]   1 0.833   0.884
## [3,]   2 0.799   0.874
## [4,]   3 0.701   0.846
## [5,]   4 0.670   0.837
## [6,]   5 0.880   0.897
## [7,]   6 0.890   0.900
## [8,]   7 0.946   0.907
## Type 2: with drift no trend 
##      lag    ADF p.value
## [1,]   0 -0.516   0.873
## [2,]   1 -0.517   0.873
## [3,]   2 -0.541   0.864
## [4,]   3 -0.598   0.844
## [5,]   4 -0.618   0.837
## [6,]   5 -0.492   0.882
## [7,]   6 -0.486   0.884
## [8,]   7 -0.455   0.895
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -2.00   0.576
## [2,]   1 -2.01   0.574
## [3,]   2 -2.04   0.562
## [4,]   3 -2.14   0.517
## [5,]   4 -2.18   0.501
## [6,]   5 -1.95   0.597
## [7,]   6 -1.94   0.602
## [8,]   7 -1.89   0.623
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01Differenced## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -40.6    0.01
## [2,]   1 -28.2    0.01
## [3,]   2 -21.9    0.01
## [4,]   3 -18.8    0.01
## [5,]   4 -19.2    0.01
## [6,]   5 -17.5    0.01
## [7,]   6 -16.5    0.01
## [8,]   7 -16.7    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -40.6    0.01
## [2,]   1 -28.2    0.01
## [3,]   2 -21.9    0.01
## [4,]   3 -18.8    0.01
## [5,]   4 -19.3    0.01
## [6,]   5 -17.6    0.01
## [7,]   6 -16.6    0.01
## [8,]   7 -16.9    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -40.6    0.01
## [2,]   1 -28.2    0.01
## [3,]   2 -21.9    0.01
## [4,]   3 -18.8    0.01
## [5,]   4 -19.3    0.01
## [6,]   5 -17.6    0.01
## [7,]   6 -16.6    0.01
## [8,]   7 -16.9    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

谷歌的增强 Dickey-Fuller 测试

Non-Differenced## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag  ADF p.value
## [1,]   0 1.21   0.940
## [2,]   1 1.15   0.933
## [3,]   2 1.23   0.943
## [4,]   3 1.16   0.935
## [5,]   4 1.20   0.939
## [6,]   5 1.28   0.950
## [7,]   6 1.36   0.956
## [8,]   7 1.34   0.954
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -1.07   0.677
## [2,]   1 -1.07   0.679
## [3,]   2 -1.06   0.680
## [4,]   3 -1.07   0.678
## [5,]   4 -1.07   0.676
## [6,]   5 -1.07   0.677
## [7,]   6 -1.04   0.687
## [8,]   7 -1.06   0.680
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -3.52  0.0404
## [2,]   1 -3.64  0.0275
## [3,]   2 -3.49  0.0430
## [4,]   3 -3.63  0.0292
## [5,]   4 -3.56  0.0364
## [6,]   5 -3.36  0.0598
## [7,]   6 -3.22  0.0851
## [8,]   7 -3.26  0.0768
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01Differenced## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -39.6    0.01
## [2,]   1 -29.5    0.01
## [3,]   2 -23.1    0.01
## [4,]   3 -20.4    0.01
## [5,]   4 -19.1    0.01
## [6,]   5 -18.0    0.01
## [7,]   6 -16.3    0.01
## [8,]   7 -16.8    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -39.6    0.01
## [2,]   1 -29.6    0.01
## [3,]   2 -23.2    0.01
## [4,]   3 -20.5    0.01
## [5,]   4 -19.2    0.01
## [6,]   5 -18.1    0.01
## [7,]   6 -16.4    0.01
## [8,]   7 -16.9    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -39.6    0.01
## [2,]   1 -29.6    0.01
## [3,]   2 -23.2    0.01
## [4,]   3 -20.5    0.01
## [5,]   4 -19.2    0.01
## [6,]   5 -18.1    0.01
## [7,]   6 -16.4    0.01
## [8,]   7 -16.9    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

格兰杰因果关系

众所周知，在观察数据中建立因果关系非常困难。格兰杰因果关系是一个较低的酒吧。它只是说，如果 X 的以前值可以预测 y 的未来值，那么 X 格兰杰原因是 y。它是通过估计 X 对 y 的滞后值的回归并执行 f 检验来执行的。如果 p 值足够小，你拒绝零假设，即 X 的所有滞后值的系数都是 0。简单地说，小 p 值表示滞后的 x 对未来的 y 有预测能力，并有相应的置信度。

下面的输出简单地说明了 FANG 的每只股票都比其他股票具有预测能力。这意味着同步建模是一个好方法。

## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F    Pr(>F)    
## 1   1574                         
## 2   1604 -30 2.9968 1.226e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F    Pr(>F)    
## 1   1574                         
## 2   1604 -30 2.4849 1.635e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F   Pr(>F)    
## 1   1574                        
## 2   1604 -30 2.5479 9.15e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F   Pr(>F)    
## 1   1574                        
## 2   1604 -30 2.5026 1.39e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df     F    Pr(>F)    
## 1   1574                        
## 2   1604 -30 3.838 2.082e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F    Pr(>F)    
## 1   1574                         
## 2   1604 -30 2.0986 0.0004837 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F   Pr(>F)    
## 1   1574                        
## 2   1604 -30 2.0299 0.000852 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F    Pr(>F)    
## 1   1574                         
## 2   1604 -30 3.3817 2.497e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F    Pr(>F)    
## 1   1574                         
## 2   1604 -30 3.0723 5.784e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F   Pr(>F)   
## 1   1574                       
## 2   1604 -30 1.7205 0.009195 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df     F   Pr(>F)   
## 1   1574                      
## 2   1604 -30 1.728 0.008712 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: j ~ Lags(j, 1:30) + Lags(i, 1:30)
## Model 2: j ~ Lags(j, 1:30)
##   Res.Df  Df      F    Pr(>F)    
## 1   1574                         
## 2   1604 -30 3.3607 3.101e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

协整检验

协整是我反复提到的一个词。这仅仅意味着我们的股票中至少有两只之间存在长期关系。协整检验被称为 Johansen 检验，以开发它的统计学家/计量经济学家的名字命名。它可以用两种方式来表述，迹检验和最大特征值检验，它们有不同的假设。细节是相当技术性的，所以 TL；它的 dr 版本是，如果长期关系矩阵的秩 r 为 0 或等于被测试的时间序列的数量，则不存在协整，不同的建模技术是合适的。如果 r 大于 0 且小于时间序列数，则存在 r 个协整关系。

从下面的输出中，我们可以非常确定存在大于 0 的协整关系，但是我们不能在任何合理的统计显著性水平上确定存在大于 1 的协整关系。因此，我们在建模时将假设 FANG 股票之间存在一种协整关系。

###################### 
# Johansen-Procedure # 
######################Test type: trace statistic , with linear trend in cointegrationEigenvalues (lambda):
[1] 2.054935e-02 1.133947e-02 6.555429e-03 4.658323e-03 -2.362375e-18Values of teststatistic and critical values of test:test 10pct 5pct 1pct
r <= 3 | 7.63 10.49 12.25 16.26
r <= 2 | 18.39 22.76 25.32 30.45
r <= 1 | 37.03 39.06 42.44 48.45
r = 0 | 70.98 59.14 62.99 70.05Eigenvectors, normalised to first column:
(These are the cointegration relations)Close.fb.l30 Close.amzn.l30 Close.nflx.l30 Close.googl.l30 trend.l30
Close.fb.l30 1.000000000 1.0000000 1.0000000 1.000000000 1.0000000
Close.amzn.l30 -0.004276658 1.0364900 0.2719707 0.102530094 -0.1327168
Close.nflx.l30 0.159547714 -1.8623141 -1.2702538 -0.746632744 0.5949550
Close.googl.l30 -0.143075556 -1.9821424 0.4921765 -0.150661845 0.2315739
trend.l30 -0.054567826 0.2916345 -0.3575875 0.005284555 -0.3248181Weights W:
(This is the loading matrix)Close.fb.l30 Close.amzn.l30 Close.nflx.l30 Close.googl.l30 trend.l30
Close.fb.d -0.00219271 -0.0006776299 -0.003233474 0.00356431 1.732653e-17
Close.amzn.d 0.18630132 -0.0084137379 -0.008556667 0.01693689 -8.853648e-16
Close.nflx.d 0.02105269 -0.0006557141 0.001122418 0.01106052 -1.271123e-16
Close.googl.d 0.10670446 0.0031632593 -0.012012380 0.01474776 -9.437844e-16

选择模型的滞后数量

为了选择包含在最终模型中的滞后数量，我选取了我的训练数据的一个小子集；从 2018 年 5 月 9 日开始到年底。这个日期在小型测试数据集中留下了 163 个观察值，这与测试数据集中的观察值数量相等。然后，我用 1 到 30 个小训练数据估计模型，然后用该模型估计小测试数据集。我计算了每只股票的平均绝对百分比误差(MAPE ),并将它们相加。然后我选择了使每只股票的 MAPE 总和最小的滞后数。最小化该调谐数据集中的总 MAPE 的滞后数量是 5。

模型有多精确？

评估模型准确性时，进行基线比较会有所帮助。对于分类，它可能是总是预测模式的误分类率。对于横截面回归，它可以是因变量的平均值。对于时间序列，我希望我的预测是前一段时间的值。这意味着我的基准准确度等级是说明天的股价就是今天的股价的准确度。

这样，脸书的基准 MAPE 为 1.36%，亚马逊为 1.20%，网飞为 1.70%，谷歌为 1.15%。对于一种不依赖复杂统计模型的预测方法来说，这还不算太糟糕。VECM 能胜过这个吗？

## [1] "Close.fb"
## [1] 1.360928
## [1] "Close.amzn"
## [1] 1.203783
## [1] "Close.nflx"
## [1] 1.697295
## [1] "Close.googl"
## [1] 1.150686

现在我们可以估计 VECM。它将具有下面的简化形式，其中α的列向量是对长期调整的速度，β的行向量包含长期关系系数。

Formulation of the VECM with L lags

是的。VECM 得出脸书的 MAPE 为 0.30%，亚马逊为 0.28%，网飞为 0.40%，谷歌为 0.29%，这还不算太差。这些 MAPE 大约是基线错误率的四分之一。

## [1] "Close.fb"
## [1] 0.2972355
## [1] "Close.amzn"
## [1] 0.277129
## [1] "Close.nflx"
## [1] 0.3978918
## [1] "Close.googl"
## [1] 0.2943474

绘制每只股票的预测收盘价和实际收盘价表明，预测股价非常接近实际股价。由于股票价格比例的差异，每只股票都有自己的图表，以便更好地查看拟合值和预测值。

关闭思路

所以这个模型表现不错。有什么意义？该模型预测收盘股价。因此，如果交易价格远低于预测的收盘价，远低于预测的收盘价是置信区间和你的风险承受能力的函数，在收盘前买入并卖出股票是有意义的；前提是没有充足的理由让股价持续下跌。如果价格远高于预测的收盘价格，卖出它是有意义的()我不是在提供财务建议。只是陈述如果我愿意，我会如何使用这个模型。)。只有当你是日内交易者时，这些行为才有意义，但是这个模型的洞察力是有限的，因为它只预测前一天。

模型的准确性告诉我们模型是适合方股票的。从这里开始，重新构建模型来预测更远的时间，以决定一个期权是否被正确定价是一个合乎逻辑的地方。记住预测越来越远的未来是困难的，这一点很重要。随着人们对未来的预测越来越远，误差带变得越来越宽。在实施交易策略之前，严格的回溯测试和过去的表现不能保证未来的表现是关键。

我希望你对协整时间序列和 VEC 建模有所了解。还有很多工作没有在这里讨论，比如解释长期关系矩阵和脉冲响应函数。我希望这个博客能激发人们对此类数据的兴趣。

我的代码和数据可以在https://github.com/jkclem/FANG-stock-prediction获得。

从赛季表现预测花样滑冰世锦赛排名

原文：https://towardsdatascience.com/predicting-figure-skating-championship-ranking-from-season-performances-fc704fa7971a?source=collection_archive---------14-----------------------

体育分析

第 1 部分:线性模型

• 要看我为这个项目写的代码，可以看看它的 Github 回购
• 对于项目的其他部分:第一部分， 第二部分第三部分第四部分第五部分第六部分****

背景

几年前我开始关注花样滑冰，现在已经成为这项运动的忠实粉丝。与许多运动类似，花样滑冰比赛分为两大类:

• 每年 10 月至次年 3 月的一系列小型赛事，包括:在 6 个不同国家(美、加、俄、法、中、日)举行的大奖赛系列赛，大奖赛系列赛 6 名最佳选手之间的大奖赛决赛比赛，欧洲运动员的欧洲锦标赛，以及其他大洲运动员的另一场比赛，适当的名称为四大洲。当然，每四年，人们熟悉的冬季奥运会会在二月举行。
• 一年一度的世界锦标赛:这是本赛季的盛大赛事，通常在三月底举行，为本赛季画上圆满的句号。全世界大约有 24 名选手将参加每个项目(男、女)的比赛。

问题

每年，在世界锦标赛之前，互联网上都充斥着对谁将登上冠军领奖台的预测。做出预测的显而易见的来源是赛季早期比赛中运动员的分数，根据这些分数对运动员进行排名的显而易见的方法是比较每个运动员在他或她之前所有比赛中的平均分数。

然而，这种方法的一个潜在问题是分数是跨不同事件平均的，没有两个事件是相同的。首先，每个项目的评判小组大不相同。另一方面，可能有一些其他因素可以系统地提高或降低某个项目的运动员分数，例如冰的情况，项目发生的时间，甚至海拔高度(运动员可能不适应高海拔的项目)。正如下面 2017 赛季男子滑冰运动员的方框图所示，每个项目的得分中心和分布可能会明显不同:

Box plot of skater scores for 2017 season events

因此，一个排名模型可以梳理出运动员效应**，每个运动员的内在能力，从事件效应，一个事件对运动员表现的影响，可能会对谁在世界锦标赛中表现更好或更差做出更好的预测。**

在这个项目中，我们将探索几个这样的模型来实现这个目标。你正在阅读的第一部分涉及简单的线性模型，而项目的后续部分将涵盖更复杂的模型。

数据

数据抓取

好消息是:从 2005 年(现行评分系统实施的那一年)开始的所有比赛的分数都可以通过国际滑冰联盟的官方网站获得。分数包含在普通 HTML 页面的表格中，所以使用 Beautiful Soup 来抓取它们就足够了。

坏消息是:每个赛事都有自己的分数网页。正如十多年来不同活动组织者建立的网页所预期的那样，它们有轻微但非常令人沮丧的不同链接模板和 HTML 格式，从 2005 年的页面和 2018 年的页面可以看出。

Left: score page in 2005. Right: score page in 2018.

因此，我们非常小心地确保所有的分数都是准确的。例如，下面是几个链接模板，用于捕捉所有网页的事件分数:

**'http://www.isuresults.com/results/gp{0}{1}/CAT00{3}RS.HTM',
'http://www.isuresults.com/results/gp{0}20{1}/CAT00{3}RS.HTM',
'http://www.isuresults.com/results/gp{0}{1}{2}/CAT00{3}RS.HTM',
'http://www.isuresults.com/results/season{1}{2}/gp{0}{1}{2}/CAT00{3}RS.HTM',
'http://www.isuresults.com/results/season{1}{2}/gp{0}{1}{2}/data0{3}90.htm',
'http://www.isuresults.com/results/season{1}{2}/gp{0}20{1}/CAT00{3}RS.HTM',
'http://www.isuresults.com/results/season{1}{2}/gp{0}20{1}/data0{3}90.htm'**

链接样板中的占位符有助于识别特定事件，包括:

• {0}:赛事类型(如'usa'代表美国大奖赛，'gpf'代表大奖赛决赛，以此类推)
• {1}:季节的开始年份(如'06')
• {2}:季节的结束年(如'07')
• {3}:对应评分页面的选手性别('1'为男性，'2'为女性)

例如，如果{0}='gpf'、{1}='06'、{2}='07'、{3}='1'，那么相应的模板就变成了http://www.isuresults.com/results/gpf0607/CAT001RS.HTM，这是 2006-2007 赛季男子滑冰大奖赛总决赛的比分网页。如果一个模板不起作用，scraper 会沿着列表向下尝试另一个模板，直到找到一个有效的链接。

数据清理

一旦所有事件的分数被清除，它们将通过以下步骤被清除:

1. 删除所有多余的空格，包括不换行空格如\xa0
2. 将溜冰者的名字从姓氏中分离出来。由于名称顺序的不一致，这是必要的:一些页面将名称存储为First Name LAST NAME，而一些存储为LAST NAME First Name。因此，将名和姓分开可以使运动员在他或她参加的所有比赛中被完全识别。值得庆幸的是，所有选手的姓氏都是大写的，只需要一些简单的正则表达式就可以进行拆分。例如，匹配运动员姓氏的正则表达式的一部分是[A-ZÄÖÜ-]+，它允许像带重音符号的PÖYKIÖ或带连字符的GILLERON-GORRY这样的姓氏被捕获。
3. 标准化备选名称拼写。由于运动员的数量足够少，通过他们的名字进行人工检查可以识别出运动员名字在多个项目中的不同拼写，例如名字的Min-Jeong和Min-Jung，或者姓氏的TUKTAMYSHEVA和TUKTAMISHEVA。从下面的列表中可以看出，可供选择的拼写大多是有多个名字(有时用连字符连接)的韩国滑冰运动员，夹杂着多音节的俄罗斯名字和充满元音的北欧名字。

Alternative spellings of skaters’ names (blue: male, pink: female)

数据分割

在 14 年的粗略分数中，随机 10 年的分数将用于训练和验证不同的排名模型，而剩余 4 年的分数将用于在最后对模型进行基准测试。这个项目的后面部分将专注于更复杂的选手排名模型，我将利用 4 年的测试时间在项目的最后部分对所有模型进行基准测试。

排名度量

在我们比较不同的模型对世界锦标赛选手的排名之前，让我们首先定义我们如何评价一个排名比另一个更好或更差。

我们知道，一个好的预测排名是一个与最终世界锦标赛结果非常接近的排名。可以衡量两个排名之间相似性的一个简单度量就是 肯德尔排名相关系数 (也叫肯德尔的τ系数)。让我们通过一个简单的例子来看看这个指标是如何工作的:

**Predicted ranking   Actual ranking  
  1\. AARON            1\. BRUCE                 
  2\. BRUCE            2\. AARON                   
  3\. CRAIG            3\. CRAIG**

上表显示了一些预测排名(左)和世界锦标赛后的实际排名(右)，每个从最好到最差，为 3 个虚构的选手。比较这两个排名，我们可以看到预测的排名错误地将 Aaron 置于 Bruce 之上。然而，克雷格的排名是正确的:他在两个排名中都低于其他两名选手。

肯德尔τ系数(𝜏)可以用下面的公式量化上述描述:

n: number of ranked skaters

该公式的工作原理如下:

• 根据AARON > BRUCE > CRAIG的预测排名，可以生成有序的选手对，其中选手对中第一个选手的排名高于第二个选手。对于 3 个溜冰者(n = 3)的排名，我们最终得到 3×(3–1)/2 = 3 个有序对。它们是:

**(AARON, BRUCE), (AARON, CRAIG), (BRUCE, CRAIG)**

• 在这 3 对有序的组合中，第一对(AARON, BRUCE)不再保持BRUCE > AARON > CRAIG的冠军排名，而后面两对仍然保持。换句话说，一致对——在两种排序中都成立的有序对——的数量是 2，而不一致对**——在一种排序中成立但在另一种排序中不成立的有序对——的数量是 1。**
• 因此，这两个排名的肯德尔τ系数为:

结果，在两个排序中存在的有序对越多，这两个排序之间的肯德尔τ系数就越高。两个相同的排名将有𝜏 =1，而两个完全相反的排名(想想:一个排名从最高到最低，而另一个从最低到最高)将有𝜏 = -1。因此，当与同赛季的世界冠军排名进行比较时，肯德尔的τ系数较高(更接近 1)的预测排名将优于系数较低的预测排名。

履行

在 Python 中，scipy库内置了[scipy.stats.kendalltau](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kendalltau.html)函数来计算系数。然而，“从零开始”实现这种计算非常容易:itertools模块中的[combinations](https://docs.python.org/2/library/itertools.html#itertools.combinations)函数将从每个排名中生成有序对，两个排名之间的一致对只不过是两组有序对的交集(通过&运算符):

**from itertools import combinations
ranking1 = ['AARON', 'BRUCE', 'CRAIG']
ranking2 = ['BRUCE', 'AARON', 'CRAIG']# Generate set of ordered pairs from each ranking
ordered_pairs1 = set(combinations(ranking1, 2))
# {('AARON', 'BRUCE'), ('AARON', 'CRAIG'), ('BRUCE', 'CRAIG')}
ordered_pairs2 = set(combinations(ranking2, 2))
# {('AARON', 'CRAIG'), ('BRUCE', 'AARON'), ('BRUCE', 'CRAIG')}# Calculate number of total, concordant & discordant pairs and tau
n_pairs = len(ordered_pairs1) # 3
n_concordant_pairs = len(ordered_pairs1 & ordered_pairs2) # 2
n_discordant_pairs = n_pairs - n_concordant_pairs # 1# Calculate Kendall's tau from pair counts
tau = (n_concordant_pairs - n_discordant_pairs) / n_pairs # 0.33**

季节平均模型(基线)

正如开始提到的，这是预测世界锦标赛上运动员排名的最简单和最明显的方法，可以通过以下方式实现:

1. 平均每位选手在本赛季参加的所有比赛中的得分
2. 根据这些赛季平均水平，从最高到最低对运动员进行排名，并使用该排名作为世界锦标赛的预测排名

以 2017 年男子滑冰运动员的成绩为例，我们可以看到，使用赛季平均值预测的排名与世界锦标赛的最终排名具有相对较强的相关性(见下面的排行榜对比)。和任何运动一样，都有很大的惊喜:Denis, TEN排名下降了 10 位，从第 6 位降至第 16 位，而Brendan, KERRY上升了 7 位，从第 22 位升至第 15 位。然而，快速浏览一下，赛季中顶级选手的表现往往比低级选手更稳定，这可以从排行榜顶部与底部相比相对较少的排名中看出。

Ranking by season average vs. world ranking (for male skaters in the 2017 season). Colors are based on the world championship ranking.

24 个滑手，2017 年世锦赛实际排名有 24×23/2 = 276 个可能的有序对。在这 276 对中，季节平均模型正确预测了 234 对(如Yuzuru, HANYU > Javier, FERNANDEZ)，剩下 42 对预测错误(如Boyang, JIN > Nathan, CHEN))。这转化为肯德尔τ分数为(234 - 42)/276 = 0.695。对于基线模型来说还不算太差！

加性模型

理论

让我们从季节平均值的基线模型转移到其他模型，这些模型有望将事件影响与运动员影响分开。一个简单的这样的模型是假设在给定的事件中给定的溜冰者的分数是以下的加**😗*

• 赛季中所有分数的基准分数
• 潜在的事件分数对于该事件是唯一的，但是对于该事件的所有运动员是相同的
• 一个潜在的滑冰者分数对于该滑冰者来说是唯一的，但是对于该滑冰者参加的所有项目来说是相同的
• 零均值高斯噪声

这个加法模型可以用数学方法表示为:

• y_event-skater:给定运动员在给定比赛项目中的得分
• θ_baseline:基线得分(整个赛季得分不变)
• θ_event:潜在事件分数(运动员之间不变)
• θ_skater:潜在选手得分(跨赛事不变)
• N(0, σ²):均值为零的高斯噪声

乍一看，似乎没有简单的方法来学习这些潜在分数，因为每对运动员都有不同的潜在分数组合。换句话说，不存在机器学习问题所预期的所有数据点的固定特征集。

但是，模型公式可以重写为以下形式:

这种形式的模型公式中的新术语是:

• θ_e:事件的潜在得分e
• I(e = event):如果事件e是给定事件，则为 1，否则为 0
• θ_s:滑手潜在得分s
• I(s = skater) : 1 如果滑冰者s是给定的滑冰者，否则为 0

以这种方式书写，很容易看出，对于给定的一对运动员，只有该运动员和该项目的潜在得分将计入该对运动员的实际得分，这与前面的公式完全相同。这是因为在上述公式中，那些潜在得分的系数I(e = event)和I(s = skater)将为 1，而其他滑手和项目的潜在得分的系数将为 0。

同样，以这种方式书写，每个运动员的分数可以简单地认为是基线分数和零均值高斯噪声之上的二元变量I(s = skater)和I(e = event)的线性组合。换句话说，问题就变成了一个简单的二元特征线性回归，其回归系数的确是所有项目和滑手的潜在得分——分别为θ_e和θ_s——回归截距不过是基线得分θ_baseline。

编码

我们通过一个小例子(使用 2017 赛季的男性得分)来看看这个线性回归模型是如何编码的:

• 清理完数据后，我们以下面的数据框season_scores结束，它代表了本赛季所有赛事选手的分数(当然，除了世界锦标赛)。对于 2017 赛季，有来自 9 项赛事的 62 名选手的 120 个分数。

• 接下来，我们使用 panda 的[get_dummies](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.get_dummies.html)函数将这个数据帧转换成包含虚拟特征的数据帧dummies，虚拟特征就是上面提到的二进制特征。

******dummies = pd.get_dummies(season_scores[['name', 'event']], prefix=['', ''], prefix_sep='', drop_first=True)******

请注意，我们必须使用函数的drop_first=True参数，因为包含所有溜冰者和事件的二元变量将导致二元特征完全共线，这对于线性回归是不允许的；这有时被称为虚拟变量陷阱。实际上，这种争论将会抛弃一个运动员和一个项目的二元特征。由于有 62 名溜冰者和 9 个事件，我们的dummies数据框将有 61 个溜冰者二进制列和 8 个事件二进制列，总共有 69 个二进制要素。可以通过下面的方法找到被丢弃的运动员和事件的身份:**

******# Get unique names of skaters and events from original data frame
unique_skaters = season_scores['name'].unique()
unique_events = season_scores['event'].unique()# Get unique names of skaters and events from dummies data frame 
dummies_skater_count = len(unique_skaters) - 1
dummies_skaters = dummies.columns[:dummies_skater_count]
dummies_events = dummies.columns[dummies_skater_count:]# Dropped skater is 'Adam, RIPPON'
dropped_skater = list(set(unique_skaters) - set(dummies_skaters))[0]# Dropped event is '4C'
dropped_event = list(set(unique_events) - set(dummies_events))[0]******

删除这两列后，得到的dummies数据表如下左图所示，其中 120 行对应于该赛季的 120 个唯一得分，69 列为二进制哑变量。所有 1 的列被附加到数据表的左侧，以对应基线得分，该基线得分出现在所有赛季得分中。因此，基线分数将只是回归系数中的一个，它现在有 70 个。

Left: feature matrix (predictors). Center: coefficient vector. Right: score vector (response). Highlighted in red is how one row of the feature matrix can make a dot product with the coefficient vector to approximate the season score of an event-skater pair (CA-Yuzuru, HANYU).

创建这个 120×70 数据表的原因是，它可以表示一个二进制特征矩阵X，当乘以表示所有 70 个回归系数(包括基线/截距)的向量θ时，它可以近似该赛季的所有得分——长度为 120 的向量y。换句话说，线性回归问题可以用矩阵形式表示为:

其中ε表示对每个赛季得分有贡献的高斯噪声。因此，回归系数的最小二乘估计可以通过求解线性回归的标准方程来轻松计算:

由于每个季节的数据都很小，如微小的 120×70 要素矩阵所示，因此在 Python 中可以很容易地求解法线方程，以返回回归系数的精确值(与梯度下降之类的东西相反)。这个的代码如下，用.T做矩阵转置，@做矩阵乘法，np.linalg.inv做矩阵求逆。

******# Create numpy feature matrix and response vector
X = dummies.values
y = season_scores['score'].values# Append column of all 1's to the left (0th-index) of feature matrix
X = np.insert(X, obj=0, values=1, axis=1)# Apply normal equation to find linear regression coefficients (including intercept)
coefs = np.linalg.inv(X.T @ X) @ (X.T @ y)******

从通过解正规方程找到的回归系数中，可以提取基线分数、溜冰者分数和项目分数。请注意，我们需要添加回运动员和事件的分数，这些分数在回归之前被删除；他们的分数为零，因为通过放弃他们，他们的分数已经被吸收到基线分数中。

******# Get baseline, skater & event scores from regression coefficients
baseline_score = coefs[0]
skater_scores = pd.Series(coefs[1:dummies_skater_count+1], index=dummies_skaters)
event_scores = pd.Series(coefs[dummies_skater_count+1:], index=dummies_events)# Add back scores of dropped skater and dropped event
skater_scores[dropped_skater] = 0
event_scores[dropped_event] = 0# Sort skater scores and event scores from highest to lowest
skater_scores.sort_values(ascending=False, inplace=True)
event_scores.sort_values(ascending=False, inplace=True)******

结果

一旦对潜在的溜冰者分数进行了排序，基于加法模型的预测排名就像将具有最高潜在分数的溜冰者排名为第一，将具有第二高潜在分数的溜冰者排名为第二一样简单，以此类推。

此外，这些潜在分数可以被解释为一个运动员与被淘汰的运动员相比好或差多少，一旦去除了事件影响(根据加法模型)，参考分数为零。

另外，请注意，在已经计算出潜在分数的 62 名滑手中，只有 24 人获得了 2017 年世界锦标赛的参赛资格。因此，我们只根据那些滑手的预测排名来评估模型。

估价

首先，让我们看看加法模型在逼近赛季中的实际得分时表现如何，因为这是回归的主要目标。下面是一个图表，显示了本赛季每对赛事选手的预测分数，每个分数都按其潜在成分(基线分数+赛事分数+选手分数)与该对选手的真实分数进行了细分。

Each tick represents a skater. Each gray bar represents the discrepancy between predicted and true season scores.

上图显示，预测得分与赛季中的真实得分非常接近，与真实得分的最大偏差(也称为残差)不到 30 分。我们可以根据预测得分绘制这些残差，如果使用了季节平均模型，还可以绘制相同的残差。此外，我们可以绘制这些残差的直方图和 QQ 图，以确认线性回归模型的高斯噪声假设是否合理。

Top: residual vs. predicted score for season average and additive models. Bottom: histogram and QQ plot of additive models’s residuals

从上面的情节可以看出:

• 来自加法模型的残差在整个预测分数范围内是相当随机的。这加强了随机噪声的模型假设，与不同分数范围的噪声不同(也称为异方差)，其在残差与预测分数图中具有特征性的“扇形”形状。
• 从残差的直方图和 QQ 图来看，残差的分布非常接近高斯分布。这证明了线性模型的零均值高斯噪声假设。
• 然而，似乎有异常大量的零残差，从直方图的高峰以及 QQ 图中间的扭结可以看出。经过一些检查，这些零残差原来只属于在赛季中只参加一个项目的运动员。这是有意义的，例如，如果梯度下降被用来寻找潜在得分:这样的滑冰运动员的潜在得分只对单个项目的预测得分有贡献，而整个赛季的其他预测得分不变。因此，在梯度下降过程中，潜在的运动员分数可以自由地自我调整，使得该项目-运动员对的真实分数和预测分数之间的差，即残差变为零。
• 最后，将加法模型的残差与平均模型的残差进行比较，结果显示它们在规模上非常相似。事实上，线性模型的均方根误差 (RMSE)为 8.84 点，比季节平均模型的 10.3 点要好一些。

其次，也是更重要的，让我们看看加法模型在预测世界锦标赛的最终排名方面做得如何(与运动员赛季平均水平的基线模型相比)。细心的读者会注意到，之前显示的加法模型得出的排名与赛季平均排名非常相似，这在将这两个排名相互比较时得到了证实，同时也与世界锦标赛排名进行了比较:

Left: predicted ranking from additive (linear) model. Center: predicted ranking from season average model. Right: actual ranking from 2017 world championship. Colors are based on the world championship ranking.

在比较这两种模型的预测排名时，有一些有趣的情况(以红色突出显示):

• 还有一些滑手按季均排名更准确，比如Moris, KVITELASHVILI(加法排名第 16，季均排名第 12，世界排名第 13)。
• 相比之下，也有一些从加性模型预测的排名更准确，比如Paul, FENTZ(加性第 20，季均第 17，世界第 20)。
• 当然，由于体育运动的不可预测性，有些运动员的两个模型都没有提供好的预测:Boyang, JIN被预测会在两个排名中名列第七，但却出人意料地在世界锦标赛中获得了第三名。

总的来说，在世界锦标赛的 276 个有序对中，有 239 个是加法模型正确得到的，这转化为肯德尔的 tau 分数为 0.732。与具有 234 个正确对的季节平均模型(肯德尔的τ= 0.695)相比，相加模型与基线模型相比似乎提供很少的预测改进。因此，让我们转向下一款有望表现更好的车型。

乘法模型

理论

我们可以将这些潜在得分乘以，再乘以零均值高斯噪声的指数，从而得出赛季得分，而不是将赛季中的项目-运动员得分视为由一些潜在得分相加而成(基线+项目+运动员)。这种乘法关系可以用数学方法表示为:****

乘法模型相对于加法模型的一个潜在优势是，潜在分数可以在该模型中相互作用。例如，在附加模型中，一个给定的溜冰者的潜在效应保持不变，不管他或她参加什么项目。相比之下，乘法模型中的溜冰者分数充当基线和事件分数之上的“乘数”:如果你愿意，可以把这想象为溜冰者在“利用”每个事件方面有多好。

我们在公式末尾乘以高斯噪声的指数的原因是因为可以在整个公式中取对数，这导致:

换句话说，与加法模型一样，乘法模型只不过是一个具有二元特征的线性回归问题。然而，不同之处在于，响应现在是赛季得分的对数，回归系数现在是潜在得分的对数。

编码

由于模型核心的线性回归保持不变，乘法模型的编码很大程度上与加法模型相同，除了在最后解正规方程时:我们用X.T @ np.log(y)代替了X.T @ y。

******coefs = np.linalg.inv(X.T @ X) @ (X.T @ np.log(y))******

另一个小的修改发生在回归系数被获得之后，我们需要取它们的指数来得到潜在的基线、事件和选手分数。关于如何从回归系数中提取baseline_score、event_scores和skater_scores，请参考加性模型。

******baseline_score = np.exp(baseline_score)
event_scores = np.exp(event_scores)
skater_scores = np.exp(skater_scores)******

结果

以下是乘法模型恢复的选手潜在得分:

最优秀的滑手，Yuzuru, HANYU“乘数”为 1.15；这意味着对于任何一个项目，在基线得分和项目得分相乘后，他的预测得分将是该数字的 1.15 倍。垫底的是排名最差的滑手Julian Zhi Jie, YEE，乘数为 0.77。

估价

尽管它们的形式不同，乘法模型为赛季中的每一对运动员提供了与加法模型非常相似的预测。当绘制两个模型的真实得分和预测得分之间的偏差以及赛季平均值的基线模型时，可以在下面看到这一点。因此，两个模型的 RMSE 实际上是不可区分的:加法模型为 8.84，乘法模型为 8.88(而季节平均模型为 10.3)。

Deviation from true score for each event-skater pair of the season for all 3 models (average, additive, multiplicative)

在预测世界锦标赛的最终排名方面，这两个模型再次几乎相同(参见它们与最终世界排名的比较):

Predicted rankings of multiplicative and additive models compared to the 2017 world ranking. Colors are based on the world championship ranking.

2017 赛季的两个排名只有两个不同，在上图中用红色突出显示:

• 乘法排名预测Nathan, CHEN排名高于Javier, FERNANDEZ，而加法排名预测相反。事实证明，在最终的世界锦标赛中，加法排名是正确的。
• 乘法排名预测Kevin, REYNOLDS排名高于Maxim, KOVTUN，而加法排名预测相反。然而，在这种情况下，乘法排名在最终的世界锦标赛中确实是正确的。

因此，乘法排序与世界排序具有完全相同的一致对数量(276 对中的 239 对)和与加法排序相同的肯德尔τ(0.732)，这仍然高于季节平均值基线模型的肯德尔τ(0.695)。换句话说，至少在 2017 赛季，两个模型在预测世界冠军的最终排名方面似乎没有什么表现差异。

典型惩罚

假设模型选择(加法或乘法)似乎不会显著影响预测性能，让我们看看如何以另一种方式改进任一模型。

具有讽刺意味的是，改进这些模型的一个常用方法是惩罚它们(也称为规则化)。对于这个问题，罚分可以理解为人为地将潜在事件和滑手分数向零收缩。一方面，这种萎缩无疑会使模型在预测赛季得分方面不太准确。另一方面，它可以防止模型过度适应赛季分数，赛季分数可能会很嘈杂，可能不会反映溜冰者的真实能力(可能他们在那天早些时候食物中毒)。因此，这可能会更好地捕捉运动员的真实能力，并对他们进行更准确的排名。****

对于这个问题，我会选择 岭回归 来惩罚两个模型，原因有二:

• 岭回归及其近亲(例如套索和弹性网回归)通常用于正则化线性回归问题。它也非常容易实现并合并到目前使用的正规方程中，如下所示。
• 但是，与类似的正则化类型(如 lasso)相比，岭回归的一个优点是回归系数永远不会完全收缩为零。换句话说，多个滑手之间永远不会有平手，都是零潜分，适合这个排名问题。

理论

让我们看看如何将岭回归纳入用于线性回归的正规方程。回想一下加法模型的正规方程的解:

• θ:回归系数(70: 1 基线+ 8 项+ 61 名选手)
• X:二元特征表(120×70)
• y:响应向量(120)。对于乘法模型，这变成了log(y)

我们可以修改岭回归的这个解决方案，使回归系数向零收缩(下面公式的推导见最后的参考资料部分):

这个等式中的额外项是:

• λ:惩罚参数(标量)
• I:单位矩阵(70×70 矩阵，其条目除了对角线上的 1 之外都是 0，每个对角线条目对应一个回归系数)。注意，通常情况下，我们不会惩罚回归的截距。因此，单位矩阵I的第一个对角元素将是 0 而不是 1。

惩罚参数λ越高，模型将变得越惩罚，即潜在事件和运动员分数进一步向零收缩。相比之下，当λ = 0 时，模型成为原始的、未受惩罚的版本。

编码

给定上面的公式，为加性模型的岭回归编写正规方程的解再简单不过了，因为我们已经为同一个方程的正规版本编写了代码。

******# Create 70x70 identity matrix
I = np.identity(n=X.shape[1])# Set first diagonal entry of identity matrix to zero
I[0, 0] = 0# Solving normal equation for ridge regression
lambda_reg = 0.01
coefs = np.linalg.inv(X.T @ X + lambda_reg * I) @ (X.T @ y)******

结果

以 2017 赛季的加性模型为例，当λ增加 10 倍时，潜在项目和滑手分数一般会向零收缩。相比之下，基线分数相对不受惩罚参数的影响(这本来就不是它的本意):

Effect of λ (penalty parameter) on baseline, skater, and event scores of additive model in 2017 season

估价

当潜在事件和运动员分数向零收缩时，模型在逼近赛季分数时自然会变差。因此，模型 RMSE 将增加(见下图):在极端情况下，当λ = 10000 时，RMSE 为 35.6，这比仅仅取每个溜冰者的季节平均值(RMSE = 10.3)要差得多，更不用说原始的、未被惩罚的模型(RMSE = 8.8)。

Effect of λ (penalty parameter) on RMSE, Kendall’s tau, and predicted ranking of additive model in 2017 season. Colors are based on world championship ranking.

然而，这种极端的萎缩会导致更好地预测世界锦标赛的排名吗？唉，没有。从上图来看，在所有 276 个有序配对中，极度惩罚模型只正确获得了 231 对，低于未惩罚模型的 239 对，甚至低于季节平均模型的 234 对！通过绘制当λ增加时预测排名如何变化的图表(见上文)，我们可以看到，直到λ = 0.1，惩罚模型的预测排名与原始模型相比没有太大变化，并且一致对的数量保持在 239。然而，在此之后，惩罚模型的排名开始有点混乱，这与和谐对数量的下降相一致(从λ = 0.1 时的 239 到λ = 1 时的 235)。

简而言之，惩罚加性模型造成了双重打击:它不仅使模型在逼近赛季得分方面变得更差，而且在预测世界冠军排名方面也变得更差。

然而，这只是单个赛季(2017 年)的结果。让我们在训练样本中对所有 10 年应用相同的模型(加法和乘法)。对于每个模型，我们从季节平均值的基线模型计算它们在 RMSE 和肯德尔τ中的差异。然后，我们使用 9 个自由度的 t 分布(因为 n=10 年)绘制这些差异的平均值，以及它们的 95%置信区间。这与执行配对 t-检验来比较每个模型与基线相比的 RMSE 和肯德尔τ的平均差异是一样的。

Effect of λ (penalty parameter) on differences (Δ) in RMSE and Kendall’s tau from baseline model, averaged over the 10-year training sample

从上面的图表中，我们可以看到:

• 与季节平均值的基线模型相比，加法和乘法模型实际上以相同的量影响 RMSE，并且平均以相似的量影响肯德尔τ。
• 就 RMSE(左图)而言，两个模型将 RMSE 从基线降低了λ < 0.1. However, this in itself is of little value since approximating season scores is not our goal for this project. Furthermore, the RMSE naturally increases as λ increases, especially for λ > 0.1 的统计显著量，而 RMSE 从季节平均值的基线模型显著上升。
• 就肯德尔τ(右图)而言，模型受到的惩罚越多(λ越大)，与基线相比，它们的肯德尔τ平均越差，不幸的是，这与我们之前从 2017 赛季得到的结果一致。更糟糕的是，从λ = 1 开始，这些模型开始从基线模型减少显著数量的肯德尔τ:注意肯德尔τ从基线的差异的 95%置信区间如何下降到水平线以下的零。因此，我们应该为这个项目选择非惩罚模型。

然而，即使使用未加惩罚的加法和乘法模型(λ = 0)，上图显示，平均而言，与季节平均值的基线模型相比，这些模型实际上降低了肯德尔τ😱当然，肯德尔τ差异的 95%置信区间表明，这种与基线的差异在统计学上不显著，因为它们包含了 0 处的水平线。但是，让我们分别绘制这三个模型每年的表现，以便我们可以更详细地比较这三个模型:

Each gray line represents a year in the training set, with 2017 highlighted in red

• 就 RMSE(左)而言，与训练样本中每年的季节平均模型相比，加法和乘法模型确实降低了 RMSE。然而，与 2017 赛季一样，这两个模型在 RMSE 方面没有太大差异，即它们同样接近赛季得分。
• 然而，就肯德尔的τ(右)而言，画面就不那么清晰了。正如 2017 赛季早些时候所示，加法和乘法模型预测最终的世界冠军排名都优于赛季平均模型。然而，有相当多的年份，这两种模型给出的预测都比季节平均模型差。平均所有 10 年的训练集，这并不奇怪，他们的综合表现与基线相比是如此黯淡。
• 请注意，到目前为止，所有的结果都来自男性选手。对于女滑手来说，情况基本保持不变，从下面的对比图可以看出。

Comparison between three models for female skaters

以下是关于两个模型(加法和乘法)的 RMSE 和肯德尔τ性能的最终报告，包括平均值、与季节平均值基线模型的平均差异，以及 10 年训练样本中这些差异的 95%置信区间:

对于女性滑冰运动员，我们可以看到肯德尔τ(与季节平均值的基线模型相比)差异的 95%置信区间包含加法和乘法模型的零。换句话说，在从基线预测世界冠军排名时，两个模型都没有提供任何统计上的显著改进。

结论

尽管 2017 赛季的结果最初很乐观，但事实证明，为了预测最终的世界冠军排名，该项目期间提出的两个模型(加法和乘法)并不比简单地根据赛季平均水平对滑手进行排名更好。因此，根据奥卡姆剃刀，没有令人信服的理由选择这些更复杂的模型而不是更简单的季节平均模型。

一个问题仍然存在:考虑到这些模型中潜在的事件效应和运动员效应之间看似合理的相互作用，为什么它们在预测世界锦标赛排名时没有显示出显著的改进？

一个答案是体育运动太不可预测了。当然，这些模型可以很好地预测过去发生的事件的分数，正如他们的小 RMSE 事件所证明的那样，但是很难说在未来的事件中，没有人会受伤，或者有时差反应，或者只是在冰上撞了一个怪胎。这种不可预测性可以通过所有三种模型每年高度可变的肯德尔τ来反映:对于男性滑冰运动员，这可以在 0.5 到 0.8 以上的范围内(或者女性为 0.4 到 0.7)。换句话说，仅仅因为某年的世界锦标赛恰好有利于某个模特，并不意味着下个赛季不会有任何时髦的事情发生。****

Ice, as they say, is slippery

资源

尽管与基线模型相比没有获得显著的改进，但这个项目让我知道，即使是最简单的模型，如线性回归，也可以以创造性的方式使用，如找到潜在的因素来解决这种情况下的排名问题。以下是一些有助于理解本文理论的资源:

• 对于法方程解的推导(对于非惩罚和惩罚模型)，麻省理工学院开放式课件上的机器学习课程的这个讲义非常详细。注释还解释了惩罚参数λ的重要性:在贝叶斯设置中，较高的λ意味着潜在得分的高斯先验更集中于零。因此，这些分数的估计值也更有可能接近于零，正如在这个项目中所看到的那样。
• 对于乘法模型，似乎可以用在时序分析中。此外，从这个讲座幻灯片来看，它似乎也可以被称为对数线性模型。然而，这些模型似乎包括取预测值的对数，而不是取感兴趣系数的对数(这正是本项目所做的)。

从赛季表现预测花样滑冰世锦赛排名

原文：https://towardsdatascience.com/predicting-figure-skating-world-championship-ranking-from-season-performances-part-2-hybrid-7d296747b15?source=collection_archive---------32-----------------------

体育分析

第 2 部分:混合模型

• 要看我为这个项目写的代码，可以看看它的 Github 回购
• 对于项目的其他部分: 第一部分 、第二部分 第三部分 、 第四部分 、 第五部分 、 第六部分

背景

在项目的第一部分，我试图根据运动员在该赛季以前的比赛项目中获得的分数来预测年度花样滑冰世界锦标赛中的排名。主要策略是将滑手效应(每个滑手的内在能力)与事件效应(一个事件对滑手表现的影响)分开，以便建立更准确的排名。

针对这个排名问题提出了两个模型，其中一对运动员的赛季得分可以近似为:

• 一个基线得分，它在所有赛季得分中保持不变
• 潜在的事件分数，该分数在所有参加该事件的选手中是不变的
• 一个潜在的滑冰者分数，该分数在该滑冰者参加的所有比赛中都是不变的

这些潜在得分可以相加在一起(加法模型)，或者相乘(乘法模型)以接近赛季得分:

在任一情况下，通过重写上述公式，可以容易地找到潜在分数，使得它们成为简单的线性模型，其回归系数是潜在分数本身。因此，这些分数可以简单地通过求解线性回归的正规方程来得到。从那里，每个模型的预测排名只是潜在滑手分数的排名，从最高到最低。

问题

将这两个模型预测的排名与根据赛季平均水平简单排名的基线模型进行比较。用于评估这些模型的度量是相对于实际世界锦标赛排名的肯德尔排名相关系数(也称为肯德尔的τ):如果模型的预测排名具有更接近 1 的肯德尔τ，则模型更好，这意味着预测排名与当年世界锦标赛中的实际排名更相关。

Baseline: season average model. The 95%-confidence intervals are built using a t-distribution with 9 degrees of freedom (since n=10 years).

从以上关于这两个模型在训练集中选择的 10 年(共 14 年)的报告来看，这两个模型比赛季平均值的基线模型更好地逼近赛季得分，因为它们与基线模型之间的均方根误差 (RMSE)差异的 95%置信区间小于零。

然而，在更重要的 Kendallτ度量中，遗憾的是，与基线模型相比，这些模型没有提供该度量的任何显著增加，因为对于男性和女性溜冰者，它们与基线之间的 Kendallτ差异的 95%置信区间包含零。

混合模型

鉴于这两个模型在排名预测方面表现平平，我提出了一个新模型，它是这两个模型的奇怪混合:

• 给定事件中给定选手的分数
• θ_baseline:基线得分(整个赛季得分不变)
• θ_event:该项目潜在得分(运动员之间不变)
• θ_skater:溜冰者潜在得分(跨事件恒定)

根据模型公式，给定溜冰者在赛季期间的给定事件中获得的分数可以近似为该溜冰者的潜在分数乘以该事件的潜在分数，类似于乘法模型。但是，该产品随后会添加到基线得分中，类似于加法模型。

由于有些年份加法模型比乘法模型具有更高的肯德尔τ，而有些年份则相反，因此混合模型可能提供比单独的任一模型更好的预测排名。你可以在这个项目的 Github repo 里找到这个 Jupyter 笔记本里关于混动车型的分析。

寻找潜在分数

由于混合模型中加法和乘法的奇怪组合，我们既不能将其重新表述为线性模型，也不能取公式的对数将其转换为线性模型(类似于我们对乘法模型所做的)。

尽管如此，寻找该模型的最佳参数的策略仍然与前两个模型相同:最小化赛季中所有事件-滑冰运动员对的预测分数和真实分数之间的平方差之和:

• J:模型的目标函数；继续以 2017 赛季的男选手为例，这将是 1 个基线分数、9 个潜在事件分数(以蓝色突出显示)和 62 个潜在选手分数(以黄色突出显示)的函数。
• ŷ_e,s:来自混合模型的项目-滑冰运动员对的预测得分
• y_e,s:在赛季中记录事件-溜冰者对的真实分数
• θ_baseline:基线得分
• θ_e:事件潜在得分e
• θ_s:滑手潜在得分s

梯度下降算法

有了上面的目标函数，我们可以使用 good ol 'gradient descent来找到模型的最佳基线、事件和选手分数。让我们看看每个参数的梯度是什么样的，首先是基准分数**(θ_baseline)😗*

根据上面的公式，目标函数相对于基线得分的梯度简单地是赛季中所有得分的真实得分和预测得分之间的差的总和(从这里开始将被称为残差)。

因此，给定基线、项目和运动员分数的一些当前值，我们总是可以计算基线分数的梯度，并根据该梯度更新基线分数(因为我们希望最小化目标函数):

与梯度相乘的额外项α被称为学习速率，并控制梯度下降算法运行的速度。低α可能会使算法收敛太慢，而高α可能会使更新太大，参数不断从一个极端摆动到另一个极端，这在我这个项目中发生过几次。

给定事件的潜在得分可以类似地更新，以俄罗斯的事件(潜在得分θ_RU)为例:****

当将目标函数与该事件的潜在得分进行微分时，所有不涉及该事件的平方差都将消失，这就只剩下涉及来自该事件的参加滑手的项，包括哈维尔、费尔南德斯(θ_FERNANDEZ)和亚历山大、马约罗夫(θ_MAJOROV)。使用链式法则，我们可以看到，梯度只不过是参与该事件的每个溜冰者的残差乘以该溜冰者各自的潜在得分，然后将这些乘积加在一起。

使用这个梯度，这个俄罗斯事件的潜在得分可以使用梯度下降来更新:

最后，对于溜冰者的潜在得分，例如 MURA 的高仁(θ_MURA)，目标函数相对于该潜在得分的梯度得分为:****

2017 赛季，他参加了 2 场比赛，一场在加拿大，一场在法国，分别有潜在分数θ_CA和θ_FR。因此，梯度只不过是每个事件的残差乘以该事件的潜在得分，然后将这些乘积相加。因此，可以使用梯度下降来更新该溜冰者的潜在得分:

简而言之，寻找基线、项目和运动员分数的梯度下降算法可以概括为以下步骤:

在第二步中，可以通过检查模型的 RMSE(即平均残差平方的平方根)是否已经稳定到最小值来监控收敛。梯度下降收敛后得到的运动员分数将被用于世界锦标赛的运动员排名。

混合模型编码

给定上面的梯度下降算法，让我们看一个如何用 Python 编码的小例子。下面，我们有一个数据框(season_scores)，包含 4 名滑冰运动员(马约罗夫、费尔南德斯、葛、MURA)和 3 项赛事(加拿大、法国、俄罗斯)的 7 个赛季成绩:

首先，我们使用 pandas 的pivot_table函数将这个长表格式转换成数据透视表格式，将运动员姓名作为行索引，将事件名称作为列。将赛季得分转换成这种数据透视表格式的原因很快就会明了。

season_pivot = pd.pivot_table(season_scores[['name', 'event', 'score']], values='score', index='name', columns='event')

下面是我们最终将得到的 4×3 的赛季得分数据透视表。注意不是所有的滑手都参加所有的项目，所以数据透视表中有缺失值，用NaN(不是数字)表示。

最后，在我们开始梯度下降之前，让我们把熊猫数据帧season_pivot转换成一个 4×3 的 numpy 矩阵true_scores，这样我们可以更容易地操作它。这实际上将删除所有的行和列名称，所以让我们存储运动员名称(行)和事件名称(列)，这样我们在运行梯度下降后仍然可以识别它们。

skater_names = list(season_pivot.index) 
# ['Alexander, MAJOROV', 'Javier, FERNANDEZ', 'Misha, GE', 'Takahito, MURA']event_names = list(season_pivot.columns) 
# ['CA', 'FR', 'RU']true_scores = season_pivot.values
# array([[   nan,    nan, 192.14],
#        [   nan, 285.38, 292.98],
#        [226.07, 229.06,    nan],
#        [222.13, 248.42,    nan]])

第一步:初始化基线、潜在项目和潜在选手得分

首先，让我们使用 numpy 的带有特定种子的[RandomState](https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.random.RandomState.html)对象初始化梯度下降的基线分数(这样每次运行代码时结果都是一致的)，并调用该对象的[random_sample](https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.random.RandomState.random_sample.html#numpy.random.RandomState.random_sample)方法返回一个介于 0 和 1 之间的随机数。初始值没有理由在这个范围内。我只是觉得这样会更好看。

random_state = np.random.RandomState(seed=42)
baseline = random_state.random_sample() # 0.866

现在，让我们初始化潜在事件和选手分数。我不再逐一初始化它们，而是将所有选手的分数初始化为大小为 4×1 ( skater_scores)的列向量，将所有事件的分数初始化为大小为 1×3 ( event_scores)的行向量，其条目也是 0 到 1 之间的随机数，再次使用random_sample方法。

skater_scores = random_state.random_sample(4, 1))
# array([[0.37454012],
#        [0.95071431],
#        [0.73199394],
#        [0.59865848]])event_scores = random_state.random_sample((1, 3))
# array([[0.15601864, 0.15599452, 0.05808361]])

以向量形式初始化这些潜在得分的原因是，通过将这些向量相乘(使用矩阵乘法运算符@)并在此基础上添加基线，可以一次性计算出每对运动员的预测得分:

predicted_scores = skater_scores @ event_scores + baseline
# array([[0.52284634, 0.42976104, 1.19802617],
#        [0.48872716, 0.41705697, 1.00857581],
#        [0.46792756, 0.40931237, 0.89308382],
#        [0.39887817, 0.38360225, 0.50967974]])

为了更清楚地说明这一点，下图显示了两个潜在向量(大小分别为 4×1 和 1×3)与添加在顶部的基线相乘，以形成每个项目-选手对的预测得分的 4×3 矩阵。用红色突出显示的是一对运动员(CA-MAJOROV ),这是通过乘以潜在向量中的相应条目计算出来的，添加的基线在顶部。

Dotted line: implied indices of skaters and events, even though they don’t exist in the numpy array

接下来，由于已经计算了每对运动员的真实分数(true_scores)和预测分数(predicted_scores),我们可以很容易地计算出赛季中每个分数的残差。

residuals = predicted_scores - true_scores

回想一下，我们的true_scores numpy 矩阵包含赛季中不存在的赛事-运动员配对的NaN值。因此，当计算残差时，这些对的相应残差也是NaN。

步骤 2a:计算基线、项目和运动员分数的梯度

残差矩阵中NaN的存在使得计算梯度变得非常方便，因为我们不必跟踪哪个运动员在赛季中存在或不存在。例如，先前导出的基线梯度是该季节中所有现有事件-滑冰运动员对{e，s}的残差之和:

因此，基线梯度可以简单地实现为:

baseline_gradient = np.nansum(residuals) # -1692

[np.nansum](https://docs.scipy.org/doc/numpy/reference/generated/numpy.nansum.html)函数有效地跳过了不存在的NaN值(通过将它们视为零)，这符合我们的目的，即只对赛季中现有的事件-选手对的残差求和。

为了计算所有潜在事件分数的梯度，np.nansum函数因其axis参数而再次大有帮助:****

event_gradients = np.nansum(residuals * skater_scores, axis=0, keepdims=True)

让我们分解这行代码，看看它是如何工作的:

• residuals * skater_scores:通过使用*操作，我们简单地将residual矩阵与步骤 1 中初始化的skater_scores逐元素相乘。

然而，由于两个数组的大小不同:矩阵的大小为 4×3 (4 个运动员乘 3 个事件)，而列向量的大小为 4×1(4 个运动员的潜在得分)，numpy 的广播开始运行。实际上，numpy 将水平复制skater_scores向量 3 次，使其与residuals矩阵的维数相同(4×3)，如下图所示:

Highlighted in red: relevant values for the Russian event (RU). **Dotted arrows:** numpy broadcasting of latent skater scores

注意，在广播之后，两个矩阵(现在大小相似)之间的逐元素乘法将在与residuals矩阵中的值相同的位置留下NaN值。这允许np.nansum来计算相关的事件梯度，但这次该函数有两个新参数:

• axis=0:这允许在residuals * skater scores矩阵的行(代表运动员)之间计算总和。实际上，这将对参加每个项目的所有选手的剩余分数和潜在分数的乘积求和，从而得出该项目潜在分数的梯度(见上图)

Highlighted in red: relevant values for the Russian event (RU)

• keepdims=True:对 4×3 residuals * skater scores矩阵的行求和将产生一个 1×3 维的行向量，它正好包含 3 个事件的梯度。从技术上讲，这是 numpy 中 shape (1，3)的一个二维数组，然而默认情况下，np.nansum将结果折叠成 shape (3，)的一个一维数组。因此，设置keepdims=True将使事件梯度保持在 2-D 中，因为稍后将从相同形状(1，3)的event_score数组中减去该梯度。

计算所有潜在溜冰者分数的梯度与计算事件梯度相同，除了我们将残差(residuals大小为 4×3)乘以事件分数(event_scores大小为 1×3，将垂直广播 4 次)。****

skater_gradients = np.nansum(residuals * event_scores, axis=1, keepdims=True)

Highlighted in red: relevant values for the skater Takahito, MURA. **Dotted arrows:** numpy broadcasting of latent skater scores

乘法的结果将通过使用带有参数axis=1的np.nansum跨列(每个选手参加的项目)求和，以获得选手梯度:

skater_gradients = np.nansum(residuals * event_scores, axis=1, keepdims=True)

Highlighted in red: relevant values for the skater Takahito, MURA

步骤 2b:使用梯度更新基线、项目和运动员分数

一旦计算出梯度，使用它们来更新基线、项目和选手的潜在得分只涉及简单的减法(使用学习率alpha=0.0005):

alpha = 0.0005
baseline = baseline - alpha * baseline_gradient
event_scores = event_scores - alpha * event_gradients
skater_scores = skater_scores - alpha * skater_gradients

使用一个简单的for循环，我们可以多次重复步骤 2，在计算梯度和使用它们更新潜在得分之间交替进行:

alpha = 0.0005for i in range(1000):
    # 2a. Calculate gradients
    predicted_scores = skater_scores @ event_scores + baseline
    residuals = predicted_scores - true_scores

    baseline_gradient = np.nansum(residuals)
    event_gradients = np.nansum(residuals * skater_scores, axis=0, keepdims=True)
    skater_gradients = np.nansum(residuals * event_scores, axis=1, keepdims=True)

    # 2b. Update latent scores using gradients
    baseline = baseline - alpha * baseline_gradient
    event_scores = event_scores - alpha * event_gradients
    skater_scores = skater_scores - alpha * skater_gradients

梯度下降算法的收敛可以通过计算每个循环后的 RMSEnp.sqrt(np.nanmean(residuals**2))来监控，并检查它与前一个循环的不同程度。经过 1000 次迭代后，模型的 RMSE 为 4.30，迭代与迭代之间的差异为 10e-9。这表明该算法已经很好地真正收敛。

最后，梯度下降收敛后的潜在选手分数被转换成熊猫系列，这允许早先存储的选手名字被添加回来。请注意，在将二维skater_scores数组(4，1)转换成序列之前，需要用[ravel](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ravel.html)将它折叠成一维数组(4，)。此外，调用sort_values(ascending=False)从最高到最低对运动员的潜在得分进行排序，因为最初的目标是通过这些潜在得分对运动员进行排名。

skater_scores = pd.Series(skater_scores.ravel(), index=skater_names).sort_values(ascending=False)
# Javier, FERNANDEZ     16.701678
# Takahito, MURA        14.101404
# Misha, GE             13.587320
# Alexander, MAJOROV    10.678634

结果

当在 2017 赛季期间将梯度下降算法应用于男性滑冰运动员时，我们可以通过每次迭代的残差和 RMSE 来监控算法的进展。此外，每次迭代中的潜在选手分数可用于建立临时排名，其与世界锦标赛排名的相关性由肯德尔的 tau 测量。在算法期间监控的所有指标(学习率α = 0.0005)显示在下面的动画仪表板中:

Left: heat map of residuals, along with baseline, event, and skater scores at each iteration. Skaters are ordered based on their predicted ranks after 150 iterations. Right (top to bottom): RMSE, Kendall’s tau, and comparison between predicted ranking and world ranking at each iteration (colors are based on world ranking).

从上面的仪表板中，我们可以看到:

• 随着梯度下降算法的运行，残差(预测得分-真实得分)一般会减少，从它们的热图越来越亮就可以看出。我第一次看的时候，觉得很酷！
• 因此，RMSE 也会随着迭代次数的增加而降低(见右上图)。事实上，RMSE 甚至在 20 次迭代后就显著降低，之后开始稳定。
• 在预测排名方面，我们还注意到，随着梯度下降的进行，它开始越来越好地与世界冠军排名相关联(见右下图):连接两个排名的线开始变得越来越不纠结。
• 然而，可以看出，在算法开始附近有一个时刻(大约 25 次迭代左右),两个等级之间的相关性“看起来”最好，这也对应于在同一次迭代中肯德尔τ的峰(中右图)。这意味着，混合模型不需要尽可能好地逼近赛季得分，以便获得运动员的良好预测排名；否则，它可能会过度拟合赛季得分，从而损害预测排名的准确性。
• 这种过度拟合的潜在证据意味着我们应该惩罚混合模型(类似于我们在前一部分对线性模型所做的)。避免过度拟合的另一个策略是当预测等级的肯德尔τ是最高的，也就是驼峰所在的位置时，过早地停止梯度下降。在机器学习的说法中，这被称为提前停止，这两种减少模型过度拟合的策略将在稍后探讨。

以下是梯度下降收敛后的潜在滑手分数；只有获得 2017 年世界锦标赛资格的 24 名选手将被显示，从最高到最低的选手分数排列如下:

Latent skater scores from hybrid model for male skaters in 2017 season (with associated rank)

模型评估

在以学习速率alpha=0.0005和 10，000 次迭代对混合模型运行梯度下降后，模型的最终 RMSE 为 8.86，迭代与迭代之间的差异为 1.7e-8，对于收敛来说足够小。

将该 RMSE 与项目第一部分中的线性模型进行比较时——加法模型为 8.84，乘法模型为 8.88——我们可以看到，混合模型在逼近赛季得分的程度上正好位于之前的模型之间。这也优于季节平均值的基线模型，其 RMSE 为 10.3。

所有型号的每对运动员与真实分数的偏差如下所示。该图中所有 3 个模型之间预测得分的接近程度证实了它们之间相似的 RMSEs。

Deviation from true score for each event-skater pair of the season for all models

就预测排名而言，混合排名再次与早期排名几乎相同(见下图)。事实上，相比于乘法排名，混合排名只是颠倒了只有一对滑手的顺序:Kevin, REYNOLDS和Maxim, KOVTUN，用红色突出显示。

不幸的是，这种逆转与最终的世界冠军结果不符:Kevin, REYNOLDS最终排名高于Maxim, KOVTUN。结果，混合模型在所有 276 个有序对中得到 238 对，比乘法模型少一对。这意味着混合模型的肯德尔τ值略低，为 0.724，而乘法模型和加法模型的肯德尔τ值均为 0.732。然而，这仍然略高于 2017 赛季男子滑冰运动员赛季平均基线模型的肯德尔τ值 0.695(或 276 对中的 234 对)。

在评估训练集中所有 10 年的混合模型之前，让我们首先看看我们如何能够像前面承诺的那样使它对赛季分数的过度拟合更少。

典型惩罚

理论

回想一下梯度下降试图最小化的混合模型的目标函数:

为了惩罚混合模型，我们需要找到一种方法来人为地使潜在的事件和运动员分数变小(同时保持基线分数不变，以便赛季分数仍然是相当接近的)。这样做的一种方法是这样修改目标函数:

The 1/2 factor for the squared terms will make later derivations easier

上述目标函数末端的附加平方项表明:除了最小化季节分数的残差平方，我们还希望最小化潜在事件和潜在选手分数的平方(分别为θ_e和θ_s)，以使这些潜在分数保持较小。潜在得分被惩罚的程度由惩罚参数λ控制:λ越高，潜在得分的平方越小，这些潜在得分就越小。

我们可以相对于这些潜在得分对修改的目标函数进行微分，以在梯度下降期间获得梯度公式(步骤 2a):

Notice that the gradient formula for baseline score stays the same

从上面的梯度公式可以看出，对原梯度下降公式唯一的修改是分别在项目和滑手梯度公式的末尾带有λ的项。

编码

由于这些简单的修改，将这些术语合并到梯度下降的 Python 实现中变得非常容易。下面是惩罚模型的梯度下降步骤 2 的代码(使用惩罚参数lambda_param=10)，与原始模型的变化以粗体突出显示。请注意，对于我们的玩具示例，我们在收敛后获得的惩罚模型的潜在溜冰者分数小于原始模型的分数，这意味着惩罚确实有效。

alpha = 0.0005
**lambda_param = 10**for i in range(1000):
    # 2a. Calculate gradients
    predicted_scores = skater_scores @ event_scores + baseline
    residuals = predicted_scores - true_scores

    baseline_gradient = np.nansum(residuals)
    event_gradients = np.nansum(residuals * skater_scores, axis=0, keepdims=True) **+ lambda_param * event_scores**
    skater_gradients = np.nansum(residuals * event_scores, axis=1, keepdims=True) **+ lambda_param * skater_scores**

    ### 2b. Update latent scores using gradients
    baseline = baseline - alpha * baseline_gradient
    event_scores = event_scores - alpha * event_gradients
    skater_scores = skater_scores - alpha * skater_gradientsskater_scores = pd.Series(skater_scores.ravel(), index=skater_names).sort_values(ascending=False)
# Javier, FERNANDEZ     16.108775
# Takahito, MURA        13.429787
# Misha, GE             12.893933
# Alexander, MAJOROV     9.784296

估价

下面是当惩罚参数λ增加(增加 10 倍)时，混合模型的 RMSE、肯德尔τ和预测排名如何变化的图表:

Learning rate α = 0.0005, Iterations = 10,000

从上图可以看出:

• 随着惩罚参数λ的增加，模型的 RMSE 如预期一样增加，特别是从λ=10 开始，RMSE 几乎是基线模型的两倍(RMSE=10.3)。
• 然而，当我们试图惩罚模型时，预期的肯德尔τ值的增加并没有实现。λ越高，模型的预测排名就变得越差，这可以从肯德尔的τ的减少中得到证明，并且随着λ的增加，预测排名与世界冠军排名相比变得越来越“凌乱”。

当通过绘制不同λ值的基线模型(季节平均值)上的 RMSE 和肯德尔τ的平均差异来检查这种惩罚对训练集中所有 10 年的影响时，出现了相同的画面:

• 如左图所示，当λ增加时，混合模型相对于基准模型的 RMSE 平均差异(10 年以上)通常会增加，尤其是在λ=10 之后。这意味着当惩罚越多，该模型越接近赛季得分。然而，这种行为是意料之中的，因为毕竟惩罚的主要目的是减少对赛季分数的过度拟合。
• 然而，一般来说，混合模型相对于基线模型的肯德尔τ的平均差异随着λ的增加而下降。换句话说，模型给出的预测排名越差，惩罚就越多。这令人失望，但并不令人惊讶，因为惩罚(通过最小化潜在得分的平方)在项目的第一部分对线性模型有同样的影响。
• 此外，与之前的模型类似，与基线模型相比，无惩罚混合模型(λ=0)给出了更差的平均排名预测。然而，这种差异微乎其微:平均而言，Kendall tau 仅下降 0.001，并且不具有统计学意义:当λ=0 时，与基线模型相比，Kendall tau 平均差异的 95%置信区间很容易包含零水平线。

有了这些黯淡的模型惩罚结果，让我们看看是否可以用另一种方式减少模型过度拟合，这有望为混合模型产生更好的预测排名。

提前停止

理论

如前所述，在这种情况下提前停止的想法非常简单:当预测排名的肯德尔τ值最高时，我们只是提前停止梯度下降。实际上，这将防止潜在分数变得太大，就好像我们在原始混合模型中让梯度下降一直收敛一样。

因此，早期停止将使模型在预测赛季得分(更高的 RMSE)方面不如收敛混合模型准确。然而，类似于模型惩罚，这可能会使模型不太适应训练数据，并更准确地对溜冰者进行排名(更高的肯德尔τ)。这在前面的动画仪表盘中可以看到，它监控梯度下降迭代中的这些指标，我在下面的静态版本中复制了这些指标。

RMSE and Kendall’s tau across iterations of gradient descent (for male skaters in 2017 season)

然而，与模型惩罚相比，提前停止提供了一些额外的优势:

• 对于提前停止，模型仅运行—并且提前停止—一次，相比之下，对于模型惩罚，模型运行多次，每次使用不同的惩罚参数λ。
• 对于模型惩罚，我们首先需要猜测惩罚参数λ的一些值来运行模型(例如当λ增加 10 倍时，如前面所示)。如果某个值λ提高了模型的肯德尔τ，我们需要测试这个λ的邻近值，看看是否有任何额外的提高。相比之下，对于早期停止，我们总是可以找到 Kendall 最高的单个迭代，并让该迭代成为停止点。换句话说，在早期停止中，没有“四处寻找”最佳模型超参数的忙乱。

编码

在我们的梯度下降的 Python 实现中，早期停止仅仅意味着将步骤 2 的循环次数限制到一个小的数目——例如 for i in range(10)。之后，可以计算 10 次迭代后的预测排名的肯德尔τ，我们继续尝试其他少量的迭代，并选择具有最高肯德尔τ的一个。

然而，执行早期停止的更实际的方法是运行大量迭代(比如 1000 次)的梯度下降，监视每次迭代的 Kendallτ(类似于我们之前对仪表板动画所做的)，并追溯性地查明 Kendallτ最高的迭代。如果在 Kendall 的 tau 中有几个最大值，我们选择最早的迭代，因为我们希望在较少的模型过度拟合方面出错。

估价

Each dot represents the stopping point (highest Kendall’s tau) for each year. For reference, curves for 2017 are highlighted in red.

当监测每年梯度下降迭代中的 RMSE 和肯德尔τ值时(见附图)，我们可以看到:

• 在几乎所有的 10 年中，如果梯度下降很早就停止，甚至在第 40 次迭代之前，梯度下降就达到最大肯德尔τ(见下图)。相比之下，RMSE 通常在这些点之后继续下降(见上图)，这是有道理的，因为梯度下降的收敛过早地停止了。
• 10 年中只有一年在 40 次迭代后出现最大肯德尔τ，即 2014 年的第 60 次迭代。然而，这可能是侥幸，因为上图显示，在第 40 次迭代之前的几次迭代中，今年的肯德尔τ已经非常接近最大值。
• 然而，尽管大多数年份的理想停止点低于 40 次迭代，但我们可以看到，在一些情况下，一年的理想停止点可能是另一年的灾难(尝试在每个购物点画一条垂直线，并在该迭代中查看其他年份的肯德尔τ)。鉴于肯德尔τ在梯度下降的早期迭代中高度不稳定，这一点尤其正确。

因此，我们必须通过在每次迭代中取平均值来平滑所有 10 年中的肯德尔τ，并找出哪个迭代具有最高的跨年度平均值。更准确地说，对于每次迭代，我们对混合模型的肯德尔τ与每年的基线模型的差异进行平均(超过 10 年)。这使我们能够直接比较早期停止对基线模型的影响，下图(棕色部分)突出显示了这一点:

Difference in Kendall’s tau of hybrid model over baseline model (season averages) at each iteration of gradient descent

从上图中，我们可以看出:

• 在梯度下降的第 47 次迭代时，混合模型相对于基线模型的平均肯德尔τ差异(超过 10 年)达到最大值。
• 然而，快速检查表明，对于大多数年来说，这个迭代已经过了他们的停止点，这意味着他们的肯德尔的 tau 已经走下坡路相当一段时间了。此外，虽然有些年份的肯德尔τ在第 47 次迭代时显示出与基线相比的显著提高(有些年份的肯德尔τ提高了近 0.1)，但仍有一些年份的肯德尔τ无论是否提前停止，与基线模型相比始终表现不佳。
• 因此，毫不奇怪，在 10 年的训练集中，即使是最好的早期停止点，肯德尔的 tau 也只比赛季平均值的基线模型提高了 0.005。更糟糕的是，这种改善在统计学上毫无意义，多年来高度可变的 Kendall tau 改善超过基线就是证明。这相当于一个非常宽的 95%置信区间，很容易包含零水平线(见上图)。

简而言之，即使提前停止似乎也无助于提高混动车型的预测排名。然而，在我们完全放弃混合模型之前，让我们看看我们是否可以以某种方式结合两种策略来减少模型过度拟合——模型惩罚和提前停止——作为改进该模型的最后努力。

结合模型惩罚和提前停止

将模型惩罚和早期停止相结合无非是在计算梯度(模型惩罚)的同时包括惩罚参数λ，同时将梯度下降的迭代次数限制为较小的数目(早期停止)。因此，除了限制梯度下降算法步骤 2 的迭代次数外，Python 代码对模型惩罚没有任何改变——例如，只有for i in range(10)。

对于模型参数的每个组合(梯度下降的迭代次数以及惩罚参数λ),我们可以测量混合模型相对于基线模型的肯德尔τ的平均差异(超过 10 年),并将其绘制在下面的热图中:

Color intensity signifies how much higher on average the Kendall’s tau of hybrid model over the baseline model is

从上面的热图可以看出:

• 即使模型被罚(λ>0)，梯度下降的理想停止点仍在 40–50 次迭代附近徘徊，因为与基线相比，这是肯德尔τ平均最高的区域。
• 然而，随着模型越来越多地被惩罚(增加λ)，当与模型惩罚相结合时，早期停止似乎越来越无效。这在λ=10 及以上的贫瘠荒地上看得最清楚，除了一个小岛在λ=100 时迭代次数接近 80 次。
• 早期停止迭代和λ之间的关系可以以不同的方式可视化，首先通过挑选每个λ的 Kendall(在基线模型上)具有最高平均差的迭代(该λ的理想停止点),并将其性能与如果没有执行早期停止(即如果允许模型完全收敛)的性能进行比较。

The small numbers on the graph show the ideal iteration for early stopping at each λ

以上是男性和女性溜冰者的对比图，显示了非常有趣的模式:

• 对于女性滑冰运动员来说，与之前讨论的男性运动员类似，单独的惩罚往往会降低她们预测排名的准确性(右图中的蓝线):随着λ的增加，肯德尔τ与基线模型的平均差异越小。
• 对于男性和女性来说，提前停止可以使肯德尔的 tau 值比基线模型有最大的提升:男性平均提升 0.005，而女性则提升了 0.009。诚然，这些改善在统计上并不显著，这可以从它们包含水平线为零的宽置信区间中看出。
• 此外，这种通过提前停止的改进适用于每一级惩罚，尤其是当模型被高度惩罚(高λ)时。相比之下，即使早期停止也仍然无法挽救一个高度惩罚的模型，正如当λ增加时早期停止模型的 Kendall tau 改进普遍降低所证明的那样(两个图中的紫色线)。

结论

从上面的分析可以清楚地看出，对于这个问题，提前停止不仅比模型惩罚(如前所述)更容易实现，而且与后者相比，它还具有明显的性能优势。因此，对于两种性别，只选择早期停车来改进混合模式。以下是将其应用于 10 年训练集时的最终结果:

• 请注意，即使早期停止模型与基线模型相比在肯德尔τ方面提供了边际改善，这些改善在很大程度上是不显著的。这可以从肯德尔τ与基线模型的平均差异中看出:这种差异的 95%置信区间对于两种性别来说仍然低于零。
• 此外，由于世界锦标赛排名本身被用于选择最佳模型参数，即梯度下降应该停止时的迭代，模型本身很可能过度拟合 10 年训练集中的数据。因此，在项目的最后部分，我将评估这些提前停止的模型(以及其他模型)与测试集中剩余的 4 年作为最终基准。

资源

在项目的这一部分开发的混合模型最初是受推荐系统中使用的各种矩阵分解技术的启发。一般来说，这些技术会将许多用户对许多商品(例如，在亚马逊上)留下的评级矩阵“分解”成特定于用户的矩阵和特定于商品的矩阵。当这两个矩阵相乘时，它们可以近似原始评级矩阵。

同样的原理在我们的混合模型中成立:评级矩阵只不过是所有事件-溜冰者对的季节分数的数据透视表，而用户特定和项目特定的矩阵只不过是溜冰者和事件潜在分数向量。因此，当乘以skater_scores @ season_scores(加上顶部的基线分数)时，我们可以近似得出赛季的true_scores。您可以查看这篇非常可读的关于推荐系统矩阵分解的文章，其中也包括梯度下降算法的推导，该算法与本项目中使用的算法非常相似，通常被称为 FunkSVD 。

也就是说，我们目前的混合模型只涉及一个单一因素:每个运动员和每个项目只由一个潜在的分数来表示。相比之下，推荐系统中的矩阵分解往往涉及多个因素。在项目的下一部分中，我将描述多因子矩阵分解如何应用于我们的排序问题。

用前馈神经网络预测食物份量

原文：https://towardsdatascience.com/predicting-food-serving-sizes-with-a-feed-forward-neural-network-2d8d40d82f72?source=collection_archive---------30-----------------------

在我在大会的沉浸式项目中，我必须完成一个顶点项目，一个完全由我自己设计的项目，除了满足几个关键标准之外，没有任何关于做什么或如何去做的指示。所以我决定专注于我非常了解并且热爱的东西——食物！我找到了一个由美国农业部维护的数据集，其中包含了该国超过 260，000 种包装和品牌食品，被恰当地命名为美国农业部品牌食品数据库。他们对数据库有相当开放的 API 访问，但幸运的是，完整的 API 也可以作为 zip 文件下载，节省了我 30 多个小时来 ping 他们的服务器。它保存了大量信息，从制造商和成分的基本信息到更具技术性的信息，如用于确定每种营养价值的特定实验室方法。引起我注意的一件事是这些产品的食用量差异很大。

Spam! Credit: Lhe3460, CC BY-SA 4.0

问题

现在，对于那些不在食品科学游戏中的人来说，每份食物的份量实际上有食品和药物管理局的指导方针，食品公司必须遵守，称为习惯食用的参考量。这是一个灰色地带，因为它们不仅仅是一刀切的(老实说，我不能说这是一个双关语)，但它们必须被用作公司的出发点。因此，开发团队仍有大量工作要做，以确定食品的最终食用量。这成了我项目的中心目标:建立一个模型来预测一种食品的份量。

清洁

第一步是合并所有的数据。因为这是一个关系数据库，所以我要找的信息被分散在多个表中。我将食品名称和它们的详细信息结合起来，比如品牌所有者、市场类别和成分。添加营养成分稍微有点棘手，因为我必须创建一个数据透视表，将它们与产品放在同一格式中。每种食物的能量含量有两个不同的单位:千卡和千焦，这有点小问题。千焦在欧洲更常见，而美国使用千卡，这是一个大写的 C 卡路里(是的，卡路里和卡路里之间有区别，不，我也不喜欢它)。我编码了一个快速转换，然后我就回到了正轨。最终的组合表如下所示:

然而，完整的营养表并没有被剔除。它必须被缩小到前 10 个最常见的。这主要是由于营养标签的一个重要事实:它们并不总是列出相同的东西。以外的常见营养素(脂肪、蛋白质、胆固醇等。)，必输项很少，大部分维生素和矿物质在大多数情况下是可选的。所以，在 97 种不同的营养成分中，绝大多数超过 90%是零的情况下，我不得不坚持强制输入。

从那以后，我还得清理许多遗留的东西。一个例子是许多卡路里计数是空的。起初，这似乎是一件大事，因为估算不是一件轻而易举的事情，但如果你对营养学很感兴趣，你应该知道实际上有一种非常简单的方法来估算卡路里，通过使用食物的三种主要营养素——脂肪、蛋白质和碳水化合物。总的来说，这些是我们每天摄入的卡路里的提供者，每种类型都有特定的量。蛋白质和碳水化合物每克都给我们的身体提供 4 卡路里，而脂肪每克提供 9 卡路里。再一次，一些快速的数学和代码是有序的，我们丢失的卡路里已经被考虑。

当我开始深入研究这些数据时，我注意到一些有趣的事实。食用量的分布有一个长的右尾，大部分食物在 0-50 克之间。由于形状的原因，我决定最好将数据转换成对数形式，以便进行建模。查看市场类别的影响也揭示了最大份量的食物类型，液体、汤类和冷冻食品构成了前 10 项。这很有意义，因为与许多其他食物相比，液体相当稠密，像冷冻主菜这样的东西可能会取代其他几种食物，本质上是多份食物(如意大利面、蔬菜和肉类)。

建模

我最终使用的功能是全营养事实，和市场类别(这是 dummied)。我曾试图对配料表执行 NLP 以进一步扩展功能，但由于课程的时间限制，我无法有效地优化它们。现在我有了更多的时间，我想回到这个话题上来。

Basic Neural Network. Credit: Glosser.ca, CC BY-SA 3.0

使用并评估了几种模型，但最终我决定使用前馈神经网络。如果你需要复习一下神经网络，可以查看一下简单设置的图表。它们本质上是一系列复杂的逻辑回归(线条)，由您决定有多少隐藏层，以及每层中有多少节点(圆圈)。有一些关于建立神经网络的指导方针，但是可能性确实是无穷无尽的。

带着为我的神经网络寻找最佳参数的艰巨任务，我决定创建一个工具来帮助我加快这个过程。在 SciKit-Learn 中，有一个非常有用的工具，GridSearchCV可以为各种各样的模型找到最佳参数。您设置您正在构建的模型，然后可以输入广泛的超参数，GridSearchCV将构建这些模型，并告诉您哪个组合是最佳的。对于 Keras 中的神经网络，没有与此功能等效的功能，所以我构建了自己的功能。它是三个不同功能的组合:一个用于接受参数字典并迭代出每个排列，一个用于基于这些参数排列构建每个不同的神经网络，另一个用于组合这些功能，执行评估并跟踪最佳模型。这是很多代码，所以我只在这里包括了最后一块，但如果你有兴趣看到完整的细节，请查看回购。

**def** nn_grid_search(
    X,
    y,
    grid_params,
    random_state=42
):
    *### this will make a series of FFNN models* 
    *### and return the one with the best score as set below*
    *### currently set to test r2 score*

    *# list of all parameter combinations*
    all_params = permutate_params(grid_params)

    *# creating vars with to update each iter*
    best_model = **None**
    best_score = 0.0 
    best_params = **None**
    best_history = **None**

    *# train/test split the data*
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=random_state)

    *# scaling data*
    ss = StandardScaler()
    X_train_sc = ss.fit_transform(X_train)
    X_test_sc = ss.transform(X_test)

    *# looping through the unpacked parameter list*
    **for** i, params **in** enumerate(all_params):

        *# keeping track of which model we're running*
        print(f"Building model {i + 1} of {len(all_params)}")

        *# bulding the model*
        model, history = build_model(
            params_dict = params,
            X_train = X_train_sc, 
            X_test = X_test_sc, 
            y_train = y_train, 
            y_test = y_test
        )

        *# making preds and scoring*
        test_preds = model.predict(X_test_sc)
        score = metrics.r2_score(y_test, test_preds)

        *# checking if the score beats the current best*
        *# updates vars if true*
        **if** score > best_score:
            print("***Good R2 found: **{:.2%}*****".format(score))
            best_score = score
            best_model = model
            best_params = params
            best_history = history

    *# loop is done, return the best model*
    **return** {
        "best_model"   : best_model,
        "best_score"   : best_score,
        "best_params"  : best_params,
        "best_history" : best_history,
        "test_preds"   : test_preds
    }

我决定给我的模型设置 3 个隐藏层，试图在避免过度拟合的同时平衡深度学习，并着手确定我的其他参数。我给了每一层一个 ReLU 激活，这通常被认为是这些网络最好的激活功能。接下来，我需要在每个隐藏层中选择多少个节点。这通常是一个基数为 2 的数字，因此这已经有助于缩小选择范围。此外，让您的第一个图层具有与输入(即您的要素)数量相似的结点是一个好主意。总共不到 200 个特征，我选择了一系列合适的数字。由于拥有一个金字塔形状的网络(每个后续层的节点更少)也是有帮助的，所以我有一个清晰的方向来建立我的完整参数字典。

在运行搜索并获得最佳模型的最终结果后，我采用这些参数，并通过增加运行的时期和应用各种正则化技术(L2 惩罚、退出和提前停止)来继续微调模型。我的模型的最终参数集如下:

'best_params': {
 'first_layer_nodes': 256,
  'first_dropout_rate': 0,
  'second_layer_nodes': 128,
  'second_dropout_rate': 0,
  'third_layer_nodes': 64,
  'third_dropout_rate': 0,
  'reg': 0.001,
  'epochs': 75,
}

评价

我评估神经网络的第一步是可视化它运行的全部 75 个时期的损失曲线，用 Adam 优化器监控均方误差作为损失。

我发现曲线非常平滑，训练集和测试集的损失值非常一致，没有任何明显的差异。

然后，我使用两个关键指标来分析模型的表现:调整后的 R 和均方根误差(RMSE)。使用调整后的 R 非常有用，因为它不仅能告诉您模型覆盖了多少数据方差，还能说明您拥有多少要素，从而确定它们是否真的有助于您的模型。另一方面，RMSE 为您提供了一个可直接解释的值，用于查看模型的预测值通常与真实值之间的误差。

这些分数彼此非常接近，也比我构建的其他模型的分数好得多，我觉得我已经在这段时间内尽可能地充分优化了我的模型。不幸的是，当我们现实地看待这些结果时，预测的食用量通常与实际食用量相差 24 克以上可能会有些问题，因为平均食用量约为 64 克(相差 38%)。

包扎

最后，我不相信这是一个完全的成功，因为我的模型有很多错误，它可能不会比仅仅遵循 FDA 的指导方针更有帮助。由于我无法将这些指导方针整合到数据中，我无法确定我的模型如何直接与它们相抗衡，所以一切都应该持保留态度(这绝对是一语双关)。

然而，它确实显示出了希望，我认为我可以做很多事情来改进这个项目。我要做的第一件事就是花些时间完善 NLP 步骤，从而整合这些成分。此外，如果这是一个生产设置，我将实现 API 访问，以跟上对数据库或目标相关食品的特定组所做的任何更改。

这绝对是一个有趣的项目，我学到了很多关于完整的从前到后的数据科学工作流程，清理数据，尤其是构建神经网络。

如果你碰巧还在读这篇文章，那么我为你坚持和我一起探索食物世界而鼓掌！请回来查看更多帖子！

你可以在这里找到全项目回购。

你也可以在 LinkedIn 上找到我。

用 XGBoost 预测足球运动员的商业价值

原文：https://towardsdatascience.com/predicting-football-players-commercial-value-with-xgboost-d0670d9e9d2e?source=collection_archive---------13-----------------------

Pic by Jannes Glas - Unsplash

商业价值真的只依赖于纯竞技场属性吗？还是另有隐情？

如果你和我一样，你会觉得有必要了解所有事物是如何工作的，如果你对数据科学感兴趣，你会觉得有冲动去预测所有可以预测的事情。在这种情况下，我会告诉你如何预测足球运动员的商业价值仅仅依靠他们的足球技术。我们将使用 XGBoost 来这样做，并在这样做的同时了解更多的库。

我使用 Jupyter 笔记本在 Python 中进行了分析。笔记本可以在我的 GitHub 页面上找到，我会不时地加入一些代码片段，让这个故事更有教育意义。我鼓励读者拿起笔记本，边阅读文章边阅读。

数据

该数据集包括超过 18000 个足球运动员的条目，按照价值排序，从最有价值到不太有价值。虽然数据集与 FIFA '19 视频游戏相关，但其玩家商业估值和玩家游戏技能评级非常准确，因此我们可以假设我们正在处理现实生活中的玩家数据。

导入

首先，让我们获取数据集。我最初是从 Kaggle 获得的，但是你也可以从我的 GitHub 页面获得。

安装 pandas 以便能够操作数据，然后直接从 URL 导入数据集:

import pandas as pdurl = '[https://raw.githubusercontent.com/estebanvillaturek/fifa19/master/fifa19_data.csv'](https://raw.githubusercontent.com/estebanvillaturek/fifa19/master/fifa19_data.csv')data = pd.read_csv(url)data.head()data.shape

正如我们所看到的，我们得到了一个包含超过 18k 行和 89 列的数据集，包括所有类型的数据格式。

探索

我们主要对“价值”一栏感兴趣，它包含了每个球员的商业价值。我们注意到一些事情:它的条目被格式化为字符串，前面有一个€符号，百万和千分别用大写的 M 和 K 表示。这意味着，我们将有一些争论和重新格式化要做，因为我们想得到所有的数值作为数值，以便建立和运行预测模型。

重新格式化“值”

我们可以使用一个简单的 lambda 函数来移除€符号。在这种情况下，我们求助于 lambda 函数(也称为匿名函数)，因为我们只需要对整个列应用一次性函数，并且不用任何东西替换€符号。

# Remove €data.Value = data.Value.apply(lambda x: x.replace('€', ''))

接下来，由于值的范围从数百万到数千，我们需要找到一种方法，首先找到以“K”结尾的条目，然后删除它，将字符串转换为浮点数(数字数据类型)，然后乘以 1000。然后，我们将对百万欧元的条目重复相同的过程。

最好的方法是使用 for 循环迭代该列。创建迭代器之前的最佳实践是使用 next()函数模拟迭代器，该函数返回第一次迭代，这样我们就可以使用结果来创建最终迭代器。注意，为了迭代 pandas 数据框列，即一个序列，我们调用。items()方法。

# preview iteration next(data.Value.items())

第一次迭代的结果是一个元组，由第一个元素索引和第二个元素值本身组成。有了这些信息，我们就可以构建第一个迭代器来重新格式化千欧元单元格。

# reformat thousands cellsfor i,j in data.Value.items():
    if 'K' in j:
        data.Value.iloc[i] = float(j.replace('K', '')) * 1000

现在，让我们对剩下的百万欧元细胞做同样的事情。

# Reformat millions cellsfor i,j in data.Value.items():
    if 'M' in j:
        data.Value.iloc[i] = float(j.replace('M', '')) * 1000000

将“值”转换为数值

既然“Value”列中的所有条目都可以表示为数字数据类型，我们就可以相应地转换它们，这样我们就可以获得回归模型的目标(没错，我们将预测一个连续的因变量，因此我们将编译一个回归模型，但稍后会详细介绍)。

# to_numericdata.Value = pd.to_numeric(data.Value, errors = 'coerce')

创建新的全数字数据框

现在我们知道我们将在练习中使用回归变量，我们将创建一个新的、更具体的数据框，该数据框将仅包含数字变量，并将由目标变量“值”和所有预测变量、每个玩家的游戏技能的排名分数(每个玩家总共 34 个排名能力)组成。

# first, the target variablescores = pd.DataFrame(data.Value) # then, slice the ranking scores for players' abilitiesdata1 = data.iloc[:, 54:88] #finally, join the the two subsetsscores = pd.concat([scores, data1], axis = 1)

描述分数&留意偏斜度

由于我们处理的是数值，我们必须检查数据是如何分布的。这是理解数据行为的基础，也是理解我们是否必须对数据进行规范化或标准化的基础。请记住，归一化将我们的数据转换为从 0 到 1 的范围，正则化将我们的数据转换为正态分布，平均值为 0，标准差为 1。

# The describe method returns an overview of the descriptive stats of all our columns.scores.describe()

我们可以从目测数据中看出，有轻微的偏态倾向，因为中位数和平均数总是分开的(在正态分布中，中位数和平均数是相同的)。在这种情况下，偏斜的一个特殊情况是“值”变量。它有相当多的离群值(梅西、罗纳尔多等)，将平均值拉向右边。我们可以使用 Seaborn 直观地看到分布，但首先我们需要检查数据中是否有缺失值。

检查 NaNs

检查缺失值、nan 等的简单方法。，就是使用。isnull()方法，然后按列对所有缺失值求和。像这样，我们可以看到每列有多少个 nan，因此，决定如何处理它们，即是否删除它们，估算它们，等等。

# sum up all non-values scores.isnull().sum()

结果是惊人的:尽管我们的预测变量每列只有 48 个缺失值，但我们的目标变量“值”有将近 6k 个缺失值。这意味着，如果我们决定删除所有的 nan，我们将丢失几乎三分之一的数据集。).那是不行的。

估算中值以避免数据丢失

或者，我们将求助于数据插补，这是一种允许我们用另一个值替换所有缺失值的技术。对于此示例，我们将使用每列的中值来估算每列中的 NaNs。我们选择中位数，因为它根本不会影响每个变量的分布，也许除了“价值”之外。让我们检查一下。

# to not lose around 6k data points, let's explore missing data imputationscores.iloc[:, :].fillna(scores.iloc[:, :].median(), inplace = True)print(scores.isnull().sum())

现在，数据集中没有丢失的值，所以我们可以开始可视化地研究分布。

检查单个分布

这里我们也可以求助于 for 循环，我们可以要求机器绘制每个变量的分布。

# Explore the distribution of each variableprovisional = scoresfor i, j in provisional.items():
    sns.kdeplot(j, shade = True)
    plt.figure()

正如我们所看到的，“价值”仍然是高度倾斜的，我们的大多数变量倾向于正态分布(除了那些相对于守门员能力排名的变量)，其中一些变量倾向于双峰。

正常化(如果需要)

首先，让我们把目标变量和预测变量分开，因为它们需要不同的处理方法。

# Separate examples and targetX = scores.iloc[:, 1:]y = scores.iloc[:, 0]

为了练习的缘故，我们将标准化和规范化(标度)我们的预测因子，以便能够评估治疗的差异。我们将使用 sklearn 库进行预处理任务。

# Now we’ll first normalize the data so that we can have it in a 0 to 1 range and then standardize itfrom sklearn import preprocessing#nrm_X = preprocessing.normalize(X)#std_X1 = preprocessing.scale(X)#std_X = preprocessing.scale(nrm_X)

现在，让我们直观地比较我们的原始“交叉”预测器及其缩放和标准化版本。

# Assess the difference visuallysns.kdeplot(std_X[:, 0])sns.kdeplot(std_X1[:, 0])

Original Crossing variable

The orange distribution is the scaled version w/o normalization. The blue distribution is the same scaled variable, only after being normalized first.

我们可以看到，虽然差别不大，但蓝色分布的方差要小一些。让我们继续这个标准化和缩放的预测数据。

对数转换目标变量

正如我们前面看到的，“价值”是相当扭曲的，因为有天文数字估计的商业价值的突出离群值:梅西，罗纳尔多，小内马尔，等等。这是有问题的，因为绝大多数球员的价值都远低于百万欧元。

要在数据中解决这个问题，我们必须对“值”变量进行对数转换。在处理随时间增长的经济或金融数据(工资、财富、GDP 等)时，对数变换非常有用。)因为它们往往是指数增长，而不是线性增长。但是要小心，不要对零进行对数转换，否则会得到 inf 值，这将使您无法绘制或继续分析您的数据！

# Now we'll have to rescale the target variable - we'll log-transform it - it's crazy skewed!plt.figure(figsize=(15,10))
sns.kdeplot(y, shade = True)

由于数据中有 252 个零，我们将再次估算中值，以便能够对变量进行对数转换。

y[y == 0].count() # we have 252 zeroes in 'Value' - we have to impute the with the mediany[y == 0] = np.median(scores.Value)y[y == 0].count() # ready, now we can log-transform the variable w/o getting any - inf values

现在，让我们对目标变量 y 进行对数变换，并对其值进行指数运算，以练习将它们恢复到原始比例。

# do the log-transformation - looks much better!y_log = np.log1p(y)plt.figure(figsize = (15,10))
sns.kdeplot(y_log, shade = True)

Log-transformed ‘Value’ variable

# now we can just do the exponential of any value of y to get back the original valuesnp.expm1(y_log[0])

构建 XGBoost 回归模型

首先，我们可以直观地评估随着每个游戏技能等级分数的增加，价值是如何增加的。为此，我们可以再次使用 For 循环，这样我们就可以创建 34 个散点图，每个散点图对应玩家的额定游戏能力。这里我们将使用未缩放的数据集，但是为了更好的可视化，我们将对图中的 x 变量进行对数变换。

# Visualize the relationships within the datasetx_col = "Value"
y_columns = scores.loc[:, scores.columns != 'Value']for y_col in y_columns:figure = plt.figure
    plt.figure(figsize = (15,10))
    ax = sns.regplot(x=x_col, y=y_col, data = scores, logx = True)
    ax.set_xlabel(x_col)
    ax.set_ylabel(y_col)
    ax.set_title("{} & {}".format(x_col, y_col))plt.show()

我们可以清楚地看到，随着每项技能的增加，价值也会增加，这是很自然的。我们还可以看到，在一大群更为普通的玩家中，具有更高商业价值的顶级玩家是如何成为显著的离群者的。

现在让我们来看看 XGBoost 模型。正如这里完美解释的，XGBoost 是 2016 年开发的一种基于决策树的卓越算法，它利用梯度下降和不同硬件相关计算技术的力量来提升模型，使其符合结构化集合的数据，这是其他算法无法做到的，至少在今天是如此。它非常适合构建分类或回归模型。

首先，我们将拆分数据，

# We will use a 20% of the data set as the test setfrom sklearn.model_selection import train_test_splitX_train, X_test, y_train, y_test = train_test_split(std_X, y_log ,test_size=0.2, random_state = 42)

然后我们将编译并拟合回归模型。

# Compile and fit the modelimport xgboostmodel = xgboost.XGBRegressor(objective='reg:squarederror')model.fit(X_train, y_train)

我们用我们的测试集来预测。

# Now we predict based on the fitted modelpreds = model.predict(X_test)

最后，我们基于适当的准确性指标来评估我们的模型。在这种情况下，我们将使用平均绝对误差(MAE ),它代表我们的预测和我们的真实测试目标之间的差异的绝对值。

# Accuracyfrom sklearn.metrics import mean_absolute_errormae = mean_absolute_error(y_test, preds)

请记住，由于我们的预测和测试目标是对数转换的，我们需要取 MAE 的指数来获得原始单位的度量。

np.expm1(mae)

这给了我们大约 0.75 €的平均相对误差，相当不错，对吧？

结论

我们刚刚看到了一种方法，只根据球员不同的足球能力评级来预测球员的商业价值，而不考虑其他因素。虽然预测本身相当精确，但我们在模型中的变量解释的方差并不太大(R2 上升到只有大约 35%)，这强调了清楚你想通过模型实现什么的重要性:要么是预测准确性，要么是通过自变量解释因变量中的方差来彻底解释现象，这是机器学习和统计学之间的基本差异之一。

奖励曲目

1. 想想职业足球历史上所有可以避免的转会失误！看看《经济学人》的这篇文章，感受一下成功(和失败)的球员在完成转会时的评估对球队表现的巨大影响。
2. 如果你想知道在预测商业价值时哪些玩家能力是最关键的，你可以探索 XGBoost 提供的特性重要性指标。不过要小心:看看这篇文章，感受一下重要性度量告诉你什么！

让我知道你的想法，别忘了留下评论！

从零开始使用快速 AI API 通过 RNNs 预测未来的医疗诊断

原文：https://towardsdatascience.com/predicting-future-medical-diagnoses-with-rnns-using-fast-ai-api-from-scratch-ecf78aaf56a2?source=collection_archive---------6-----------------------

全面 pytorch 实现医生 AI 论文使用电子健康档案

在本教程的第一篇第一部分中，我们创建了 Edward Choi 等人的医生 AI:通过递归神经网络预测临床事件论文(2016) 的粗略模板。在本教程中，我们使用 Fast.ai 自下而上方法对其进行了进一步处理。该代码功能齐全，有关数据处理的详细信息可在第一部分中获取。

详细代码: Github

加载数据

关于数据集:

本研究将利用 MIMIC III 电子健康记录(EHR)数据集，该数据集包含 38，645 名成人和 7，875 名新生儿的 58，000 多份住院记录。该数据集是从 2001 年 6 月至 2012 年 10 月在贝斯以色列女执事医疗中心的去识别重症监护病房住院病人的集合。在第一部分中可以找到所使用的数据预处理步骤的详细演练。

数据预处理数据集将被加载，并按75%:15%:10%比率分成训练、测试和验证集。

Data loading function

填充序列:处理可变长度序列

使用在第一部分中创建的人工 EHR 数据，我们将序列填充到每个小批量中最长序列的长度。为了帮助更深入地解释这一点，让我们看一看在第一部分中创建的Artificial EHR data。

使用人工生成的 EHR 数据进行详细解释

在这里，您可以看到我们有一个包含两个列表数组，每个列表代表一个独特的患者。现在，在每个列表中有一系列列表，每个列表代表一次独特的访问。最后，编码的数字代表每次就诊时分配的诊断代码。值得注意的是，鉴于每位患者病情的独特性，指定的就诊和诊断代码都有variable length序列。因为 EHR 数据本质上是纵向，我们通常对了解患者的风险或随时间的进展感兴趣。当使用表格数据处理时，这些嵌套的依赖于时间的variable length序列会很快变得复杂。回想第一部分的下图，详细描述了每次就诊日期与就诊期间指定的诊断代码之间的映射关系。

Patient Sequence Encodings

Python Pickled List of List containing patient visits and encoded diagnosis codes

那么我们到底用这个嵌套列表填充什么呢？

让我们分解填充函数:

1. lenghts = np.array([len(seq) for seq in seqs]) - 1这里神秘地从长度中减去 1，在作者的笔记中，他提到visit和label文件必须匹配，因为算法会考虑推理时间的时间延迟。

这是什么意思？考虑到数据的结构，每个患者记录中的最后一次就诊将被删除。如此处所示:

Removing the last visit for inference

旁白:在字符级 RNN 中处理可变长度序列

如果这是一个角色级别的问题，我们就说[ Sparkle、Dorian、Deep、Learning。这些序列首先按长度降序排列，并用零(红色)填充，其中每个字母代表一个令牌。如下所示:

Variable length sequence padding

EHR 数据:

然而，对于这种形式的 EHR 数据，给出了我们当前的问题，而不是每个编码的诊断代码代表一个唯一的令牌。在这种情况下，每次访问代表一个令牌/序列。因此，使用与字符级 RNNs 相同的方法，我们首先按照患者就诊降序排列每个小批量。在这种情况下，患者 1 具有最长的就诊历史，共有两次就诊，而患者 2 的就诊将被填充到最大长度 2，因为它是最长的序列。如下所示:

Padding EHR data

现在，我们已经解决了可变长度问题，我们可以继续对我们的序列进行多一热编码。这将产生所需的 S x B x I 尺寸(序列长度、批量大小、输入尺寸/vocab)。

在这里，我们可以很容易地看到，序列将代表每个小批量中就诊历史最长的患者，而所有其他人将被填充到这个长度(红色)。根据所需的批次大小，批次大小将代表每个时间步输入多少患者序列。最后，内部列表将被编码为词汇表的长度，在本例中是整个数据集中唯一诊断代码的数量。

Multi-one hot encoded sequences

标签

以确保标签被移位一个序列，以便算法可以准确地预测下一个时间步长。作者通过确保训练数据排除每个患者历史中的最后一次就诊来解决这一问题，使用这种逻辑for xvec, subseq in zip(x[:, idx, :], seq[:-1]):，其中我们采用每个患者就诊记录seq[:-1]中除最后一次就诊之外的所有就诊。对于标签，这意味着序列将从患者的第二次就诊开始，或者按照 python 的索引风格，第一个索引for yvec, subseq in zip(y[:, idx, :], label[1:])，其中标签label[1:]移动一位。

Label time step lag

什么是掩蔽，它有什么作用？

屏蔽允许算法知道真正的序列在 one-hot 编码数据中的位置，简单地说就是忽略/过滤掉填充值，在我们的例子中填充值为零。这使我们能够轻松处理 RNNs 中的可变长度序列，它需要固定长度的输入。是怎么做到的？还记得lengths变量吗？该变量以降序存储每个患者序列的有效长度(回忆:在移除每个记录中的最后一个序列以进行推断后，例如，患者 1 有 3 次就诊，但长度将仅反映 2 次)。然后，代码mask[:lengths[idx], idx] = 1.中的逻辑用 1 沿着行填充我们的归零张量，以匹配从最大到最小的每个患者序列的长度。

lenghts_artificial → array([2, 1])

mask_artificial → tensor([[1., 1.], [1., 0.]])

数据加载器和采样器

Dataset类是一个抽象类，表示 x 和 y 对中的数据。

Sampler类随机打乱训练集的顺序(验证集不会被随机化)。此外，它保留创建完整批次所需的准确序列数量。

DataLoader类结合了数据集和数据采样器，后者遍历数据集并抓取批处理。

嵌入层

在将输入数据呈现给 GRU 之前，Custom_Embedding类用于将高维多热点编码向量投影到低维空间。在这一步中，作者使用了两种方法

1. 随机初始化，然后在反向推进期间学习适当的 W(emb)W(emb)权重

2.使用 Skip-gram 算法初始化预训练嵌入，然后在 back-prop 期间优化权重

在本文的实现中，我们使用了第一种方法。因此，创建了Custom Embedding类来在嵌入层上应用 tanh 激活。

Custom embedding Layer

脱落层

在本文中，作者使用了由 Srivastava (2014) 首次引入的辍学的简单应用。虽然这种方法效果很好，但它影响了 RNNs 保持长期相关性的能力，因为我们没有在每个时间步长保持相同的掩码。为什么这很重要？很简单，如果我们在每个时间步随机采样一个新的掩码，它会干扰我们的 RNNs 连接，使网络难以确定哪些信息可能是长期相关的。在这种方法中，我测试了 Gal & Ghahramani (2016)提出并由 Merity (2017) 进一步开发的 LSTMs 技术。在这里，他们提出通过在 LSTMs 中的多个时间步长上使用相同的漏失掩码来克服上述与随机采样相关的问题。在这里，我将应用相同的方法在每层(两层)之间的 GRU 上。

Dropout Layer

艾医生:通过递归神经网络预测临床事件

尽管 LSTMs 很受欢迎，也很受青睐。本文使用了 GRU 架构，因为它简单并且能够获得与 LSTMs 相似的性能。本文中使用的数据集包含263, 706 patients，而我们的数据集(MIMIC III)总共包含7537 patients。然而，作者证明了在一个医院系统缺乏训练像 AI 博士这样的深度学习模型所需的大规模数据集的情况下，迁移学习可能是一个可行的选择。使用以下架构，我的兴趣在于对患者未来诊断代码的预测。然而，人们可以很容易地推断出该算法来预测诊断和就诊间隔时间。

Model Architecture

GRU 层:

这个类使用了EHR_GRU单元格类，并允许在期望的层数上迭代。

损失函数:

用于评估模型性能的损失函数包含交叉熵的组合。每个小批量的预测损失被标准化为序列长度。最后，L2 范数正则化应用于所有的权重矩阵。

模型参数:

这里使用的参数选自 AI 博士论文中使用的参数。这种方法和我在这里介绍的方法之间的主要区别是，我对 RNNs 使用了更新的 drop out 方法。

num class = 4894
input dimsize = 4894
embSize = 200
hiddenDimSize = 200
batch size = 100 num layers = 2

加载数据:

需要注意的是，您希望将序列和标签的同一个文件传递到load_data函数中，因为模型会在内部负责调整预测的时间步长。

训练和验证循环

我的实现与论文算法的比较:

我在论文的算法上运行了相同的序列，这是用 theano 和 python 2.7 编写的，这里可以看到 10 个时期后的最佳交叉熵分数约为 86.79，而 my 为 107。虽然，通过一些超参数调整和优化，我并没有表现得更好，但算法肯定会表现得更好。

Dr. Algorithm results for comparison

观察结果:

正如你所看到的，我们的训练和验证损失几乎是一样的，实际论文中使用的数据是如此之少。如果不过度拟合，可能很难获得更好的性能。然而，本教程的目的是提供一个如何使用 EHR 数据驱动洞察力的详细演练！

完整脚本

后续步骤:

1. 使用 Fast 添加回调。人工智能的回调方法来跟踪训练数据
2. 尝试不同的初始化方法

致谢:

1. Fast.ai(雷切尔·托马斯、杰瑞米·霍华德和令人惊叹的 fast.ai 社区)
2. 多里安·普勒里

使用随机森林模型预测鱼类样本的地理来源

原文：https://towardsdatascience.com/predicting-geographic-origin-for-fish-samples-using-random-forest-models-3a0f791d0f5c?source=collection_archive---------24-----------------------

机器学习概念如何支持渔业管理

问题

我试图展示一种分析的效用，这种分析根据鱼的耳骨形状，对来自特定物种的鱼样本进行分组。基本概念是，特定种类的鱼，比如鳕鱼，有独特的耳骨形状，可以用来识别它们来自哪个地理区域。我想展示我设计的特征能多好地预测每个样本的地理来源；本质上，这项技术有合适的预测能力吗？

数据收集

你可能会问，如何量化耳骨的形状？我用离散小波变换(类似于常用的傅立叶变换)用余弦波来描述每块骨头的轮廓形状(图 1)。这种转换给了我 64 个系数来描述每条鱼耳骨的形态差异，这是我们的特征。

Figure 1: Ear bone and the derived outline using a Wavelet transformation

回应:地理起源

特征:描述骨骼形状的小波系数

数据处理

我在不同的人群中进行了不平等的抽样调查；在不涉及生物学细节的情况下，我需要在没有人口统计学细节(如长度、年龄等)的情况下，分离出地理来源对骨骼形状的影响。)促成关系。我使用重复的 ANCOVAs 来消除人口统计变量随地理区域显著变化的特征，并应用 Bonferonni 调整来最小化重复分析中第一类错误的累积。

# function that will derive the p value for covariates
ancova.pvalues <- function(resp, pred, cov){
  return(unlist(summary(aov(resp~pred*cov)))['Pr(>F)3'])
}# apply to each feature (given predictor of region (pops) and covariate of length)
p.values<-0 # stores values
for (i in 1:length(colnames(Wavecoefs))){
  p.values[i]<-ancova.pvalues(unlist(Wavecoefs[i]), pops, length_cm)
}
which(p.values<0.05) # which features should we omit

有些特征是偏斜的，所以我对那些偏斜度大于 0.75 的特征应用了 box cox 变换。

很少功能包含 NAs 我用每个特性的平均值代替了 NAs。

我使用方便的方法“train_test_split”将数据集分为训练集和测试集，分别使用 80%和 20%的数据。

from sklearn.model_selection import train_test_split
train_X, test_X, train_y, test_y = train_test_split(df, resp, random_state = 0, test_size=.2)

我的响应变量——地理区域——是不均衡抽样的(图 2)。我想展示我设计的功能可以预测来自多个不同地区的鱼的原产地，并且我想最小化样本大小的变化对模型预测的影响。为此，我随机对最常见的类(海湾，蓝色)进行欠采样。

Figure 2: Sample size for each geographic region (1-Gulf of Mexico, 2-West Atlantic, 0-East Atlantic)

包不平衡学习为此提供了一些有用的方法；我使用“RandomUnderSampler”来创建一个更加平衡的训练数据集，以便在其上拟合我的模型(图 3)。

from imblearn.under_sampling import RandomUnderSamplerrus = RandomUnderSampler(return_indices=True)
X_rus, y_rus, id_rus = rus.fit_sample(train_X, train_y)

Figure 3: Sample size for each geographic region after undersampling

建模

鉴于描述骨骼形状的特征，我使用随机森林分类器来预测样本来自的区域。首先，我确定了最佳超参数值:

max_features :每次分割要考虑的最大特征数

max_depth :任意树的最大分割数

min_samples_split :分割一个节点所需的最小样本数

min_samples_leaf :每个叶节点所需的最小样本数

自举:数据集是自举还是整个数据集用于每棵树

标准:用于评估每次分割质量的函数

max_features = [‘auto’, ‘sqrt’, ‘log2’]
max_depth = [int(x) for x in np.linspace(10, 110, num = 11)]
max_depth.append(None)
min_samples_split = [2, 5, 10]
min_samples_leaf = [1, 2, 4]
bootstrap = [True, False]
criterion= ['gini', 'entropy']
grid_param = {'max_features': max_features,
               'max_depth': max_depth,
               'min_samples_split': min_samples_split,
               'min_samples_leaf': min_samples_leaf,
               'bootstrap': bootstrap,
               'criterion':criterion }

sci-kit 学习模块有一个简便方法“GridSearchCV ”,可以通过交叉验证找到最佳的超参数值。我用了 5 折的 k 折交叉验证。

from sklearn.model_selection import GridSearchCVgd_sr = GridSearchCV(estimator=RFC, param_grid=grid_param, scoring=’accuracy’, cv=5,n_jobs=-1)gd_sr.fit(X_rus, y_rus)  
print(gd_sr.best_params_)

在确定了最佳超参数之后，我用我的机器可以在相对较短的时间内计算的最大数量的树来拟合模型。

Best_RFC=RandomForestClassifier(n_estimators=8000,max_features=
’auto’, max_depth=20,min_samples_split=5, min_samples_leaf=1, bootstrap=True, criterion=’gini’)# fit best model to training dataset
Best_RFC.fit(X_rus, y_rus)

最后，我从测试集中预测了鱼样本的来源，并计算了模型的准确性:

# predict test Y values
ypred=Best_RFC.predict(test_X)from sklearn import metricsprint(“Accuracy:”,metrics.accuracy_score(test_y, ypred))

该模型以 89%的准确率预测了测试集中每个样本的地理来源。我的预测准确率高于对相似鱼类物种起源进行分类的研究。这项工作受到样本量小的限制；我研究的物种不常被捕获。

分类矩阵让我们深入了解模型预测如何与观察到的类别相关联。

随机森林模型预测来自墨西哥湾的鱼类样本比来自东大西洋和西大西洋的样本更准确。这表明墨西哥湾的骨骼形状比其他地区更独特。

这个练习展示了机器学习概念在渔业科学中的价值，以及它们使用我建议的技术预测鱼类样本来源的能力。鉴于样本量较小，我认为我设计的特征提供了很强的预测能力。

我感谢任何反馈和建设性的批评。与该分析相关的代码可在github.com/njermain上找到

使用随机森林预测 H-1B 状态

原文：https://towardsdatascience.com/predicting-h-1b-status-using-random-forest-dc199a6d254c?source=collection_archive---------20-----------------------

当我在做这个项目的时候，我很惊讶我数据科学课程的很多同学都不知道什么是 H-1B 签证。我的简单解释是——H-1B 签证项目允许外国工人在有限的时间内为美国的公司工作。为了让一些从未办理过签证的技术发烧友更感兴趣，我补充道，塞特亚·纳德拉、桑德尔·皮帅和埃隆·马斯克曾持 H-1B 签证移居美国。我希望这个项目能给美国人带来更多的意识，这只是签证的一种。

项目背后的故事

我经常想，是什么让我在其他数据科学家中成为一名优秀的数据科学家。我意识到我的优势在于我的多元化背景，我可以提出正确的问题，并对不同背景、性别和种族的人深表同情。这是我想带到现场的东西，这样数据可以代表每个人，对每个人都有益。

我做这个项目的主要原因是，三年前我持 H-1B 签证从哈萨克斯坦来到西雅图。这是一个极其不确定、漫长且昂贵的过程，我甚至不打算在这里解释所有的步骤。相反，我决定采取行动帮助其他 H-1B 申请人。我用机器学习(ML)预测了签证被接受的概率。利用创建的 ML 模型，我创建了一个应用程序，申请者可以输入信息并立即获得结果。

这不是你的模型，这是你的数据集

第一个数据集

这个项目最困难的部分是选择数据集。我使用的初始数据集存储在 Kaggle 中。

让我们回顾一下我是如何发现我的第一个数据集不起作用的几个步骤。

**我的目标:**确定申请是否被拒绝/接受。

我的第一个数据集包括以下特征:

• 提交的完整日期(包括日/月/年)
• 工资
• 雇主信息
• 行业信息
• 依赖性(丈夫/妻子)
• 违规历史

我收集了 2016 年至 2018 年的 H-1B 申请。

EDA

在开始预测建模之前，我做了一些有趣的 EDA，这与技术行业的数据科学家有关。

Graph 1

图 1 :数据科学的工作岗位主要有三种:数据科学家、数据分析师、数据工程师。根据这三个职位的 H-1B 申请人的工资分布，数据工程师的年收入最高。

Graph 2

图表 2:上面的柱状图代表了 2018 年雇佣数据科学家最多的十家公司。这是一个很好的求职数据科学家指南。

特征工程

我的数据集有许多分类变量，因此在将我的数据拟合到模型中之前，我总结了几个步骤:

• 创建虚拟变量
• 创建工资类别(非常低、低、中、高、非常高)
• 至于工业，我划分了科技和非科技工业
• 分开的年、日和月

分类建模

在探索了逻辑回归、XGBoost 分类和随机森林之后。我的首选是随机森林，原因如下:

• 我的数据集有很多例子和特征
• 与其他模型相比，随机森林通常会给出更准确的结果
• 对于不平衡类学习有一个包-平衡随机森林分类，自动平衡类。
• 我能够通过调整最大深度和估计器的数量来改进我的模型。

AUC = 0.76 for Random Forest

随机森林模型给我的 AUC 值是 0.76。 AUC 的取值范围从 0 到 1。AUC 值越高，模型越准确。

除了高 AUC 分数，即使在使用平衡随机森林分类器后，我对该数据集的精度和召回率也是不平衡的。

我试图改变接受和拒绝类的阈值，但它似乎不起作用。但这是我拿到的最高 f1 分。

在仔细评估我的模型后，我意识到问题出在我的特征上，我需要更多关于每个申请人的信息。

第二数据集

在对一个数据集做了更多的研究后，我发现在 https://www.foreignlaborcert.doleta.gov/performancedata.cfm上，OFLC(外国劳工认证办公室)提供了 H-1B 申请人的完整披露数据。

我添加了以下功能:

• 职称
• 教育专业
• 教育水平
• 国籍
• 多年的经验

我使用与第一个数据集中相同的参数运行了平衡随机森林分类器，结果显著改善:

AUC 评分也有所改善:

汇总

使用第二个数据集，我能够提高我的分数，但是被接受的类的精度和 F1 分数仍然没有被接受的类高。我在 Flask 应用程序中实现了该模型，并将其部署到 Heroku 中。我有一个演示版本，我还在努力使它对人们更有用。你可以通过输入你的信息来尝试一下。

赫尔库 App:https://afternoon-sea-49094.herokuapp.com/

GitHub 库:https://github.com/AisOmar/H1B_predict

预测心脏病死亡率

原文：https://towardsdatascience.com/predicting-heart-disease-mortality-2320fa759f81?source=collection_archive---------29-----------------------

构建可以识别 2019 年高危状态的机器学习模型。

根据疾病控制中心的数据，“美国每年大约有 610，000 人死于心脏病，也就是说每 4 个人中就有一个。”阅读这篇文章的人不可能没有受到这种疾病的影响。今年早些时候，年仅 57 岁的我失去了一位家人。其原因已被很好地记录和理解，但它仍然是美国的主要死亡原因。在这方面，公共政策的改变有可能有助于拯救生命吗？

通过建立一个机器学习模型来预测各州的心脏病死亡率，我们应该能够确定哪些州已经有效地降低了这些死亡率。如果是这样的话，我们就有希望将这些政策原则推广到其他国家。

史料

自 1999 年以来，疾病控制和预防中心每年都会公布心脏死亡率和其他 T2 主要死因的数据。通过跟踪每个州相对于人口的死亡人数，我们可以看到死亡率在过去的 20 年中一直呈下降趋势。

正如我们所见，这一趋势似乎在 2011 年趋于平稳，可能是由于人口老龄化。

心脏病的主要原因

在进行了一些一般性研究后，确定了对心脏病有重大影响的四个主要因素:

• 年龄
• 性别
• 酒精消耗
• 吸烟率

通过对每个预测变量与目标进行简单的线性回归，我们可以开始了解哪些因素具有最显著的影响。

从这里，我们可以看到人口统计学似乎对心脏病有最大的影响，75-79 岁年龄段是主要的预测因素。也许更令人惊讶的是，随着男性在人口中所占比例的增加，心脏病死亡率下降了。这是否意味着男性不易患心脏病？不一定，但它似乎表明女性更有可能死于它。疾病预防控制中心认为意识是这方面的首要因素，

尽管在过去几十年中人们的意识有所提高，但只有大约一半(56%)的女性认识到心脏病是她们的头号杀手。

此外，女性往往会经历更广泛的症状，这可能会导致她们意识不到自己患有心脏病。

另一个令人惊讶的结果是，饮酒似乎与心脏病死亡率没有预期的那么大关联。事实上，对数据的初步审查似乎表明，葡萄酒消费量的增加与死亡率的降低相关。然而，这里有一个警告，葡萄酒消费也被发现与吸烟率呈负相关。由于吸烟确实与心脏病有很强的正相关性，因此有理由认为不是酒，而是较低的吸烟率导致了较低的死亡率。

定义变量

自然，对心脏病的预测需要建立在前期数据的基础上，并有充分的准备时间。换句话说，如果我们想预测 2019 年的结果，需要基于 2018 年及之前的数据。为了分析的目的，使用了 3-5 年的时滞。换句话说，如果预测 2016 年，我们使用 2011-2013 年的数据。

风险水平是基于历史心脏病发病率的总和。如果预测的比率在所有观察到的比率的 33%以下，我们将其标记为低风险。反之，如果是在 33%以上，就是高风险。

构建机器学习模型

有了这些信息，我们现在可以开始构建模型了。总的来说，六种模型架构的 1413 个变体应用于使用 5 重交叉验证的数据的 80%子集。将每个模型应用到测试数据中，我们可以测量准确性并比较不同架构的有效性。

这里，我们可以看到支持向量机模型能够实现 91.7%的 f1 得分。下面是一张显示 2019 年预测的地图，其中九个州被确定为高风险。

最引人注目的是该国中部高风险州的聚集，尽管不能确定这是否是巧合。

值得注意的是，该模型预测，从 2016 年到 2019 年，六个州的风险水平将发生变化，如下图所示。

正如我们所看到的，俄克拉荷马州、密歇根州和宾夕法尼亚州不再被预测为高风险州。然而，田纳西州却反其道而行之，从中等风险到高风险。但是有一种状态特别有趣。

俄克拉何马州

2016 年，俄克拉荷马州观察到的心脏病死亡率将其归入高风险类别，但我们的模型预测它在 2019 年是低风险的。当我第一次看到这个的时候，我认为这是一个错误。一个国家如何能在如此短的时间内如此大幅度地降低其风险水平？事实证明，有强有力的证据支持这一预测。

经过更仔细的检查，很明显这一变化的原因是吸烟率从 2010 年的 25.5%下降到 2015 年的 19%左右，创历史新低。

该州卫生与公众服务部部长特里·克莱恩将这一下降归因于该州的禁烟令。结合认证健康俄克拉荷马计划，激励商业地产做同样的事情，结果成为焦点。

推荐

这项研究的发现为想要降低心脏病死亡率的州指出了两个简单的策略。

• 效仿俄克拉荷马州，在公共场所实施禁烟令。
• 开展针对女性和退休人员的心脏病宣传活动。

城市中的提高认识运动可能最有效，因为妇女往往占人口的较高比例。虽然还有更多的事情可以做，但我希望这一分析能够为今后的一些初步步骤提供一些思路。

资源

这篇文章改编自我几个月前写的一篇关于这个话题的技术文章。如果你想知道更多关于准则和策略的细节，请参考它。

该分析是使用 Jupyter 笔记本和相关库进行的。包含所有代码和数据的库可以在 GitHub 上找到。

从零开始利用深度学习和 Keras 预测医院再入院

原文：https://towardsdatascience.com/predicting-hospital-readmission-with-deep-learning-from-scratch-and-with-keras-309efc0f75fc?source=collection_archive---------27-----------------------

让我们使用深度学习来识别有再次入院风险的患者！

介绍

最近，我看了我的朋友 Eric Ma 关于深度学习基础的视频。为了教授深度学习，他将其分为 3 个关键部分:模型、损失函数和优化例程。在整个教程中，他使用了一个自动微分工具箱。然而，我发现自己做导数是非常令人满意的(至少对于简单的情况)。今天，我想按照 Eric 的方法从头开始构建一个 2 层神经网络，但使用代数导数(来自吴恩达的 Coursera class )，然后使用 Keras(一个深度学习框架)再次实现它。

数据集

对于这个项目，我们将使用与我上一篇关于预测医院入院的文章相同的数据集，该数据集来自 UCI 的糖尿病医院数据集(https://archive . ics . UCI . edu/ml/datasets/diabetes+130-us+hospitals+for+years+1999-2008)。关于这个项目和功能工程的回顾，请参见我之前的帖子https://towards data science . com/predicting-hospital-re-admission-for-patients-with-diabetes-using-scikit-learn-a 2e 359 b 15 f 0

项目定义

预测糖尿病患者是否会在 30 天内再次入院。

特征工程

我们将开始这篇文章，就好像我们已经完成了我上一篇文章的特性工程部分，其中包括创建数字、分类(一键编码)和顺序特性。这些功能方便地保存在以前的笔记本中，并包含在我的 github repo 中。

对于深度学习来说，重要的是填充缺失值并对数据进行归一化。我们将从 scikit-learn 中使用 SimpleImputer 和 StandardScaler 来做这件事。

从头做起

这里，我们将使用 Eric Ma 介绍的结构从头开始构建一个简单的两层神经网络:

• 模型
• 损失函数
• 优化程序

我会尽量坚持吴恩达在他的 Coursera specialization(https://www.coursera.org/specializations/deep-learning)中介绍的符号。

模型

我们将使用的模型是一个双层神经网络，如下所示:

这里我们将有 n_x 个输入变量、n_1 个隐藏节点和一个包括 m 个样本的输出节点。对于这个模型，我们将对隐藏层节点和输出层中的激活函数使用逻辑回归。

这里，我们的激活函数将具有以下形式

为了更有效地计算，我们将使用矢量化记号。在这种表示法中，X 的第一列将是第一个样本的所有特征(注意，这与 Python 中此时加载的内容相反，因此我们需要转置 X 矩阵)。

我们将用于该模型的参数将具有以下维度

为了简单起见，让我们选择 n1 = 64 个节点来近似地将输入变量的数量减半。按照 Eric 的符号，让我们将所有这些参数保存在一个字典中。我们将随机初始化这些，因为设置为 0 将不起作用，因为节点与节点之间的权重都是相同的。

为了计算我们的m示例中 y_hat 的估计值，我们可以使用下面的等式通过模型向前馈送信息(注意维度在花括号中)。

我们可以用下面的前馈函数来编写代码。这里我们将隐藏层的激活函数作为函数的变量。

损失函数

现在我们有了一个给定一些参数计算 y_hat 的方法，我们需要找到“最佳”参数。为了定义“最佳”，我们需要一个成本函数来定义参数有多好。我们用于二元分类的损失函数是:

很明显这是从哪里来的，对吗？

我更愿意看看这是从哪里来的，所以让我们绕一小段路，推导出这个方程。如果我们把我们模型的输出看作给定 x 的 y 的概率，我们可以为一个例子写如下:

可以更巧妙的写成:

如果我们假设所有样本都是独立的，那么看到所有数据的可能性就是个体概率的乘积:

现在我们要做的就是找到最大化这种可能性的参数。这听起来很复杂，因为有产品术语。幸运的是，最大化似然函数的对数也最大化了似然函数(因为对数是单调增加的)。在我们这样做之前，让我们提醒自己关于日志的属性:

应用于我们的似然函数，我们得到:

这已经接近我们的成本函数了！唯一的区别是我们乘以-1，然后除以 m(样本数)。负乘法将它从最大化问题转换到最小化问题。

对于我们的优化程序，我们将需要这个成本函数J的导数。Eric 用 python 包 jax.grad 完成了这个任务

dlogistic_loss = grad(logistic_loss)

但是我想把它明确地写出来，以便更好地理解数学。

为了做这个导数，我们实际上在我们的神经网络中从右向左工作，这个过程被称为反向传播。我们可以用导数的基本原理做到这一点:链式法则！

在我们深入研究这个之前，让我们先来看一些我们将要用到的函数导数

和线性

由于我们的成本是对每个样本求和的，现在让我们从计算中去掉每个样本的符号(和 1/m 乘数)。这里我们将使用吴恩达的简写符号和定义

现在，我们可以对输出层中的参数求导(并适当考虑矩阵数学):

在网络中向后移动并考虑任何激活功能g^[layer](Z[layer])。这里*代表元素级乘法，因为我们对单个样本使用链规则，所以它开始起作用。

正如你在这里看到的，有一个清晰的模式，允许我们将其扩展到任意数量的隐藏层。

现在让我们编写反向传播函数，它将参数、激活函数的导数函数、前馈值和输出值作为输入。

这里d_logistic是逻辑函数相对于z的导数。

此时，最好验证我们的参数和 d_params 对于每个参数集具有相同的形状。

优化程序

我们将使用梯度下降来更新我们的参数。梯度下降基本上通过将梯度的相反方向移动学习量alpha来迭代更新参数。我们可以在 for 循环中运行这个函数，并记录损失。注意 tqdmn 允许我们用一个进度条来查看进度(有点整洁，谢谢 Eric！).

我们可以验证损耗随着迭代而减少:

然后，我们可以计算训练集和验证集的预测值:

使用 scikit-learn 指标，我们可以绘制 ROC 曲线:

这种临时模型的训练非常缓慢。让我们使用 Keras，它有额外的更有效的优化例程，如 Adam。Keras 也非常适合用几行代码构建更复杂的网络。

克拉斯

首先让我们导入一些包

我们需要稍微调整一下 Keras 的输出标签，使每个标签都有一列:

现在我们可以使用序列来构建我们的模型。在这里，我将使用 ReLu 激活函数，而不是逻辑函数，因为 ReLu 往往工作得更好。我还会添加辍学，这是一种正规化的形式，有助于减少过度拟合。

这里，最终的输出层有两个节点(每个标签一个)。使用 softmax 函数将这些输出标准化，以将分数转换为概率。

然后，我们用指定的损失函数和优化器来编译模型。这里‘分类交叉熵’是上面定义的损失函数的公式。这也可以扩展到包括任意数量的结果

现在我们用这个模型

这里有两个输入参数batch_size和epochs包含在拟合中。批量大小表示每次迭代中使用的样本数量。在我们的临时实现中，我们在每次迭代中包含了所有的样本，这需要更多的时间来计算。如果您用较小的批处理运行计算，您就能够更快地迭代。Epoch 被定义为在将整个数据集分成更小的批次后迭代整个数据集的次数。

与 scikit-learn 类似，我们使用 predict proba 获得预测，并只获取第二列:

我在隐藏层、丢弃率和附加层中试验了不同数量的节点。最终模型的验证 AUC = 0.66，如以下 ROC 所示:

不幸的是，这与我们在前一篇文章中训练的所有其他模型具有相同的性能！

结论

在这篇文章中，我们用 Keras 从头开始训练了一个 2 层神经网络。如果你有任何问题，请在下面评论。代码在我的 github 回购上:【https://github.com/andrewwlong/diabetes_readmission_deep