Модель Диксона-Коулза: коррекция низких счетов и временное затухание в футболе

1. Введение: Структурный кризис независимой модели Пуассона

В количественном моделировании футбольных ассоциаций классическая модель независимого распределения Пуассона, предложенная в прикладном спортивном анализе М. Дж. Махером (M. J. Maher) в 1982 году, долгое время служила главным фундаментом для оценки ожидаемых голов и расчета вероятностей исходов матчей. Базовая аксиома Махера утверждает, что количество голов, забитых командой хозяев ($X$) и командой гостей ($Y$) за стандартные 90 минут основного времени, подчиняется взаимно независимым пуассоновским процессам с математическими ожиданиями $\lambda$ и $\mu$ соответственно:

Независимость является модельным допущением, которое нужно проверять на определённом наборе данных. Состояние матча и тактическое взаимодействие могут создавать зависимость, а дополнительная поправка способна подстроиться под шум. В сравнении нужно различать изменение зависимости и изменение интенсивностей голов. Расхождение моделей само по себе не подтверждает ошибку котировки или знак математического ожидания у букмекера.

В 1997 году в журнале Journal of the Royal Statistical Society была опубликована фундаментальная работа британских статистиков Марка Диксона и Стюарта Коулза (Mark J. Dixon & Stuart G. Coles): «Modelling Association Football Scores and Inefficiencies in the Football Betting Market». Модель Диксона-Коулза решила эту проблему с помощью введения матричного поправочного фактора зависимости ${\tau }_{x,y}(\lambda ,\mu ,\rho )$, избирательно перераспределяющего вероятности низких счетов с сохранением маргинальных распределений Пуассона, а также экспоненциального временного затухания формы команд.

В статье выведена поправка четырёх ячеек с сохранением пуассоновских маргинальных распределений, задано взвешенное по времени правдоподобие и приведён пример Python с ограничениями. В числовом сравнении используются явные фиксированные параметры. Подгонку и вероятности можно воспроизводить отдельно от утверждений о текущих коэффициентах букмекеров и прибыльности.

2. Двумерный корректирующий фактор Диксона-Коулза ${\tau }_{x,y}$

Rho является параметром зависимости низких счетов, а не коэффициентом корреляции Пирсона между числами голов. Допустимый диапазон зависит от обеих интенсивностей. Реализация правдоподобия ниже требует строго положительных факторов tau. При lambda=10, mu=2 и rho=-0.35 получается tau(0,1)=-2.5. Набор параметров нужно отклонять целиком, даже при расчёте высокого счёта. Обнуление этой ячейки и перенормировка дали бы другое распределение.

Чтобы ликвидировать расхождение на низких счетах, не разрушая при этом пуассоновскую структуру распределения, Диксон и Коулз ввели модифицированную функцию совместной плотности вероятности:

Корректирующий фактор ${\tau }_{x,y}(\lambda ,\mu ,\rho )$ представляет собой кусочно-дискретный скалярный множитель, воздействующий на счетные комбинации $\{0,1\}\times \{0,1\}$ и зависящий от параметра взаимосвязи $\rho$:

Математическое доказательство сохранения полной вероятности

Для корректности статистической модели модификация обязана строго удовлетворять условию нормировки: сумма всех элементов матрицы совместных вероятностей должна быть равна единице: ${\sum }_{x=0}^{\infty }{\sum }_{y=0}^{\infty }P(X=x,Y=y)=1$.

Разобьем двойную сумму на исходный пуассоновский базис и сумму возмущений по четырем модифицированным ячейкам:

Распишем поправку $\mathrm{\Delta }$ по состояниям $(\mathrm{0,0}),(\mathrm{0,1}),(\mathrm{1,0}),(\mathrm{1,1})$:

Подставляя вероятности $P_0(\lambda )=e^{-\lambda }$, $P_1(\lambda )=\lambda e^{-\lambda }$, $P_0(\mu )=e^{-\mu }$ и $P_1(\mu )=\mu e^{-\mu }$, выносим общий множитель $e^{-\lambda -\mu }$ за скобки:

Сумма возмущений строго равна нулю! Это означает, что фактор Диксона-Коулза лишь перераспределяет вероятностную массу между четырьмя базовыми состояниями, не нарушая маргинальных пуассоновских интегралов и сохраняя тотальную вероятность на уровне ровно $1.00$ ($100\%$).

Физический смысл отрицательного параметра $\rho$

В примере подгонки выбрана область неположительного rho вместе с ограничениями положительности, зависящими от интенсивностей. Это выбор модели, а не универсальный знак, установленный по футбольным данным. Всё семейство допускает и положительный rho: при lambda=1.65, mu=1.15 и rho=0.10 четыре фактора равны 0.81025, 1.165, 1.115 и 0.9, то есть положительны.

