Capstone - Complete Analytics System

import pandas as pd
import numpy as np
from statsbombpy import sb
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.metrics import roc_auc_score, brier_score_loss, calibration_curve
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

# Load data from multiple competitions
competitions = sb.competitions()
target_comps = competitions[competitions['competition_name'].isin(['La Liga', 'Premier League'])]

all_shots = []
for _, comp in target_comps.head(2).iterrows():
    matches = sb.matches(competition_id=comp['competition_id'],
                         season_id=comp['season_id'])
    for _, match in matches.head(20).iterrows():
        events = sb.events(match_id=match['match_id'])
        shots = events[events['type'] == 'Shot'].copy()
        if len(shots) > 0:
            all_shots.append(shots)

shots_df = pd.concat(all_shots, ignore_index=True)
print(f"Total shots loaded: {len(shots_df)}")

# Feature engineering
def engineer_features(df):
    features = pd.DataFrame()

    # Extract location coordinates
    features['x'] = df['location'].apply(lambda loc: loc[0] if isinstance(loc, list) else np.nan)
    features['y'] = df['location'].apply(lambda loc: loc[1] if isinstance(loc, list) else np.nan)

    # Geometric features
    features['distance'] = np.sqrt((120 - features['x'])**2 + (40 - features['y'])**2)
    features['angle'] = np.arctan2(np.abs(features['y'] - 40), 120 - features['x']) * 180 / np.pi

    # Body part features
    features['is_header'] = (df['shot_body_part'] == 'Head').astype(int)
    features['is_right_foot'] = (df['shot_body_part'] == 'Right Foot').astype(int)
    features['is_left_foot'] = (df['shot_body_part'] == 'Left Foot').astype(int)

    # Shot type features
    features['is_open_play'] = (df['shot_type'] == 'Open Play').astype(int)
    features['is_free_kick'] = (df['shot_type'] == 'Free Kick').astype(int)
    features['is_corner'] = (df['shot_type'] == 'Corner').astype(int)

    # Game state features
    features['minute_normalized'] = df['minute'] / 90
    features['is_first_half'] = (df['period'] == 1).astype(int)

    # Shot context
    features['under_pressure'] = df['under_pressure'].fillna(False).astype(int)
    features['first_time'] = df['shot_first_time'].fillna(False).astype(int)

    # Target variable
    features['is_goal'] = (df['shot_outcome'] == 'Goal').astype(int)

    # Provider xG for comparison
    features['provider_xg'] = df['shot_statsbomb_xg']

    return features

# Engineer features and clean data
shots_features = engineer_features(shots_df)
shots_features = shots_features[shots_df['shot_type'] != 'Penalty']  # Exclude penalties
shots_clean = shots_features.dropna()

print(f"Shots for modeling: {len(shots_clean)}")
print(f"Goal conversion rate: {shots_clean['is_goal'].mean()*100:.2f}%")

# Define feature columns
feature_cols = ['distance', 'angle', 'is_header', 'is_right_foot', 'is_left_foot',
                'is_open_play', 'is_free_kick', 'is_corner', 'minute_normalized',
                'is_first_half', 'under_pressure', 'first_time']

X = shots_clean[feature_cols]
y = shots_clean['is_goal']
provider_xg = shots_clean['provider_xg']

# Split data
X_train, X_test, y_train, y_test, provider_train, provider_test = train_test_split(
    X, y, provider_xg, test_size=0.2, random_state=42, stratify=y
)

# Scale features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Model 1: Logistic Regression
print("\n=== LOGISTIC REGRESSION ===")
lr_model = LogisticRegression(max_iter=1000, random_state=42)
lr_model.fit(X_train_scaled, y_train)
lr_pred = lr_model.predict_proba(X_test_scaled)[:, 1]

lr_auc = roc_auc_score(y_test, lr_pred)
lr_brier = brier_score_loss(y_test, lr_pred)
print(f"AUC: {lr_auc:.4f}")
print(f"Brier Score: {lr_brier:.4f}")

