数学建模社区-数学中国
标题: 用Python预测足球：逻辑斯蒂克模型前篇 [打印本页]
作者: 浅夏110 时间: 2020-5-16 11:18
标题: 用Python预测足球：逻辑斯蒂克模型前篇
这个系列的整体架构如下
第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)) b# u; ?) b. g. H5 Z
上面的代码是最基础的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 params5 x8 t5 k6 w: @% l$ h
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:]
# G7 {7 W" C) G$ T7 K* y! sBTL模型虽然是使用最广泛的体育比赛统计模型，但其只能处理二分类问题，也即是只有胜负的体育比赛，例如NBA，网球，围棋等，而足球比赛存在平局，所以BTL模型并不适用于足球赛事。有人可能会想到用多分类的logistic模型处理足球赛事，但是足球比赛的胜平负是有顺序关系的，胜的等级明显高于平，而平又高于负，多分类模型不适用这种等级变量。所幸，有专门的logistic模型处理等级变量，那就是ordered logistic模型。
ordered logistic模型ordered logistic模型又可以称之为proportional odds模型（统计学很多名词其实表达的意思一样），是专门处理等级数据的logistic模型。其数学形式为 
，而 
是logistic模型。想详细了解ordered logistic，请看维基百科Ordered logit 对于ordered logistic模型，直接用胜负平这种数据是不合适的，我们需要转换数据
例如
['胜', '平', '胜', '负']
+ S0 g4 Z) P% I4 u( o0 e转换为
[[1, 0, 0],[0, 1, 0],[1, 0, 0],[0, 0, 1]]
: S$ j# P) U6 S; x! t# O! T所以本质上，我们还是将等级数据转换成了二分类问题
由于之前用的是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)0 M) C0 b) ]9 r0 [
拟牛顿迭代（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()]
( r7 h- p! l) o' [( K- z' H; nBFGS虽然不需要求解二阶导数了，一阶雅克比矩阵还是需要的。雅克比矩阵如下
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))9 r% y2 n0 n# x& E
接下来便是构造模拟黑塞矩阵，并且求解参数啦
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/ ?7 A0 e+ y3 h. g5 L2 r- ~
以上代码中，像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)))
! i6 X+ a; x+ o4 }2 H" T, ^牛顿迭代
纯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% ~+ T# K4 A) g! z+ J! O5 \ y2 |
看到这密密麻麻的代码，笔者也是头疼，但是黑塞矩阵就是这么复杂，接下来我们呈现步长横为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
& b3 [8 \; R5 ?/ H0 t( o测试代码方面，只需要把BFGS测试代码的最后一行换成
res = newton(y, np.array([1, 0]), x, np.zeros((1, 10)))
9 w3 T. d* k: h$ j4 `( M. d) H我们的快速ordered logistic参数估代码大功告成
4 |5 f% Z1 k% i" a7 A
7 w; K& ^; `" {$ m& D1 w" ~8 \
| 欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) |
Powered by Discuz! X2.5 |