用Python预测足球：逻辑斯蒂克模型前篇

浅夏110

这个系列的整体架构如下

第1部分：逻辑斯迪克模型。主要介绍logistic模型和ordered logistic模型，ordered logistic可以用来预测足球的胜负平概率。

第2部分：泊松模型。主要介绍最基础的poisson模型，poisson模型可以预测足球比赛的进球数的概率，从而依据进球数预测胜负平的概率。

第3部分：零膨胀泊松模型及联立泊松模型。由于足球比赛中存在大量的某一方0进球比赛（比如英超17-18赛季的水晶宫，截止10月8日一球未进），传统的泊松模型可能误差较大，所以我们需要零膨胀泊松模型。更近一步，我们认为主队进球和客队进球并不是相互独立，这时我们需要联立泊松模型，而EM算法可以很好的估计联立泊松模型。

第4部分：状态空间模型和时间序列。足球比赛本质上是一个离散时间序列，我们可以在逻辑斯迪克模型或泊松模型上添加马尔科夫过程，通过卡尔曼滤波或mcmc估计进化方差。

第5部分：岭回归，lasso和弹性网。我们可以为我们的模型添加很多特征变量，例如是否参加欧冠，是否夺冠关键战，当特征变量太多的时候我们会需要进行特征选择，岭回归，lasso和弹性网是我们的好帮手，如何开发快速的坐标下降算法，是个挑战。

第6部分：支持向量机。终于来到机器学习的世界，支持向量机特别擅长处理小样本高维空间的特征选择，对于足球比赛这种变量狂多样本又较少的预测，效果卓越。

第7部分：神经网络。笔者不确定神经网络的效果好于glm和svm，反正看着办吧。

现在我们进入正题。

logistic模型

logistic模型大概是使用最频繁的二分类模型。对于体育比赛来说，logistic模型又可以称之为TBL(bradley-terry-lucy)模型（某种程度上），该因为BTL的形式为 ，其中参数 可以看做为体育比赛对正双方的实力。对于了解在线游戏配对系统的人来说，BTL和elo rating很像，那是因为elo rating用的对数底是10，而BTL是自然对数，BTL更多的是用来进行参数估计。

既然BTL模型只是特殊的logistic模型，我们可以用普通logistic进行求解，只需要稍微改下输入数据的形式。我们先看看如何求解logistic模型。

from __future__ import print_function, divisionimport numpy as npimport warningsdef z(threshold, x, b):    # z函数    return threshold + np.dot(x, b)def logistic(z_val):    # logistic函数或sigmoid函数    z_val[z_val < -700] = -700    return 1 / (1 + np.exp(-z_val))
1 r: v! ^  X7 X

上面的代码是最基础的logistic模型，其中z函数中threshold是阈值，x是数据，b是权重，logistic函数中的z_val是z函数的值。

逻辑斯迪克模型的参数估计大致有三种，极大似然估计，MCMC估计和迭代加权最小二乘估计。迭代加权最小二乘估计（irls）是目前最稳健也是相对较快，实现起来也很简单的算法。想了解迭代加权最小二乘估计，可以查看维基百科Iteratively reweighted least squares。对于logistic模型来说，irls的误差为 ，权重为 ，其中 是响应变量， 是logistic函数值。irls的代码如下

def irls(y, x, max_iter=1000):    # 待估计参数初值    params = np.zeros(x.shape[1] + 1)    # 加入了阈值哑变量的新的数据值    new_x = np.hstack((np.ones_like(y), x))    for i in range(max_iter):        z_val = z(params[:1], x, params[1:])        logistic_val = logistic(z_val)        # 残差        e = (y[:, 0] - logistic_val) / (logistic_val * (1 - logistic_val))        # 权重        _w = logistic_val * (1 - logistic_val)        w = np.diag(_w ** 0.5)        wa = np.dot(w, new_x)        temp1 = np.dot(wa.transpose(), w)        temp2 = np.linalg.inv(np.dot(wa.transpose(), wa))        params_temp = np.dot(np.dot(temp2, temp1), (z_val + e))        if np.max(np.abs(params_temp - params)) <= 1e-5:            return params_temp        params = params_temp    warnings.warn('no coverage')    return params
) D$ ]: r" j: K5 W' d5 f