В допустимой области отрицательный rho даёт следующие алгебраические изменения относительно независимости при тех же интенсивностях:

• Для счета 0:0: ${\tau }_{\mathrm{0,0}}=1-\lambda \mu \rho$. Так как $\rho <0$, произведение $-\lambda \mu \rho >0$, откуда ${\tau }_{\mathrm{0,0}}>1$. Вероятность безголевой ничьей возрастает относительно наивного Пуассона.
• Для счета 1:1: ${\tau }_{\mathrm{1,1}}=1-\rho$. Так как $\rho <0$, значение $1-\rho =1+|\rho |>1$. Вероятность результативной ничьей 1:1 также увеличивается.
• Для счетов 1:0 и 0:1: ${\tau }_{\mathrm{0,1}}=1+\lambda \rho <1$ и ${\tau }_{\mathrm{1,0}}=1+\mu \rho <1$. Вероятности минимальных побед 1:0 понижаются, идеально балансируя прирост ничьих.

3. Параметризация сил команд: атака, оборона и фактор поля

В модели Диксона-Коулза ожидаемые голы хозяев ${\lambda }_{i,j}$ и гостей ${\mu }_{i,j}$ выражаются через скрытые латентные рейтинги силы клубов:

Где:

• ${\alpha }_i>0$: Рейтинг атакующей мощи команды $i$. При ${\alpha }_i=1.30$ команда генерирует на 30% больше опасных моментов и голов, чем среднестатистический клуб лиги.
• ${\beta }_j>0$: Коэффициент оборонительной уязвимости команды $j$. При ${\beta }_j=0.80$ клуб пропускает на 20% меньше голов относительно среднего уровня чемпионата.
• $\gamma >0$: Лига-фактор преимущества своего поля. В высших дивизионах Европы эмпирическое значение $\gamma$ колеблется в коридоре от $1.20$ до $1.35$.

Ограничение идентифицируемости параметров

Если умножить все параметры атаки ${\alpha }_k$ на произвольную константу $c>0$, а все параметры обороны ${\beta }_k$ разделить на $c$, величины ожидаемых голов ${\lambda }_{i,j}=(c{\alpha }_i)({\beta }_j/c)\gamma ={\alpha }_i{\beta }_j\gamma$ останутся абсолютно неизменными. Чтобы ликвидировать эту масштабную неопределенность, Диксон и Коулз зафиксировали строгое ограничение на среднее арифметическое атакующих параметров лиги из $N$ команд:

Средняя атака, равная единице, устраняет симметрию масштаба между атакой и защитой. Это не доказывает единственность глобального максимума правдоподобия, достаточность информации в сыгранных парах и предсказательную точность. Рейтинги относительны к подогнанной модели и её данным. Одна нормализация не делает рейтинги разных турниров напрямую сопоставимыми.

4. Экспоненциальное затухание формы: параметр памяти $\xi$

Игры двухлетней давности не должны иметь равный вес с матчами, сыгранными на прошлой неделе. Составы команд обновляются в трансферные окна, тренеры меняют схемы, ключевые игроки выбывают из-за травм.

Для непрерывного дисконтирования давности Диксон и Коулз разработали функцию экспоненциального затухания весов $w_k(\xi )$ для каждого исторического матча $k$, сыгранного в момент $t_k$, относительно текущего дня $t$:

Где $(t-t_k)$ измеряется в днях, а $\xi >0$ — непрерывный гиперпараметр затухания памяти модели.

5. Оптимизация методом максимального правдоподобия (MLE)

Для выборки из $M$ завершившихся матчей с результатами $(x_k,y_k)$ и датами $t_k$ логарифмическая функция правдоподобия с весовыми коэффициентами формулируется следующим образом:

Опуская независимые факториалы, минимизируется взвешенное отрицательное логарифмическое правдоподобие:

Пример использует SLSQP для равенства масштаба атаки и нелинейных неравенств tau. Перед сохранением параметров проверяются успех оптимизации и допустимость результата. Успешное численное завершение не подтверждает глобальный оптимум или качество прогноза будущих матчей.

6. Числовое сравнение при фиксированных интенсивностях

В расчёте заданы постоянные $\lambda =1.65$, $\mu =1.15$ и $\rho =-0.12$. Вероятности суммируются до 45 голов. Коэффициенты равны обратным модельным вероятностям. Для воспроизведения примера не нужны наблюдаемые котировки закрытия. Различие двух моделей само по себе не подтверждает рыночный перевес.

