"""Fit your own xG: reproduction of Lab 1. Teaching artifact by Ahmed Beddi, on StatsBomb open data (2018 FIFA World Cup). Run next to wc2018_shots.csv, wc2018_own_goals.csv and wc2018_matches.csv."""

# # Fit your own xG — a reproduction of Lab 1 (Soccer analytics: maths and statistics)
# 
# **Teaching artifact by Ahmed Beddi, on public data.** Shots from the 2018 FIFA World Cup, StatsBomb open data (competition 43, season 3), exported to `wc2018_shots.csv`. Data provided by StatsBomb under its open-data user agreement (https://github.com/statsbomb/open-data).
# 
# This notebook reproduces the headline numbers of the course's Lab 1 (Sessions 3–5): a two-feature logistic expected-goals (xG) model built from shot geometry, judged against StatsBomb's own xG, and the penalty-spot gap that shows where geometry stops working. The interactive xG explorer on the course page uses the same formulas and the coefficients fitted here.

import numpy as np, pandas as pd, statsmodels.api as sm
import matplotlib.pyplot as plt
shots = pd.read_csv('wc2018_shots.csv')
own = pd.read_csv('wc2018_own_goals.csv')
matches = pd.read_csv('wc2018_matches.csv')
print(len(shots), 'shots;', int(shots.goal.sum()), 'recorded as goals')

# ## Lab 0's mystery: 183 goals in the shot file, 169 in the results
# Shoot-out kicks (period 5) are recorded as shots; own goals are not shots at all.

in_game = shots[shots.period < 5]
print('in-game shots:', len(in_game), '| in-game goals from shots:', int(in_game.goal.sum()))
print('own goals:', len(own), '| total from results:', int((matches.home_score + matches.away_score).sum()))
print('base rate (in-game):', round(in_game.goal.mean(), 3))

# ## Geometry: distance and the angle the goal gives you
# StatsBomb coordinates: the pitch is 120 × 80 yards, the attacking team shoots towards x = 120, and the posts sit at y = 36 and y = 44.

def geometry(df):
    d = np.hypot(120 - df.x, 40 - df.y)
    ax, ay = 120 - df.x, 36 - df.y
    bx, by = 120 - df.x, 44 - df.y
    cos = (ax * bx + ay * by) / (np.hypot(ax, ay) * np.hypot(bx, by))
    return d, np.arccos(np.clip(cos, -1, 1))
in_game = in_game.assign(dist=geometry(in_game)[0], angle=geometry(in_game)[1])
print('median shot distance:', round(in_game.dist.median(), 1), 'yd; median goal distance:', round(in_game[in_game.goal == 1].dist.median(), 1), 'yd')

# ## Fit the model on open-play shots
# $z = b_0 + b_1\,\text{angle} + b_2\,\text{distance}$, $\;xG = 1/(1+e^{-z})$. Unpenalised maximum likelihood (statsmodels `Logit`).

op = in_game[in_game.shot_type == 'Open Play']
X = sm.add_constant(op[['angle', 'dist']])
fit = sm.Logit(op.goal, X).fit(disp=0)
print(fit.params.round(4))
p = fit.predict(X)
print('open-play shots:', len(op), '| goals:', int(op.goal.sum()), '| sum of fitted xG:', round(p.sum(), 1))
print('Brier, our model:', round(((p - op.goal) ** 2).mean(), 4))
print('Brier, StatsBomb xG:', round(((op.shot_statsbomb_xg - op.goal) ** 2).mean(), 4))

# ## Calibration by distance band
# Fitted mean xG against the actual conversion rate, and the illustrative teaching coefficients used in Session 3 ($b_0=-0.6$, $b_1=2.2$, $b_2=-0.1$) for comparison. The teaching model has the right shape and the wrong level: that is the point of Session 5.

sig = lambda z: 1 / (1 + np.exp(-z))
op = op.assign(fitted=p, teaching=sig(-0.6 + 2.2 * op.angle - 0.1 * op.dist))
bands = pd.cut(op.dist, [0, 6, 12, 18, 24, 30, 120], right=False)
tab = op.groupby(bands, observed=True).agg(shots=('goal', 'size'), goals=('goal', 'sum'), actual=('goal', 'mean'), fitted=('fitted', 'mean'), teaching=('teaching', 'mean')).round(3)
print(tab)
ax = tab[['actual', 'fitted', 'teaching']].plot.bar(figsize=(8, 4), color=['#1f5f63', '#a3471c', '#b0862b'])
ax.set_ylabel('conversion / mean xG'); ax.set_xlabel('distance band (yards)'); ax.set_title('2018 World Cup open-play shots')
plt.tight_layout(); plt.savefig("xg_calibration.png", dpi=120)

# ## Where geometry stops working: the penalty spot

pens = in_game[in_game.shot_type == 'Penalty']
pp = fit.predict(sm.add_constant(pens[['angle', 'dist']], has_constant='add'))
print('in-game penalties:', len(pens), '| scored:', int(pens.goal.sum()), '| actual rate:', round(pens.goal.mean(), 2), '| geometry model predicts:', round(pp.mean(), 2))

# ## Limits
# In-sample scores flatter every model; the course's Session 6 refits on the group stage and tests on the knockout rounds. A geometry-only model does not know about the goalkeeper, defenders, body part or shot type, which is why the penalty gap appears. One tournament is a small sample: 1,556 open-play shots and 129 goals.