logistic模型的测试代码如下

x = np.random.uniform(-3, 3, (1000, 10))b = np.random.uniform(0, 1, (10, 1))threshold = np.random.uniform(-2, 2)z_val = z(threshold, x, b)logistic_val = logistic(z_val)y = np.random.binomial(1, logistic_val)est_params = irls(y, x)est_threshold = est_params[0]est_b = est_params[1:]

BTL模型虽然是使用最广泛的体育比赛统计模型，但其只能处理二分类问题，也即是只有胜负的体育比赛，例如NBA，网球，围棋等，而足球比赛存在平局，所以BTL模型并不适用于足球赛事。有人可能会想到用多分类的logistic模型处理足球赛事，但是足球比赛的胜平负是有顺序关系的，胜的等级明显高于平，而平又高于负，多分类模型不适用这种等级变量。所幸，有专门的logistic模型处理等级变量，那就是ordered logistic模型。

ordered logistic模型

ordered logistic模型又可以称之为proportional odds模型（统计学很多名词其实表达的意思一样），是专门处理等级数据的logistic模型。其数学形式为 ，而 是logistic模型。想详细了解ordered logistic，请看维基百科Ordered logit

对于ordered logistic模型，直接用胜负平这种数据是不合适的，我们需要转换数据

例如

['胜', '平', '胜', '负']
' Z' y7 P# w2 G$ i9 ^% J9 y

转换为

[[1, 0, 0],[0, 1, 0],[1, 0, 0],[0, 0, 1]]

所以本质上，我们还是将等级数据转换成了二分类问题

由于之前用的是irls算法，我们这次对ordered logistic用拟牛顿（BFGS）和牛顿迭代进行参数估计。基础代码如下