Тотал меньше 2.5 содержит все четыре изменённые ячейки, поправки которых сокращаются. Вероятность равна 46.945368% в обеих моделях, коэффициент около 2.1301. Тотал больше 2.5 тоже не меняется. Зависимость может изменить тотал, включающий лишь часть этих ячеек, например меньше 1.5, а также ничьи и «обе забьют». Повторная оценка интенсивностей способна изменить тотал 2.5 через новые интенсивности. Это отдельный эффект, отличный от применения rho при неизменных интенсивностях.

7. Промышленный код на Python: DixonColesModel

Пример ограничивает все четыре фактора низких счетов для каждой упорядоченной пары известных команд, фиксирует среднюю атаку на единице и проверяет успех и допустимость результата оптимизации. Неверные вероятности отклоняются, а не обнуляются. По умолчанию сетка использует квантили Пуассона с отброшенной массой примерно меньше 2e-12, сообщает эту массу и сохраняет исходные вероятности без перенормировки. Xi задан как условный дневной коэффициент затухания; его выбор требует отдельной проверки.

import numpy as np
from scipy.optimize import minimize
from scipy.stats import poisson


class DixonColesModel:
    def __init__(self, xi=0.0065):
        # Illustrative daily decay rate, not a validated universal optimum.
        if not np.isfinite(xi) or xi < 0:
            raise ValueError("xi must be a finite nonnegative daily rate")
        self.xi = xi
        self.teams = []
        self.team_indices = {}
        self.params = None
        self.gamma = 1.25
        self.rho = -0.11

    @staticmethod
    def tau(x, y, lambda_, mu, rho):
        # Check the whole joint model, including unobserved low-score cells.
        if not np.all(np.isfinite([lambda_, mu, rho])) or min(lambda_, mu) <= 0:
            raise ValueError("Intensities must be finite and positive")
        factors = [1 - lambda_ * mu * rho, 1 + lambda_ * rho,
                   1 + mu * rho, 1 - rho]
        if min(factors) <= 0:
            raise ValueError("rho is outside the strictly positive tau domain")
        return {(0, 0): factors[0], (0, 1): factors[1],
                (1, 0): factors[2], (1, 1): factors[3]}.get((x, y), 1.0)

    def bivariate_prob(self, x, y, lambda_, mu, rho):
        if not np.all(np.isfinite([x, y])) or min(x, y) < 0 or int(x) != x or int(y) != y:
            raise ValueError("Goals must be nonnegative integers")
        factor = self.tau(x, y, lambda_, mu, rho)
        return factor * poisson.pmf(x, lambda_) * poisson.pmf(y, mu)

    def _unpack_params(self, params, n_teams):
        return (params[:n_teams], params[n_teams:2*n_teams],
                params[2*n_teams], params[2*n_teams + 1])

    def fit(self, matches):
        self.params = None
        if not matches:
            raise ValueError("At least one match is required")
        for m in matches:
            goals = [m['home_goals'], m['away_goals']]
            age = m.get('days_ago', 0.0)
            if (m['home'] == m['away'] or not np.all(np.isfinite(goals + [age]))
                    or min(goals) < 0 or any(int(v) != v for v in goals) or age < 0):
                raise ValueError("Invalid teams, goals or match age")
        self.teams = sorted({m['home'] for m in matches} | {m['away'] for m in matches})
        self.team_indices = {team: i for i, team in enumerate(self.teams)}
        n_teams = len(self.teams)
        init = np.concatenate([np.ones(n_teams), np.ones(n_teams), [1.25, -0.10]])

        def tau_constraints(params):
            alpha, beta, gamma, rho = self._unpack_params(params, n_teams)
            values = []
            # Every ordered pairing of known teams must define a valid model.
            for h in range(n_teams):
                for a in range(n_teams):
                    if h == a:
                        continue
                    lam, mu = alpha[h] * beta[a] * gamma, alpha[a] * beta[h]
                    values.extend([1 - lam * mu * rho, 1 + lam * rho,
                                   1 + mu * rho, 1 - rho])
            return np.asarray(values) - 1e-8

        def loss_fn(params):
            alpha, beta, gamma, rho = self._unpack_params(params, n_teams)
            nll = 0.0
            for m in matches:
                h, a = self.team_indices[m['home']], self.team_indices[m['away']]
                lam, mu = alpha[h] * beta[a] * gamma, alpha[a] * beta[h]
                try:
                    factor = self.tau(m['home_goals'], m['away_goals'], lam, mu, rho)
                except ValueError:
                    return 1e12  # Reject an infeasible optimizer trial, never clip a probability.
                weight = np.exp(-self.xi * m.get('days_ago', 0.0))
                log_prob = (np.log(factor) + poisson.logpmf(m['home_goals'], lam)
                            + poisson.logpmf(m['away_goals'], mu))
                nll -= weight * log_prob
            return float(nll)

        constraints = [
            {'type': 'ineq', 'fun': tau_constraints},
            {'type': 'eq', 'fun': lambda p: np.mean(p[:n_teams]) - 1.0},
        ]
        bounds = [(0.05, 3.5)] * (2*n_teams) + [(1.0, 1.8), (-0.35, 0.0)]
        res = minimize(loss_fn, init, method='SLSQP', bounds=bounds,
                       constraints=constraints, options={'maxiter': 1000, 'ftol': 1e-9})
        if (not res.success or not np.all(np.isfinite(res.x)) or not np.isfinite(res.fun)
                or np.min(tau_constraints(res.x)) < -1e-9
                or abs(np.mean(res.x[:n_teams]) - 1) > 1e-8):
            raise RuntimeError(f"No valid converged fit: {res.message}")
        self.params = res.x
        _, _, self.gamma, self.rho = self._unpack_params(res.x, n_teams)
        return res

    def predict_score_matrix(self, home_team, away_team, max_goals=None):
        if self.params is None:
            raise ValueError("Fit the model first")
        alpha, beta, gamma, rho = self._unpack_params(self.params, len(self.teams))
        h, a = self.team_indices[home_team], self.team_indices[away_team]
        lam, mu = alpha[h] * beta[a] * gamma, alpha[a] * beta[h]
        self.tau(0, 0, lam, mu, rho)
        if max_goals is None:
            max_goals = int(max(poisson.ppf(1 - 1e-12, lam), poisson.ppf(1 - 1e-12, mu)))
        if not isinstance(max_goals, (int, np.integer)) or max_goals < 1:
            raise ValueError("max_goals must be an integer of at least one")
        goals = np.arange(max_goals + 1)
        matrix = np.outer(poisson.pmf(goals, lam), poisson.pmf(goals, mu))
        for x, y in [(0, 0), (0, 1), (1, 0), (1, 1)]:
            matrix[x, y] *= self.tau(x, y, lam, mu, rho)
        # Keep the original probability mass; disclose the omitted grid tail.
        self.omitted_mass = float(max(0.0, 1 - matrix.sum()))
        return matrix, lam, mu

    def predict_1x2(self, home_team, away_team):
        matrix, lam, mu = self.predict_score_matrix(home_team, away_team)
        home, draw, away = np.tril(matrix, -1).sum(), np.trace(matrix), np.triu(matrix, 1).sum()
        return {'home_win': float(home), 'draw': float(draw), 'away_win': float(away),
                'fair_odds_home': float(1/home), 'fair_odds_draw': float(1/draw),
                'fair_odds_away': float(1/away), 'lambda': float(lam), 'mu': float(mu),
                'omitted_probability_mass': self.omitted_mass}