# Feature importance
lr_importance = pd.DataFrame({
    'feature': feature_cols,
    'coefficient': lr_model.coef_[0]
}).sort_values('coefficient', key=abs, ascending=False)
print("\nTop Features (Logistic Regression):")
print(lr_importance.head())

# Model 2: Random Forest
print("\n=== RANDOM FOREST ===")
rf_model = RandomForestClassifier(n_estimators=100, max_depth=6, random_state=42)
rf_model.fit(X_train, y_train)
rf_pred = rf_model.predict_proba(X_test)[:, 1]

rf_auc = roc_auc_score(y_test, rf_pred)
rf_brier = brier_score_loss(y_test, rf_pred)
print(f"AUC: {rf_auc:.4f}")
print(f"Brier Score: {rf_brier:.4f}")

# Model 3: Gradient Boosting
print("\n=== GRADIENT BOOSTING ===")
gb_model = GradientBoostingClassifier(n_estimators=100, max_depth=4,
                                       learning_rate=0.1, random_state=42)
gb_model.fit(X_train, y_train)
gb_pred = gb_model.predict_proba(X_test)[:, 1]

gb_auc = roc_auc_score(y_test, gb_pred)
gb_brier = brier_score_loss(y_test, gb_pred)
print(f"AUC: {gb_auc:.4f}")
print(f"Brier Score: {gb_brier:.4f}")

# Feature importance for Gradient Boosting
gb_importance = pd.DataFrame({
    'feature': feature_cols,
    'importance': gb_model.feature_importances_
}).sort_values('importance', ascending=False)
print("\nTop Features (Gradient Boosting):")
print(gb_importance.head())

# Model Comparison
print("\n=== MODEL COMPARISON ===")
provider_auc = roc_auc_score(y_test, provider_test)
provider_brier = brier_score_loss(y_test, provider_test)

comparison = pd.DataFrame({
    'Model': ['Logistic Regression', 'Random Forest', 'Gradient Boosting', 'Provider xG'],
    'AUC': [lr_auc, rf_auc, gb_auc, provider_auc],
    'Brier Score': [lr_brier, rf_brier, gb_brier, provider_brier]
})
print(comparison.to_string(index=False))

# Calibration Plots
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

models = [
    ('Logistic Regression', lr_pred, 'blue'),
    ('Gradient Boosting', gb_pred, 'green'),
    ('Provider xG', provider_test.values, 'purple')
]

for ax, (name, pred, color) in zip(axes.flatten()[:3], models):
    prob_true, prob_pred = calibration_curve(y_test, pred, n_bins=10)
    ax.plot(prob_pred, prob_true, 'o-', color=color, label=name)
    ax.plot([0, 1], [0, 1], 'r--', label='Perfect calibration')
    ax.set_xlabel('Predicted xG')
    ax.set_ylabel('Actual Goal Rate')
    ax.set_title(f'{name} Calibration')
    ax.legend()
    ax.set_xlim(0, 0.6)
    ax.set_ylim(0, 0.6)

# All models comparison
ax = axes[1, 1]
for name, pred, color in models:
    prob_true, prob_pred = calibration_curve(y_test, pred, n_bins=10)
    ax.plot(prob_pred, prob_true, 'o-', color=color, label=name)
ax.plot([0, 1], [0, 1], 'r--', label='Perfect')
ax.set_xlabel('Predicted')
ax.set_ylabel('Actual')
ax.set_title('All Models Comparison')
ax.legend()
ax.set_xlim(0, 0.6)
ax.set_ylim(0, 0.6)