# coding=utf-8from __future__ import print_function, divisionimport numpy as npdef z(thresholds, x, b):    # z函数    _z_val = thresholds + np.dot(x, b.transpose())    _z_val[_z_val < -700] = -700    return _z_valdef logistic(z_val):    # logistic函数    return 1 / (1 + np.exp(-z_val))def pq(p_val):    # p * (1 -p)值    return p_val * (1 - p_val)
+ V% y, V( j7 u+ C/ x3 u& p

拟牛顿迭代（BFGS）算法

相比牛顿迭代，BFGS算法的优势是不用求解复杂的黑塞矩阵，劣势是迭代次数较多，必须要进行步长的一维搜索，如果不用scipy写好的C代码进行步长一维搜索，那速度还真是比牛顿迭代慢不少。想了解BFGS算法，还是请看维基百科Broyden-Fletcher-Goldfarb-Shanno algorithm

实现BFGS算法，我们首先创建一维搜索所需的对数似然函数（我们的一维搜索采取暴力搜索方法，考虑到python的龟速for循环，用基于numpy广播机制的暴力搜索不一定慢），我们的步长设定是0.001到0.1，每一步递增0.001。所以对数似然函数和暴力搜索函数如下

def loglik(y, thresholds, x, b):    # 用于暴力搜索的对数似然函数    loglik_val = 0    rep_len = y.shape[1]    p_val = {}    for i in range(thresholds.shape[1]):        z_val = z(thresholds[:, i], x, b)        p_val = logistic(z_val)    for i in range(rep_len):        p_pre = np.ones((y.shape[0], 99)) if i == 0 else p_val[i - 1]        p = np.zeros((y.shape[0], 99)) if i == rep_len - 1 else p_val        # 防止对数运算溢出        p_pre[p_pre < p] = p[p_pre < p]        loglik_val += np.dot(y[:, i], np.log(p_pre - p + 1e-200))    return loglik_valdef line_search(y, thresholds, x, b, p):    # 暴力一维搜索    steps = np.arange(0.001, 0.1, step=0.001)    loglik_val = loglik(y, thresholds - steps[:, np.newaxis] * p[:2], x, b - steps[:, np.newaxis] * p[2:])    return steps[loglik_val.argmax()]
/ z4 F! J  k9 q/ b7 V" i  k

BFGS虽然不需要求解二阶导数了，一阶雅克比矩阵还是需要的。雅克比矩阵如下

def dloglik(y, x, len_threshold, rep_len, p_val, pq_val):    # 一阶导数    slop_dloglik_val = 0    thresholds_dloglik_val = np.zeros(len_threshold)    for i in range(rep_len):        p_pre, dp_pre = (np.ones(y.shape[0]), np.zeros(y.shape[0])) if i == 0 else (p_val[:, i - 1], pq_val[:, i - 1])        p, dp = (np.zeros(y.shape[0]), np.zeros(y.shape[0])) if i == rep_len - 1 else (p_val[:, i], pq_val[:, i])        slop_dloglik_val += np.sum((x * (1 - p_pre[:, np.newaxis] - p[:, np.newaxis]))[y[:, i] != 0, :], axis=0)        if i < rep_len - 1:            thresholds_dloglik_val += np.sum((dp / (p - p_pre + 1e-200))[y[:, i] != 0])        if i > 0:            thresholds_dloglik_val[i - 1] += np.sum((dp_pre / (p_pre - p + 1e-200))[y[:, i] != 0])    return np.concatenate((thresholds_dloglik_val, slop_dloglik_val))
2 l( [1 J! q# M( I: `

接下来便是构造模拟黑塞矩阵，并且求解参数啦

def bfgs(y, thresholds, x, b, max_iter=2000):    # BFGS拟牛顿算法    len_threshold = len(thresholds)    rep_len = y.shape[1]    d = np.eye(b.shape[1] + len(thresholds))    dloglik_val = _get_dloglik(b, len_threshold, rep_len, thresholds, x, y)    for i in range(max_iter):        _delta = -np.dot(d, dloglik_val)        step = line_search(y, thresholds, x, b, _delta)        delta = step * _delta        temp_b = b - delta[len_threshold:]        temp_thresholds = thresholds - delta[:len_threshold]        if np.max(np.abs(delta)) < 1e-5:            return temp_thresholds, temp_b        b = temp_b        thresholds = temp_thresholds        temp_dloglik_val = _get_dloglik(b, len_threshold, rep_len, thresholds, x, y)        dloglik_delta = temp_dloglik_val - dloglik_val        params = np.concatenate((thresholds[np.newaxis, :], b), axis=1)        d = get_fake_hess(dloglik_delta, d, delta, params)        dloglik_val = temp_dloglik_val    return thresholds, bdef _get_dloglik(b, len_threshold, rep_len, thresholds, x, y):    z_val = z(thresholds, x, b)    p_val = logistic(z_val)    pq_val = pq(p_val)    dloglik_val = dloglik(y, x, len_threshold, rep_len, p_val, pq_val)    return dloglik_valdef get_fake_hess(dloglik_delta, d, delta, params):    # bfgs的模拟黑塞矩阵    eye = np.eye(params.shape[1])    temp0 = np.dot(dloglik_delta.transpose(), delta)    temp1 = eye - np.dot(delta[:, np.newaxis], dloglik_delta[:, np.newaxis].transpose()) / temp0    temp2 = eye - np.dot(dloglik_delta[:, np.newaxis], delta[:, np.newaxis].transpose()) / temp0    d = np.dot(temp1, np.dot(d, temp2)) + np.dot(delta[:, np.newaxis], delta[:, np.newaxis].transpose()) / temp0    return d
* N4 m, T) d" N" j+ a6 s7 C

以上代码中，像y代表响应变量，x代表自变量，b代表权重或斜率，thresholds代表阈值，都是通用表示方法。测试代码如下

x = np.random.normal(0, 2, (1000, 10))b = np.random.uniform(1, 2, (1, 10))z1_val = z(-2, x, b)sigmoid1_val = logistic(z1_val)z2_val = z(1, x, b)sigmoid2_val = logistic(z2_val)p1 = sigmoid1_valp2 = sigmoid2_val - sigmoid1_valp3 = 1 - sigmoid2_valp = np.concatenate((p3, p2, p1), axis=1)y = np.zeros((1000, 3))for i, _ in enumerate(p):    y = np.random.multinomial(1, _)res = bfgs(y, np.array([1, 0]), x, np.zeros((1, 10)))

牛顿迭代

纯python代码的BFGS算法实在太慢，我们再来编写能快速收敛的牛顿迭代算法。牛顿迭代需要求解黑塞矩阵。黑塞矩阵代码如下

def ddloglik(y, x, len_threshold, len_b, rep_len, p_val, pq_val, b_idx):    # 二阶导数得期望值（fisher score）    ddloglik_val = np.zeros((len_threshold + len_b, len_threshold + len_b))    for i in range(rep_len):        p, dp= (np.zeros(y.shape[0]), np.zeros(y.shape[0])) if i == rep_len - 1 else (p_val[:, i], pq_val[:, i])        p_pre, dp_pre= (np.ones(y.shape[0]), np.zeros(y.shape[0])) if i == 0 else (p_val[:, i - 1], pq_val[:, i - 1])        b_base_ddloglik = -(dp_pre - dp) * (1 - p_pre - p)        if i < rep_len - 1:            ddloglik_val[len_threshold:, i] += np.sum(x * ((1 - p_pre - p) * dp)[:, np.newaxis], axis=0)            ddloglik_val[i, i] += -np.sum(dp ** 2 / (p_pre - p + 1e-200))        if i > 0:            ddloglik_val[i - 1, i - 1] += -np.sum(dp_pre ** 2 / (p_pre - p + 1e-200))            ddloglik_val[len_threshold:, i - 1] += -np.sum(x * (dp_pre * (1 - p_pre - p))[:, np.newaxis], axis=0)        if 0 < i < rep_len - 1:            ddloglik_val[i, i - 1] = np.sum(dp * dp_pre / (p_pre - p + 1e-200))        for j in range(len_b):            up_idx = np.arange(j, len_b)            down_idx = b_idx[:len_b - j]            new_left = x[:, down_idx] * x[:, up_idx] * b_base_ddloglik[:, np.newaxis]            temp_ddloglik_val = np.sum(new_left, axis=0)            ddloglik_val[len_threshold:, len_threshold:][[up_idx], [down_idx]] += temp_ddloglik_val    ddloglik_val += ddloglik_val.transpose() - np.diag(ddloglik_val.diagonal())    return ddloglik_val

看到这密密麻麻的代码，笔者也是头疼，但是黑塞矩阵就是这么复杂，接下来我们呈现步长横为1的牛顿迭代代码，实际上，根据研究，不计算步长的牛顿迭代不保证100%收敛全局最大值，但是实际应用效果还不错。代码如下

def newton(y, thresholds, x, b, max_iter=2000):    # 牛顿迭代    len_b = b.shape[1]    len_threshold = len(thresholds)    rep_len = y.shape[1]    b_idx = np.arange(len_b)    for i in range(max_iter):        z_val = z(thresholds, x, b)        p_val = logistic(z_val)        pq_val = pq(p_val)        dloglik_val = dloglik(y, x, len_threshold, rep_len, p_val, pq_val)        ddloglik_val = ddloglik(y, x, len_threshold, len_b, rep_len, p_val, pq_val, b_idx)        delta = np.linalg.solve(ddloglik_val, dloglik_val)        temp_b = b - delta[len_threshold:]        temp_thresholds = thresholds - delta[:len_threshold]        if np.max(np.abs(delta)) < 1e-5:            return temp_thresholds, temp_b        b = temp_b        thresholds = temp_thresholds    return thresholds, b

测试代码方面，只需要把BFGS测试代码的最后一行换成

res = newton(y, np.array([1, 0]), x, np.zeros((1, 10)))
% J" `' Y/ e1 K; D: K

我们的快速ordered logistic参数估代码大功告成

% L) s2 R9 U' b  P/ V