8. Современные расширения: копулы и двумерный Пуассон Карлиса-Нцофраса

Хотя спецификация Диксона-Коулза остается общепризнанным стандартом индустрии для учета аномалий низких счетов, количественные исследователи разработали альтернативные методы моделирования взаимной зависимости голов:

Двумерный Пуассон с общим фактором шока (Karlis & Ntzoufras, 2003)

В модели Карлиса и Нцофраса голы команд моделируются через три независимые пуассоновские случайные величины: $X=X^{\ast }+Z$ и $Y=Y^{\ast }+Z$, где $Z\sim \text{Pois}({\lambda }_3)$ представляет собой скрытый фактор «общего шока» (общий темп матча, судейская строгость, погодные катаклизмы). Ковариация в такой системе строго неотрицательна: $\text{Cov}(X,Y)={\lambda }_3\ge 0$. Эта спецификация превосходно описывает открытые результативные матчи с высокими тоталами, но не способна адекватно смоделировать отрицательную корреляцию и защитную блокировку при счетах 0:0 и 1:1 в прагматичном позиционном футболе.

Копулы Архимеда для дискретных распределений

Современные алгоритмические фонды комбинируют непрерывные копулы (копулу Франка или Клейтона) с дискретными отрицательными биномиальными маргиналами, учитывающими эффект сверхдисперсии (overdispersion). Копулы позволяют раздельно калибровать хвостовую зависимость для экстремальных счетов. Тем не менее, по вычислительной скорости, интерпретируемости и мгновенной скорости пересчета в лайве фактор ${\tau }_{x,y}$ Диксона-Коулза по-прежнему остается золотым стандартом для маркетмейкеров и профессиональных синдикатов.

9. Аналитическое заключение и практические выводы

Математический анализ показывает: отказ от наивного предположения о независимости забитых голов — обязательный шаг для выхода на устойчивый плюс (+EV). Коррекция Диксона-Коулза ликвидирует структурное занижение вероятности ничьих, выравнивает цены на тоталы меньше и формирует математический перевес над рекреационными линиями букмекеров.