plt.tight_layout()
plt.savefig('xg_model_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

print("\nAnalysis complete! Best model determined by lowest Brier score.")

import pandas as pd
import numpy as np
from statsbombpy import sb
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import roc_auc_score, brier_score_loss
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

# Load shot data
competitions = sb.competitions()
la_liga = competitions[competitions['competition_name'] == 'La Liga'].iloc[0]
matches = sb.matches(competition_id=la_liga['competition_id'],
                     season_id=la_liga['season_id'])

all_shots = []
for _, match in matches.head(50).iterrows():
    events = sb.events(match_id=match['match_id'])
    shots = events[events['type'] == 'Shot'].copy()
    if len(shots) > 0:
        all_shots.append(shots)

shots_df = pd.concat(all_shots, ignore_index=True)
shots_df = shots_df[shots_df['shot_type'] != 'Penalty']

# Engineer base features
def engineer_base_features(df):
    features = pd.DataFrame()
    features['x'] = df['location'].apply(lambda loc: loc[0] if isinstance(loc, list) else np.nan)
    features['y'] = df['location'].apply(lambda loc: loc[1] if isinstance(loc, list) else np.nan)
    features['distance'] = np.sqrt((120 - features['x'])**2 + (40 - features['y'])**2)
    features['angle'] = np.arctan2(np.abs(features['y'] - 40), 120 - features['x']) * 180 / np.pi
    features['is_header'] = (df['shot_body_part'] == 'Head')
    features['is_goal'] = (df['shot_outcome'] == 'Goal')
    features['under_pressure'] = df['under_pressure'].fillna(False)
    features['first_time'] = df['shot_first_time'].fillna(False)
    features['is_open_play'] = (df['shot_type'] == 'Open Play')
    features['is_corner'] = (df['shot_type'] == 'Corner')
    features['is_free_kick'] = (df['shot_type'] == 'Free Kick')
    features['minute_norm'] = df['minute'] / 90
    features['provider_xg'] = df['shot_statsbomb_xg']
    return features

shots_features = engineer_base_features(shots_df).dropna()

# Split into headers and foot shots
headers = shots_features[shots_features['is_header'] == True].copy()
foot_shots = shots_features[shots_features['is_header'] == False].copy()

print(f"Total shots: {len(shots_features)}")
print(f"Headers: {len(headers)} ({headers['is_goal'].mean()*100:.2f}% conversion)")
print(f"Foot shots: {len(foot_shots)} ({foot_shots['is_goal'].mean()*100:.2f}% conversion)")

# Engineer header-specific features
headers['from_corner'] = headers['is_corner']
headers['central_position'] = (np.abs(headers['y'] - 40) < 10)
headers['close_range'] = (headers['distance'] < 12)
headers['six_yard_box'] = (headers['x'] > 114) & (np.abs(headers['y'] - 40) < 10)

# Engineer foot shot-specific features
foot_shots['is_right_foot_zone'] = (foot_shots['y'] < 40)
foot_shots['likely_one_on_one'] = (foot_shots['distance'] < 15) & foot_shots['is_open_play']
foot_shots['long_range'] = (foot_shots['distance'] > 20)
foot_shots['tight_angle'] = (foot_shots['angle'] > 30)

# Train header model
print("\n=== HEADER MODEL ===")
header_features = ['distance', 'angle', 'under_pressure', 'first_time',
                   'from_corner', 'central_position', 'close_range', 'six_yard_box']

X_header = headers[header_features].astype(float)
y_header = headers['is_goal'].astype(int)

X_header_train, X_header_test, y_header_train, y_header_test = train_test_split(
    X_header, y_header, test_size=0.2, random_state=42, stratify=y_header
)

header_model = GradientBoostingClassifier(n_estimators=100, max_depth=3,
                                          learning_rate=0.1, random_state=42)
header_model.fit(X_header_train, y_header_train)
header_pred = header_model.predict_proba(X_header_test)[:, 1]

header_auc = roc_auc_score(y_header_test, header_pred)
header_brier = brier_score_loss(y_header_test, header_pred)
print(f"Header Model AUC: {header_auc:.4f}")
print(f"Header Model Brier: {header_brier:.4f}")

# Train foot shot model
print("\n=== FOOT SHOT MODEL ===")
foot_features = ['distance', 'angle', 'under_pressure', 'first_time',
                 'is_open_play', 'is_right_foot_zone', 'likely_one_on_one',
                 'long_range', 'tight_angle']

X_foot = foot_shots[foot_features].astype(float)
y_foot = foot_shots['is_goal'].astype(int)

X_foot_train, X_foot_test, y_foot_train, y_foot_test = train_test_split(
    X_foot, y_foot, test_size=0.2, random_state=42, stratify=y_foot
)

foot_model = GradientBoostingClassifier(n_estimators=100, max_depth=4,
                                        learning_rate=0.1, random_state=42)
foot_model.fit(X_foot_train, y_foot_train)
foot_pred = foot_model.predict_proba(X_foot_test)[:, 1]

foot_auc = roc_auc_score(y_foot_test, foot_pred)
foot_brier = brier_score_loss(y_foot_test, foot_pred)
print(f"Foot Shot Model AUC: {foot_auc:.4f}")
print(f"Foot Shot Model Brier: {foot_brier:.4f}")

# Train unified model for comparison
print("\n=== UNIFIED MODEL ===")
unified_features = ['distance', 'angle', 'under_pressure', 'first_time',
                    'is_open_play', 'is_corner', 'is_header']

X_unified = shots_features[unified_features].astype(float)
y_unified = shots_features['is_goal'].astype(int)

X_unified_train, X_unified_test, y_unified_train, y_unified_test = train_test_split(
    X_unified, y_unified, test_size=0.2, random_state=42, stratify=y_unified
)

# Store test indices for ensemble evaluation
test_indices = X_unified_test.index

unified_model = GradientBoostingClassifier(n_estimators=100, max_depth=4,
                                           learning_rate=0.1, random_state=42)
unified_model.fit(X_unified_train, y_unified_train)
unified_pred = unified_model.predict_proba(X_unified_test)[:, 1]

unified_auc = roc_auc_score(y_unified_test, unified_pred)
unified_brier = brier_score_loss(y_unified_test, unified_pred)
print(f"Unified Model AUC: {unified_auc:.4f}")
print(f"Unified Model Brier: {unified_brier:.4f}")

# Create ensemble prediction
print("\n=== ENSEMBLE MODEL ===")
def ensemble_predict(row, header_model, foot_model, header_features, foot_features):
    if row['is_header']:
        row_features = {
            'distance': row['distance'], 'angle': row['angle'],
            'under_pressure': row['under_pressure'], 'first_time': row['first_time'],
            'from_corner': row['is_corner'],
            'central_position': abs(row['y'] - 40) < 10,
            'close_range': row['distance'] < 12,
            'six_yard_box': row['x'] > 114 and abs(row['y'] - 40) < 10
        }
        X = pd.DataFrame([row_features])[header_features].astype(float)
        return header_model.predict_proba(X)[0, 1]
    else:
        row_features = {
            'distance': row['distance'], 'angle': row['angle'],
            'under_pressure': row['under_pressure'], 'first_time': row['first_time'],
            'is_open_play': row['is_open_play'],
            'is_right_foot_zone': row['y'] < 40,
            'likely_one_on_one': row['distance'] < 15 and row['is_open_play'],
            'long_range': row['distance'] > 20,
            'tight_angle': row['angle'] > 30
        }
        X = pd.DataFrame([row_features])[foot_features].astype(float)
        return foot_model.predict_proba(X)[0, 1]

# Apply ensemble to test set
test_data = shots_features.loc[test_indices].copy()
ensemble_preds = test_data.apply(
    lambda row: ensemble_predict(row, header_model, foot_model,
                                  header_features, foot_features), axis=1
)

ensemble_auc = roc_auc_score(y_unified_test, ensemble_preds)
ensemble_brier = brier_score_loss(y_unified_test, ensemble_preds)
print(f"Ensemble Model AUC: {ensemble_auc:.4f}")
print(f"Ensemble Model Brier: {ensemble_brier:.4f}")

# Final Comparison
print("\n=== FINAL COMPARISON ===")
results = pd.DataFrame({
    'Model': ['Unified Model', 'Ensemble (Header + Foot)',
              'Header-Only Test', 'Foot-Only Test'],
    'AUC': [unified_auc, ensemble_auc, header_auc, foot_auc],
    'Brier Score': [unified_brier, ensemble_brier, header_brier, foot_brier]
})
print(results.to_string(index=False))

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Header conversion by distance
dist_bins = pd.cut(headers['distance'], bins=range(0, 35, 5))
header_conv = headers.groupby(dist_bins)['is_goal'].mean() * 100
axes[0].bar(range(len(header_conv)), header_conv.values, color='coral')
axes[0].set_xticks(range(len(header_conv)))
axes[0].set_xticklabels([f'{i*5}-{(i+1)*5}' for i in range(len(header_conv))])
axes[0].set_xlabel('Distance (m)')
axes[0].set_ylabel('Conversion %')
axes[0].set_title('Header Conversion by Distance')

# Foot shot conversion by distance
dist_bins = pd.cut(foot_shots['distance'], bins=range(0, 45, 5))
foot_conv = foot_shots.groupby(dist_bins)['is_goal'].mean() * 100
axes[1].bar(range(len(foot_conv)), foot_conv.values, color='steelblue')
axes[1].set_xticks(range(len(foot_conv)))
axes[1].set_xticklabels([f'{i*5}-{(i+1)*5}' for i in range(len(foot_conv))])
axes[1].set_xlabel('Distance (m)')
axes[1].set_ylabel('Conversion %')
axes[1].set_title('Foot Shot Conversion by Distance')

plt.tight_layout()
plt.savefig('body_part_xg_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

print("\nConclusion: Ensemble model routes shots to specialized models for better performance.")

import pandas as pd
import numpy as np
from statsbombpy import sb
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import roc_auc_score, brier_score_loss
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

# Load shot data
competitions = sb.competitions()
la_liga = competitions[competitions['competition_name'] == 'La Liga'].iloc[0]
matches = sb.matches(competition_id=la_liga['competition_id'],
                     season_id=la_liga['season_id'])

all_shots = []
for _, match in matches.head(50).iterrows():
    events = sb.events(match_id=match['match_id'])
    shots = events[events['type'] == 'Shot'].copy()
    if len(shots) > 0:
        all_shots.append(shots)

shots_df = pd.concat(all_shots, ignore_index=True)
shots_df = shots_df[shots_df['shot_type'] != 'Penalty']

# Extract shot end locations and engineer features
def extract_psxg_features(df):
    features = pd.DataFrame()

    # Basic location
    features['x'] = df['location'].apply(lambda loc: loc[0] if isinstance(loc, list) else np.nan)
    features['y'] = df['location'].apply(lambda loc: loc[1] if isinstance(loc, list) else np.nan)
    features['distance'] = np.sqrt((120 - features['x'])**2 + (40 - features['y'])**2)
    features['angle'] = np.arctan2(np.abs(features['y'] - 40), 120 - features['x']) * 180 / np.pi

    # Shot end location (goal frame)
    features['end_y'] = df['shot_end_location'].apply(
        lambda loc: loc[1] if isinstance(loc, list) and len(loc) >= 2 else np.nan
    )
    features['end_z'] = df['shot_end_location'].apply(
        lambda loc: loc[2] if isinstance(loc, list) and len(loc) >= 3 else np.nan
    )

    # Normalize to goal frame (goal: 36.6 to 43.4 in y, 0 to 2.44 in z)
    features['goal_y_normalized'] = (features['end_y'] - 36.6) / 7.32
    features['goal_z_normalized'] = features['end_z'] / 2.44

    # Clip to valid range
    features['goal_y_normalized'] = features['goal_y_normalized'].clip(0, 1)
    features['goal_z_normalized'] = features['goal_z_normalized'].clip(0, 1)

    # Shot placement features
    features['dist_from_center'] = np.sqrt(
        (features['goal_y_normalized'] - 0.5)**2 +
        (features['goal_z_normalized'] - 0.5)**2
    )

    features['is_corner_shot'] = (
        ((features['goal_y_normalized'] < 0.15) | (features['goal_y_normalized'] > 0.85)) &
        ((features['goal_z_normalized'] > 0.7) | (features['goal_z_normalized'] < 0.15))
    ).astype(int)

    features['near_post'] = (
        ((features['y'] < 40) & (features['goal_y_normalized'] < 0.3)) |
        ((features['y'] >= 40) & (features['goal_y_normalized'] > 0.7))
    ).astype(int)

    # Other features
    features['is_header'] = (df['shot_body_part'] == 'Head').astype(int)
    features['under_pressure'] = df['under_pressure'].fillna(False).astype(int)

    # Target and comparison
    features['is_goal'] = (df['shot_outcome'] == 'Goal').astype(int)
    features['is_on_target'] = df['shot_outcome'].isin(['Goal', 'Saved', 'Saved to Post']).astype(int)
    features['pre_shot_xg'] = df['shot_statsbomb_xg']

    # Player info
    features['player_name'] = df['player']
    features['team_name'] = df['team']

    return features

# Engineer features
shots_features = extract_psxg_features(shots_df)

# Filter to on-target shots with valid placement data
on_target = shots_features[
    (shots_features['is_on_target'] == 1) &
    shots_features['goal_y_normalized'].notna() &
    shots_features['goal_z_normalized'].notna()
].copy()

print(f"On-target shots with placement data: {len(on_target)}")
print(f"Goal conversion on target: {on_target['is_goal'].mean()*100:.2f}%")

# PSxG model features
psxg_features = ['distance', 'angle', 'goal_y_normalized', 'goal_z_normalized',
                 'dist_from_center', 'is_corner_shot', 'near_post',
                 'is_header', 'under_pressure']

# Prepare clean data
model_data = on_target.dropna(subset=psxg_features + ['is_goal', 'pre_shot_xg'])
print(f"Shots for modeling: {len(model_data)}")

# Split data
X = model_data[psxg_features]
y = model_data['is_goal']

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

# Train PSxG model
psxg_model = GradientBoostingClassifier(
    n_estimators=150, max_depth=4, learning_rate=0.1, random_state=42
)
psxg_model.fit(X_train, y_train)

# Predictions
psxg_pred = psxg_model.predict_proba(X_test)[:, 1]
pre_xg_test = model_data.loc[X_test.index, 'pre_shot_xg']

# Model evaluation
print("\n=== MODEL COMPARISON ===")
psxg_auc = roc_auc_score(y_test, psxg_pred)
psxg_brier = brier_score_loss(y_test, psxg_pred)
pre_xg_auc = roc_auc_score(y_test, pre_xg_test)
pre_xg_brier = brier_score_loss(y_test, pre_xg_test)

print(f"PSxG Model AUC: {psxg_auc:.4f}")
print(f"PSxG Model Brier: {psxg_brier:.4f}")
print(f"Pre-shot xG AUC: {pre_xg_auc:.4f}")
print(f"Pre-shot xG Brier: {pre_xg_brier:.4f}")

# Apply to full dataset
model_data['psxg'] = psxg_model.predict_proba(model_data[psxg_features])[:, 1]

# Finisher analysis
print("\n=== TOP FINISHERS (Goals Above PSxG) ===")
finisher_analysis = model_data.groupby('player_name').agg({
    'is_goal': ['count', 'sum'],
    'psxg': 'sum',
    'pre_shot_xg': 'sum'
}).reset_index()

finisher_analysis.columns = ['player_name', 'shots_on_target', 'goals',
                              'total_psxg', 'total_pre_xg']
finisher_analysis['goals_minus_psxg'] = finisher_analysis['goals'] - finisher_analysis['total_psxg']
finisher_analysis['conversion_rate'] = (finisher_analysis['goals'] /
                                         finisher_analysis['shots_on_target'] * 100)

# Filter players with enough shots
finisher_analysis = finisher_analysis[finisher_analysis['shots_on_target'] >= 10]
finisher_analysis = finisher_analysis.sort_values('goals_minus_psxg', ascending=False)

print("\nBest Finishers (overperforming PSxG):")
print(finisher_analysis.head(10).to_string(index=False))

print("\nWorst Finishers (underperforming PSxG):")
print(finisher_analysis.tail(10).to_string(index=False))

# Goalkeeper analysis (simulated by opponent team)
print("\n=== GOALKEEPER ANALYSIS (by opponent team) ===")
gk_analysis = model_data.groupby('team_name').agg({
    'is_goal': ['count', 'sum'],
    'psxg': 'sum'
}).reset_index()

gk_analysis.columns = ['team', 'shots_faced', 'goals_conceded', 'total_psxg_faced']
gk_analysis['gsax'] = gk_analysis['total_psxg_faced'] - gk_analysis['goals_conceded']
gk_analysis['save_pct'] = ((gk_analysis['shots_faced'] - gk_analysis['goals_conceded']) /
                           gk_analysis['shots_faced'] * 100)

gk_analysis = gk_analysis[gk_analysis['shots_faced'] >= 15]
gk_analysis = gk_analysis.sort_values('gsax', ascending=False)

print("\nBest Goalkeeping (highest GSAx):")
print(gk_analysis.head(5).to_string(index=False))

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Shot placement on goal frame
ax1 = axes[0]
goals = model_data[model_data['is_goal'] == 1]
saves = model_data[model_data['is_goal'] == 0]

ax1.scatter(saves['goal_y_normalized'], saves['goal_z_normalized'],
            c='blue', alpha=0.3, s=30, label='Save')
ax1.scatter(goals['goal_y_normalized'], goals['goal_z_normalized'],
            c='red', alpha=0.6, s=50, label='Goal')

# Draw goal frame
ax1.axhline(y=0, color='black', linewidth=2)
ax1.axhline(y=1, color='black', linewidth=2)
ax1.axvline(x=0, color='black', linewidth=2)
ax1.axvline(x=1, color='black', linewidth=2)

ax1.set_xlim(-0.1, 1.1)
ax1.set_ylim(-0.1, 1.1)
ax1.set_xlabel('Goal Width (Left Post to Right Post)')
ax1.set_ylabel('Goal Height (Ground to Crossbar)')
ax1.set_title('Shot Placement (On Target Shots)')
ax1.legend()

# PSxG vs Pre-shot xG
ax2 = axes[1]
test_goals = y_test == 1
ax2.scatter(pre_xg_test[~test_goals], psxg_pred[~test_goals],
            c='blue', alpha=0.4, s=30, label='No Goal')
ax2.scatter(pre_xg_test[test_goals], psxg_pred[test_goals],
            c='red', alpha=0.6, s=50, label='Goal')
ax2.plot([0, 1], [0, 1], 'k--', alpha=0.5, label='Equal')

ax2.set_xlabel('Pre-shot xG')
ax2.set_ylabel('Post-shot xG (PSxG)')
ax2.set_title('PSxG vs Pre-shot xG')
ax2.legend()
ax2.set_xlim(0, 0.8)
ax2.set_ylim(0, 1)

plt.tight_layout()
plt.savefig('psxg_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

# Shot placement conversion heatmap
fig, ax = plt.subplots(figsize=(8, 6))

# Create bins for heatmap
y_bins = pd.cut(model_data['goal_y_normalized'], bins=5, labels=['Far Left', 'Left', 'Center', 'Right', 'Far Right'])
z_bins = pd.cut(model_data['goal_z_normalized'], bins=3, labels=['Low', 'Mid', 'High'])

placement_matrix = model_data.groupby([z_bins, y_bins])['is_goal'].mean().unstack() * 100

im = ax.imshow(placement_matrix.values, cmap='RdYlGn', aspect='auto', vmin=30, vmax=90)
ax.set_xticks(range(5))
ax.set_xticklabels(['Far Left', 'Left', 'Center', 'Right', 'Far Right'])
ax.set_yticks(range(3))
ax.set_yticklabels(['Low', 'Mid', 'High'])
ax.set_xlabel('Horizontal Position')
ax.set_ylabel('Vertical Position')
ax.set_title('Goal Conversion Rate by Shot Placement (%)')

# Add text annotations
for i in range(3):
    for j in range(5):
        val = placement_matrix.values[i, j]
        if not np.isnan(val):
            ax.text(j, i, f'{val:.1f}%', ha='center', va='center',
                   color='white' if val > 60 else 'black', fontweight='bold')

plt.colorbar(im, label='Conversion %')
plt.tight_layout()
plt.savefig('shot_placement_heatmap.png', dpi=150, bbox_inches='tight')
plt.show()

print("\nPSxG analysis complete! Use for goalkeeper evaluation (GSAx) and finishing quality.")