2015 SEASON OVERVIEW

drivers_standings_2015 = recent_drivers_standings[recent_drivers_standings['year'] == 2015]
plt.figure(figsize=(18, 12))

for fullName in drivers_standings_2015['fullName'].unique():
    driver_data = drivers_standings_2015[drivers_standings_2015['fullName'] == fullName].sort_values('round')
    plt.plot(driver_data['name'], driver_data['points'], 
             marker='o', linewidth=2.5, markersize=6,
             label=f'{fullName}', alpha=0.8)

plt.ylabel('Cumulative Championship Points', fontsize=14)
plt.xlabel('Grand Prix', fontsize=14)
plt.title("Driver Championship Points During 2015 Season", fontsize=16, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
            

constructor_standings_2015 = recent_constructor_standings[recent_constructor_standings['year'] == 2015]
plt.figure(figsize=(15, 12))
race_names = constructor_standings_2015.groupby('round')['name_y'].first().sort_index()

for constructor in constructor_standings_2015['name_x'].unique():
    constructor_data = constructor_standings_2015[constructor_standings_2015['name_x'] == constructor].sort_values('round')
    plt.plot(constructor_data['name_y'], constructor_data['points'], linewidth=2.5, markersize=6,
             label=f'{constructor}', alpha=0.8)

plt.ylabel('Cumulative Championship Points', fontsize=14)
plt.xlabel('Grand Prix', fontsize=14)
plt.title("Constructor Championship Points During 2015 Season", fontsize=16, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
            

def race_by_race_analysis(self):
        """Analyze qualifying impact for each race in 2015"""
        print(f"\n{'='*80}")
        print("RACE-BY-RACE QUALIFYING IMPACT ANALYSIS")
        print("="*80)
    
        race_analysis = []
    
        for race_id in sorted(self.analysis_data['raceId'].unique()):
            race_data = self.analysis_data[self.analysis_data['raceId'] == race_id]
            finished_race_data = race_data.dropna(subset=['position'])
        
            if len(finished_race_data) >= 10:  # Minimum drivers to calculate correlation
                race_info = race_data.iloc[0]
                correlation = finished_race_data['quali_position'].corr(finished_race_data['position'])
                avg_position_change = race_data['position_change'].mean()
                pole_winner = (finished_race_data['quali_position'] == 1) & (finished_race_data['position'] == 1)
                pole_won = pole_winner.any()
            
                race_analysis.append({
                    'round': race_info['round'],
                    'race_name': race_info['race_name'],
                    'circuit': race_info['circuit_name'],
                    'correlation': correlation,
                    'avg_position_change': avg_position_change,
                    'pole_winner': pole_won,
                    'finishers': len(finished_race_data),
                    'predictability': 'High' if correlation > 0.7 else 'Medium' if correlation > 0.5 else 'Low'
                })
    
        race_df = pd.DataFrame(race_analysis)
        race_df = race_df.sort_values('correlation', ascending = False)
    
        print(f"{'Round':5} {'Race':25} {'Predictability':13} {'Correlation':12} {'Avg Pos Change':15} {'Pole Winner'}")
        print("-" * 80)
    
        for _, row in race_df.iterrows():
            pole_symbol = "✓" if row['pole_winner'] else "✗"
            print(f"{row['round']:5} {row['race_name'][:24]:25} {row['predictability']:13} {row['correlation']:12.3f} "
                f"{row['avg_position_change']:15.2f} {pole_symbol}")
    
        print(f"\nRace Analysis Summary:")
        print(f"Average correlation across races: {race_df['correlation'].mean():.3f}")
        print(f"Races where pole position won: {race_df['pole_winner'].sum()}/{len(race_df)}")
        print(f"Most predictable race (highest correlation): {race_df.loc[race_df['correlation'].idxmax(), 'race_name']}")
        print(f"Most unpredictable race (lowest correlation): {race_df.loc[race_df['correlation'].idxmin(), 'race_name']}")
    
        return race_df


analyzer = F1_2015_QualifyingAnalysis()
results = analyzer.run_complete_analysis()
correlation = results['correlation']
driver_stats = results['driver_stats']
full_dataset = results['full_data']
analyzer.race_by_race_analysis()
analyzer.driver_analysis()
            

class F1_2015_QualifyingAnalysis:
    def __init__(self):
        """Initialize 2015 F1 Qualifying Impact Analysis"""
        self.analysis_data = None
        self.prepare_analysis_data()
        
    def prepare_analysis_data(self):
        """Merge qualifying data with race results for 2015"""
        print("Preparing 2015 F1 Season Qualifying Impact Analysis...")
        quali_data = qualifying_2015[['raceId', 'driverId', 'position', 'q1', 'q2', 'q3']].copy()
        quali_data.rename(columns={'position': 'quali_position'}, inplace=True)
        self.analysis_data = performance_data.merge(
            quali_data, 
            on=['raceId', 'driverId'], 
            how='inner'
        )
        
        # Calculate key metrics
        self.analysis_data['position_change'] = (
            self.analysis_data['quali_position'] - self.analysis_data['position']
        )
        self.analysis_data['finished_race'] = ~self.analysis_data['position'].isna()
        self.analysis_data['points_scored'] = self.analysis_data['points'] > 0
        self.analysis_data['top10_finish'] = self.analysis_data['position'] <= 10
        
        # Create qualifying groups
        self.analysis_data['quali_group'] = pd.cut(
            self.analysis_data['quali_position'],
            bins=[0, 3, 10, 20, float('inf')],
            labels=['Top 3', '4th-10th', '11th-20th', 'Back of Grid']
        )
        
        print(f"✓ Dataset prepared: {len(self.analysis_data)} driver-race combinations")
        print(f"✓ Races analyzed: {len(self.analysis_data['round'].unique())} races")
        print(f"✓ Drivers included: {len(self.analysis_data['driverId'].unique())} drivers")
        
    def overall_correlation_analysis(self):
        """Analyze overall qualifying vs race position correlation for 2015"""
        print("\n" + "="*70)
        print("2015 F1 SEASON: QUALIFYING vs RACE POSITION CORRELATION")
        print("="*70)
        
        finished_races = self.analysis_data.dropna(subset=['position'])
        pearson_corr = finished_races['quali_position'].corr(finished_races['position'])
        spearman_corr, spearman_p = stats.spearmanr(
            finished_races['quali_position'], 
            finished_races['position']
        )
        
        print(f"Pearson Correlation: {pearson_corr:.4f}")
        print(f"Spearman Correlation: {spearman_corr:.4f} (p-value: {spearman_p:.2e})")
        print(f"R² (Variance Explained): {pearson_corr**2:.1%}")
        
        if pearson_corr > 0.7:
            strength = "Strong"
        elif pearson_corr > 0.5:
            strength = "Moderate"
        else:
            strength = "Weak"
            
        print(f"Correlation Strength: {strength}")
        
        return pearson_corr, spearman_corr
        
    def race_by_race_analysis(self):
        """Analyze qualifying impact for each race in 2015"""
        print(f"\n{'='*80}")
        print("RACE-BY-RACE QUALIFYING IMPACT ANALYSIS")
        print("="*80)
    
        race_analysis = []
    
        for race_id in sorted(self.analysis_data['raceId'].unique()):
            race_data = self.analysis_data[self.analysis_data['raceId'] == race_id]
            finished_race_data = race_data.dropna(subset=['position'])
        
            if len(finished_race_data) >= 10:  # Minimum drivers to calculate correlation
                race_info = race_data.iloc[0]
                correlation = finished_race_data['quali_position'].corr(finished_race_data['position'])
                avg_position_change = race_data['position_change'].mean()
                pole_winner = (finished_race_data['quali_position'] == 1) & (finished_race_data['position'] == 1)
                pole_won = pole_winner.any()
            
                race_analysis.append({
                    'round': race_info['round'],
                    'race_name': race_info['race_name'],
                    'circuit': race_info['circuit_name'],
                    'correlation': correlation,
                    'avg_position_change': avg_position_change,
                    'pole_winner': pole_won,
                    'finishers': len(finished_race_data),
                    'predictability': 'High' if correlation > 0.7 else 'Medium' if correlation > 0.5 else 'Low'
                })
    
        race_df = pd.DataFrame(race_analysis)
        race_df = race_df.sort_values('correlation', ascending = False)
    
        print(f"{'Round':<5} {'Race':<25} {'Predictability':<13} {'Correlation':<12} {'Avg Pos Change':<15} {'Pole Winner'}")
        print("-" * 80)
    
        for _, row in race_df.iterrows():
            pole_symbol = "✓" if row['pole_winner'] else "✗"
            print(f"{row['round']:<5} {row['race_name'][:24]:<25} {row['predictability']:<13} {row['correlation']:<12.3f} "
                f"{row['avg_position_change']:<15.2f} {pole_symbol}")
    
    # Summary statistics
        print(f"\nRace Analysis Summary:")
        print(f"Average correlation across races: {race_df['correlation'].mean():.3f}")
        print(f"Races where pole position won: {race_df['pole_winner'].sum()}/{len(race_df)}")
        print(f"Most predictable race (highest correlation): {race_df.loc[race_df['correlation'].idxmax(), 'race_name']}")
        print(f"Most unpredictable race (lowest correlation): {race_df.loc[race_df['correlation'].idxmin(), 'race_name']}")
    
        return race_df

        
    def position_change_analysis(self):
        """Analyze how positions change from qualifying to race"""
        print(f"\n{'='*70}")
        print("POSITION CHANGE ANALYSIS (QUALIFYING → RACE)")
        print("="*70)
        
        pos_changes = self.analysis_data['position_change'].dropna()
        
        print(f"Total driver-race combinations: {len(pos_changes)}")
        print(f"Average position change: {pos_changes.mean():.2f}")
        print(f"Median position change: {pos_changes.median():.2f}")
        print(f"Standard deviation: {pos_changes.std():.2f}")
        
        print(f"\nPosition Change Distribution:")
        gained = (pos_changes > 0).sum()
        lost = (pos_changes < 0).sum()
        stayed = (pos_changes == 0).sum()
        
        print(f"Gained positions: {gained} ({gained/len(pos_changes)*100:.1f}%)")
        print(f"Lost positions: {lost} ({lost/len(pos_changes)*100:.1f}%)")
        print(f"No change: {stayed} ({stayed/len(pos_changes)*100:.1f}%)")
        
        print(f"\nExtreme Cases:")
        print(f"Biggest gain: +{pos_changes.max():.0f} positions")
        print(f"Biggest loss: {pos_changes.min():.0f} positions")
        if pos_changes.max() > 10:
            big_gain = self.analysis_data[self.analysis_data['position_change'] == pos_changes.max()].iloc[0]
            print(f"Biggest gain by: {big_gain['fullName']} ({big_gain['race_name']})")
            
        if pos_changes.min() < -10:
            big_loss = self.analysis_data[self.analysis_data['position_change'] == pos_changes.min()].iloc[0]
            print(f"Biggest loss by: {big_loss['fullName']} ({big_loss['race_name']})")
            
        return pos_changes
        
    def qualifying_group_analysis(self):
        """Analyze performance by qualifying position groups"""
        print(f"\n{'='*70}")
        print("PERFORMANCE BY QUALIFYING POSITION GROUPS")
        print("="*70)
        
        group_stats = self.analysis_data.groupby('quali_group').agg({
            'position': ['count', 'mean', 'median'],
            'position_change': ['mean', 'std'],
            'points': ['mean', 'sum'],
            'points_scored': 'mean',
            'top10_finish': 'mean'
        }).round(2)
        
        group_stats.columns = ['races', 'avg_finish', 'median_finish', 'avg_change', 'change_std', 
                              'avg_points', 'total_points', 'points_rate', 'top10_rate']
        
        print(f"{'Group':<15} {'Races':<8} {'Avg Finish':<12} {'Avg Change':<12} {'Points Rate':<12} {'Top10 Rate'}")
        print("-" * 75)
        
        for group, row in group_stats.iterrows():
            print(f"{group:<15} {row['races']:<8.0f} {row['avg_finish']:<12.2f} {row['avg_change']:<12.2f} "
                  f"{row['points_rate']:<12.1%} {row['top10_rate']:<12.1%}")
        
        return group_stats
        
    def driver_analysis(self):
        """Analyze individual driver performance vs qualifying"""
        print(f"\n{'='*70}")
        print("DRIVER QUALIFYING vs RACE PERFORMANCE (2015)")
        print("="*70)
        
        driver_stats = self.analysis_data.groupby(['driverId', 'fullName', 'constructor_name']).agg({
            'quali_position': 'mean',
            'position': 'mean',
            'position_change': ['mean', 'std'],
            'points': 'sum',
            'raceId': 'count'
        }).round(2)
        
        driver_stats.columns = ['avg_quali', 'avg_finish', 'avg_change', 'change_consistency', 'total_points', 'races']
        driver_stats['quali_vs_finish_diff'] = driver_stats['avg_finish'] - driver_stats['avg_quali']
        
        # Sort by total points (championship order)
        driver_stats = driver_stats.sort_values('total_points', ascending=False)
        
        print(f"{'Driver':<20} {'Team':<15} {'Avg Quali':<10} {'Avg Finish':<10} {'Avg Change':<10} {'Points'}")
        print("-" * 85)
        
        for (driver_id, name, team), row in driver_stats.head(15).iterrows():
            print(f"{name[:19]:<20} {team[:14]:<15} {row['avg_quali']:<10.1f} {row['avg_finish']:<10.1f} "
                  f"{row['avg_change']:<10.2f} {row['total_points']:<6.0f}")
        
        return driver_stats
        
    def circuit_analysis(self):
        """Analyze qualifying impact by circuit"""
        print(f"\n{'='*70}")
        print("CIRCUIT-SPECIFIC QUALIFYING IMPACT")
        print("="*70)
        
        circuit_stats = []
        
        for circuit_id in self.analysis_data['circuitId'].unique():
            circuit_data = self.analysis_data[self.analysis_data['circuitId'] == circuit_id]
            finished_data = circuit_data.dropna(subset=['position'])
            
            if len(finished_data) >= 10:
                circuit_info = circuit_data.iloc[0]
                correlation = finished_data['quali_position'].corr(finished_data['position'])
                avg_change = circuit_data['position_change'].mean()
                change_std = circuit_data['position_change'].std()
                
                circuit_stats.append({
                    'circuit_name': circuit_info['circuit_name'],
                    'correlation': correlation,
                    'avg_position_change': avg_change,
                    'position_change_std': change_std,
                    'predictability': 'High' if correlation > 0.7 else 'Medium' if correlation > 0.5 else 'Low'
                })
        
        circuit_df = pd.DataFrame(circuit_stats).sort_values('correlation', ascending=False)
        
        print(f"{'Circuit':<25} {'Correlation':<12} {'Avg Change':<12} {'Predictability'}")
        print("-" * 65)
        
        for _, row in circuit_df.iterrows():
            print(f"{row['circuit_name'][:24]:<25} {row['correlation']:<12.3f} "
                  f"{row['avg_position_change']:<12.2f} {row['predictability']}")
        
        return circuit_df
        
    def create_visualizations(self):
        """Create comprehensive visualizations for 2015 analysis"""
        fig, axes = plt.subplots(2, 3, figsize=(20, 12))
        fig.suptitle('2015 F1 Season: Qualifying Impact Analysis', fontsize=16, fontweight='bold')
        
        # 1. Qualifying vs Race Position Scatter
        finished_data = self.analysis_data.dropna(subset=['position'])
        axes[0, 0].scatter(finished_data['quali_position'], finished_data['position'], 
                          alpha=0.6, s=30, color='red')
        axes[0, 0].plot([1, 22], [1, 22], 'k--', alpha=0.8, linewidth=2, label='Perfect correlation')
        axes[0, 0].set_xlabel('Qualifying Position')
        axes[0, 0].set_ylabel('Race Finish Position')
        axes[0, 0].set_title('Qualifying vs Race Position')
        axes[0, 0].legend()
        axes[0, 0].grid(True, alpha=0.3)
        axes[0, 0].set_xlim(0, 23)
        axes[0, 0].set_ylim(0, 23)
        
        # 2. Position Change Distribution
        pos_changes = self.analysis_data['position_change'].dropna()
        axes[0, 1].hist(pos_changes, bins=30, alpha=0.7, color='blue', edgecolor='black')
        axes[0, 1].axvline(0, color='red', linestyle='--', linewidth=2, label='No change')
        axes[0, 1].axvline(pos_changes.mean(), color='green', linestyle='-', linewidth=2, 
                          label=f'Mean: {pos_changes.mean():.1f}')
        axes[0, 1].set_xlabel('Position Change (Quali → Race)')
        axes[0, 1].set_ylabel('Frequency')
        axes[0, 1].set_title('Distribution of Position Changes')
        axes[0, 1].legend()
        axes[0, 1].grid(True, alpha=0.3)
        
        # 3. Points by Qualifying Position
        quali_points = self.analysis_data.groupby('quali_position')['points'].mean()
        axes[0, 2].bar(quali_points.index, quali_points.values, color='gold', alpha=0.8, edgecolor='black')
        axes[0, 2].set_xlabel('Qualifying Position')
        axes[0, 2].set_xticks(quali_points.index)
        axes[0, 2].set_ylabel('Average Points per Race')
        axes[0, 2].set_title('Average Points by Qualifying Position')
        axes[0, 2].grid(True, alpha=0.3, axis='y')
        
        # 4. Performance by Qualifying Groups
        group_data = []
        group_labels = []
        for group in ['Top 3', '4th-10th', '11th-20th', 'Back of Grid']:
            group_positions = self.analysis_data[self.analysis_data['quali_group'] == group]['position'].dropna()
            if len(group_positions) > 0:
                group_data.append(group_positions)
                group_labels.append(group)
        
        axes[1, 0].boxplot(group_data, labels=group_labels)
        axes[1, 0].set_ylabel('Race Finish Position')
        axes[1, 0].set_title('Race Results by Qualifying Groups')
        axes[1, 0].grid(True, alpha=0.3, axis='y')
        
        # 5. Constructor Performance
        constructor_perf = self.analysis_data.groupby('constructor_name').agg({
            'quali_position': 'mean',
            'position': 'mean',
            'points': 'sum'
        }).sort_values('points', ascending=False).head(10)
        
        x_pos = np.arange(len(constructor_perf))
        axes[1, 1].scatter(constructor_perf['quali_position'], constructor_perf['position'], 
                          s=constructor_perf['points']*2, alpha=0.7, c='red')
        
        for i, (idx, row) in enumerate(constructor_perf.iterrows()):
            axes[1, 1].annotate(idx[:8], (row['quali_position'], row['position']), 
                               xytext=(5, 5), textcoords='offset points', fontsize=8)
        
        axes[1, 1].plot([1, 20], [1, 20], 'k--', alpha=0.5)
        axes[1, 1].set_xlabel('Average Qualifying Position')
        axes[1, 1].set_ylabel('Average Race Position')
        axes[1, 1].set_title('Constructor Performance (Size = Total Points)')
        axes[1, 1].grid(True, alpha=0.3)
        
        # 6. Race-by-Race Correlation
        race_correlations = []
        race_names = []
        
        for race_id in sorted(self.analysis_data['raceId'].unique()):
            race_data = self.analysis_data[self.analysis_data['raceId'] == race_id]
            finished_data = race_data.dropna(subset=['position'])
            
            if len(finished_data) >= 10:
                correlation = finished_data['quali_position'].corr(finished_data['position'])
                race_correlations.append(correlation)
                race_names.append(race_data.iloc[0]['race_name'][:10])
        
        axes[1, 2].bar(range(len(race_correlations)), race_correlations, color='purple', alpha=0.7)
        axes[1, 2].set_xlabel('Race')
        axes[1, 2].set_ylabel('Correlation')
        axes[1, 2].set_title('Qualifying-Race Correlation by Race')
        axes[1, 2].set_xticks(range(len(race_names)))
        axes[1, 2].set_xticklabels(race_names, rotation=45, ha='right')
        axes[1, 2].grid(True, alpha=0.3, axis='y')
        axes[1, 2].axhline(y=0.7, color='red', linestyle='--', alpha=0.7, label='Strong correlation')
        axes[1, 2].legend()
        
        plt.tight_layout()
        plt.show()
            

def plot_race_positions(race_data, race_name):
    """
    Plot position changes for a specific race
    """
    plt.figure(figsize=(12,8))
    
    for driver in race_data['fullName'].unique():
        driver_data = race_data[race_data['fullName'] == driver]
        plt.plot(driver_data['lap'], driver_data['position'], 
                linewidth=2, label=driver, alpha=0.8)
    
    plt.xlabel('Lap Number', fontsize=12)
    plt.ylabel('Position', fontsize=12)
    plt.title(f'Position Changes Throughout {race_name}', fontsize=14)
    plt.gca().invert_yaxis()
    plt.grid(True, alpha=0.3)
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.tight_layout()
    plt.show()

# Dictionary of all race datasets
race_datasets = {
    'Bahrain Grand Prix': new_drivers_2015_BHR,
    'Saudi Arabian Grand Prix': new_drivers_2015_SAU,
    'Australian Grand Prix': new_drivers_2015_AUS,
    'Emilia Romagna Grand Prix': new_drivers_2015_EMI,
    'Miami Grand Prix': new_drivers_2015_MIA,
    'Spanish Grand Prix': new_drivers_2015_ESP,
    'Monaco Grand Prix': new_drivers_2015_MCO,
    'Azerbaijan Grand Prix': new_drivers_2015_AZE,
    'Canadian Grand Prix': new_drivers_2015_CAN,
    'British Grand Prix': new_drivers_2015_GBR,
    'Austrian Grand Prix': new_drivers_2015_AUT,
    'French Grand Prix': new_drivers_2015_FRA,
    'Hungarian Grand Prix': new_drivers_2015_HUN,
    'Belgian Grand Prix': new_drivers_2015_BEL,
    'Dutch Grand Prix': new_drivers_2015_DUT,
    'Italian Grand Prix': new_drivers_2015_ITA,
    'Singapore Grand Prix': new_drivers_2015_SGP,
    'Japanese Grand Prix': new_drivers_2015_JPN,
    'United States Grand Prix': new_drivers_2015_USA,
    'Mexico City Grand Prix': new_drivers_2015_MEX,
    'São Paulo Grand Prix': new_drivers_2015_BRA,
    'Abu Dhabi Grand Prix': new_drivers_2015_ARE
}

# Plot all races
for race_name, race_data in race_datasets.items():
    if len(race_data) > 0:  # Only plot if race has data
        plot_race_positions(race_data, race_name)
    else:
        print(f"No data available for {race_name}")
            

heatmap_data = position_changes_named.pivot_table(
    values='positions_gained', 
    index='fullName', 
    columns='name', 
    fill_value=0
)

plt.figure(figsize=(15, 10))
sns.heatmap(heatmap_data, annot=True, cmap='RdYlGn', center=0, 
            fmt='.0f', cbar_kws={'label': 'Positions Gained/Lost'})
plt.title('Position Changes by Driver and Grand Prix')
plt.xlabel('Grand Prix')
plt.ylabel('Driver')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()
            

def driverstyle_dataframe(df):
    df_display = df.copy()
    df_display = df_display.rename(columns={
        'fullName': 'Name',
        'total_positions_gained': 'Total Gained',
        'avg_positions_per_race': 'Avg per Race',
        'consistency': 'Consistency',
        'best_single_race': 'Best Single Race',
        'worst_single_race': 'Worst Single Race',
        'races_completed': 'Races Completed'
    })
    
    df_sorted = df_display.sort_values('Total Gained', ascending=False)
    
    styled = df_sorted.style.format({
        'Total Gained': '{:+.0f}',
        'Avg per Race': '{:.2f}',
        'Consistency': '{:.2f}',
        'Best Single Race': '{:+.0f}',
        'Worst Single Race': '{:+.0f}',
        'Races Completed': '{:.0f}'
    }).set_caption(
        "Driver Position Change Summary (Sorted by Average Change)"
    )
    
    return styled
driver_styled_table = driverstyle_dataframe(driver_summary)
driver_styled_table
            

races_2015 = [
    'Australian Grand Prix',
    'Malaysian Grand Prix', 
    'Chinese Grand Prix',
    'Bahrain Grand Prix',
    'Spanish Grand Prix',
    'Monaco Grand Prix',
    'Canadian Grand Prix',
    'Austrian Grand Prix',
    'British Grand Prix',
    'Hungarian Grand Prix',
    'Belgian Grand Prix',
    'Italian Grand Prix',
    'Singapore Grand Prix',
    'Japanese Grand Prix',
    'Russian Grand Prix',
    'United States Grand Prix',
    'Mexican Grand Prix',
    'Brazilian Grand Prix',
    'Abu Dhabi Grand Prix'
]

for race_name in races_2015:
    race_changes = position_changes_2015[position_changes_2015['name'] == race_name]
    plt.figure(figsize=(14, 8))
    plt.plot(race_changes['lap'], race_changes['total_position_changes'], 
            marker='o', linewidth=2, markersize=6, color='#E10600')  # F1 red color
    plt.title(f'Total Position Changes Per Lap - 2015 {race_name}', 
                fontsize=16, fontweight='bold')
    plt.xlabel('Lap Number', fontsize=12)
    plt.ylabel('Total Position Changes', fontsize=12)
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()
            

def GP_Position_Change(df):
    df_display = df.copy()
    df_display = df_display.rename(columns={
        'name': 'Grand Prix',
        'total_positions_gained': 'Total Gained',
        'avg_positions_per_driver': 'Avg per Driver',
        'consistency': 'Consistency',
        'biggest_gain': 'Best Gain',
        'biggest_loss': 'Worst Loss',
        'drivers_count': 'Drivers'
    })
    
    df_sorted = df_display.sort_values('Avg per Driver', ascending=False)
    
    styled = df_sorted.style.format({
        'Total Gained': '{:+.0f}',
        'Avg per Driver': '{:+.2f}',
        'Consistency': '{:.2f}',
        'Best Gain': '{:+.0f}',
        'Worst Loss': '{:+.0f}',
        'Drivers': '{:.0f}'
    }).set_caption(
        "Grand Prix Position Change Summary (Sorted by Average Change)")
    
    return styled

styled_table = GP_Position_Change(gp_summary)
styled_table
            

plt.figure(figsize=(12,8))
sns.boxplot(data = hamilton_2015, x ='name', y = 'milliseconds')
plt.xticks(rotation = 45)
plt.xlabel('Grand Prix')
plt.ylabel('Lap Time (milliseconds)')
plt.title('Lap Time for Hamilton in 2015 Season')
plt.show()

            

def calculate_z_scores(data, metric_cols):
    z_scores = {}
    for col in metric_cols:
        if col in data.columns:
            mean_val = data[col].mean()
            std_val = data[col].std()
            z_scores[f'{col}_zscore'] = (data[col] - mean_val) / std_val
            z_scores[f'{col}_mean'] = mean_val
            z_scores[f'{col}_std'] = std_val
    
    return pd.DataFrame(z_scores, index=data.index)
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('F1 2015 Z-Score Distributions - Hamilton vs Field', fontsize=16, fontweight='bold')

# Colors for Hamilton
hamilton_color = '#00D2BE'  # Mercedes teal
field_color = '#2E86AB'     # Blue
highlight_color = '#F18F01' # Orange

# 1. Overall Z-Score Distribution
ax1 = axes[0, 0]
n, bins, patches = ax1.hist(season_stats_2015['overall_zscore'], bins=15, alpha=0.7, 
                           color=field_color, edgecolor='black', label='All Drivers')

# Highlight Hamilton's position
if len(hamilton_data) > 0:
    hamilton_overall = hamilton_row['overall_zscore']
    ax1.axvline(hamilton_overall, color=hamilton_color, linewidth=3, 
                label=f'Hamilton ({hamilton_overall:.3f})')
    
    # Add percentile text
    percentile = stats.percentileofscore(season_stats_2015['overall_zscore'].dropna(), hamilton_overall)
    ax1.text(hamilton_overall + 0.1, ax1.get_ylim()[1] * 0.8, 
             f'{percentile:.1f}th', fontsize=10, fontweight='bold')

ax1.set_xlabel('Overall Z-Score')
ax1.set_ylabel('Number of Drivers')
ax1.set_title('Overall Performance Z-Score Distribution')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Points Z-Score Distribution
ax2 = axes[0, 1]
ax2.hist(season_stats_2015['points_zscore'], bins=15, alpha=0.7, 
         color=field_color, edgecolor='black', label='All Drivers')

if len(hamilton_data) > 0:
    hamilton_points_z = hamilton_row['points_zscore']
    ax2.axvline(hamilton_points_z, color=hamilton_color, linewidth=3, 
                label=f'Hamilton ({hamilton_points_z:.3f})')
    
    percentile = stats.percentileofscore(season_stats_2015['points_zscore'], hamilton_points_z)
    ax2.text(hamilton_points_z + 0.1, ax2.get_ylim()[1] * 0.8, 
             f'{percentile:.1f}th', fontsize=10, fontweight='bold')

ax2.set_xlabel('Points Z-Score')
ax2.set_ylabel('Number of Drivers')
ax2.set_title('Championship Points Z-Score Distribution')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 3. Position Z-Score Distribution
ax3 = axes[0, 2]
ax3.hist(season_stats_2015['position_zscore'], bins=15, alpha=0.7, 
         color=field_color, edgecolor='black', label='All Drivers')

if len(hamilton_data) > 0:
    hamilton_pos_z = hamilton_row['position_zscore']
    ax3.axvline(hamilton_pos_z, color=hamilton_color, linewidth=3, 
                label=f'Hamilton ({hamilton_pos_z:.3f})')
    
    percentile = stats.percentileofscore(season_stats_2015['position_zscore'].dropna(), hamilton_pos_z)
    ax3.text(hamilton_pos_z + 0.1, ax3.get_ylim()[1] * 0.8, 
             f'{percentile:.1f}th', fontsize=10, fontweight='bold')

ax3.set_xlabel('Position Z-Score')
ax3.set_ylabel('Number of Drivers')
ax3.set_title('Average Finishing Position Z-Score Distribution')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Speed Z-Score Distribution
ax4 = axes[1, 0]
speed_data = season_stats_2015['speed_zscore'].dropna()
ax4.hist(speed_data, bins=15, alpha=0.7, 
         color=field_color, edgecolor='black', label='All Drivers')

if len(hamilton_data) > 0 and not pd.isna(hamilton_row['speed_zscore']):
    hamilton_speed_z = hamilton_row['speed_zscore']
    ax4.axvline(hamilton_speed_z, color=hamilton_color, linewidth=3, 
                label=f'Hamilton ({hamilton_speed_z:.3f})')
    
    percentile = stats.percentileofscore(speed_data, hamilton_speed_z)
    ax4.text(hamilton_speed_z + 0.1, ax4.get_ylim()[1] * 0.8, 
             f'{percentile:.1f}th', fontsize=10, fontweight='bold')

ax4.set_xlabel('Speed Z-Score')
ax4.set_ylabel('Number of Drivers')
ax4.set_title('Fastest Lap Speed Z-Score Distribution')
ax4.legend()
ax4.grid(True, alpha=0.3)

# 5. Grid Position Z-Score Distribution
ax5 = axes[1, 1]
ax5.hist(season_stats_2015['grid_zscore'], bins=15, alpha=0.7, 
         color=field_color, edgecolor='black', label='All Drivers')

if len(hamilton_data) > 0:
    hamilton_grid_z = hamilton_row['grid_zscore']
    ax5.axvline(hamilton_grid_z, color=hamilton_color, linewidth=3, 
                label=f'Hamilton ({hamilton_grid_z:.3f})')
    
    percentile = stats.percentileofscore(season_stats_2015['grid_zscore'], hamilton_grid_z)
    ax5.text(hamilton_grid_z + 0.1, ax5.get_ylim()[1] * 0.8, 
             f'{percentile:.1f}th', fontsize=10, fontweight='bold')

ax5.set_xlabel('Grid Position Z-Score')
ax5.set_ylabel('Number of Drivers')
ax5.set_title('Average Grid Position Z-Score Distribution')
ax5.legend()
ax5.grid(True, alpha=0.3)

# 6. Hamilton's Z-Score Profile (Radar Chart Style)
ax6 = axes[1, 2]
if len(hamilton_data) > 0:
    categories = ['Points', 'Position', 'Speed', 'Grid', 'Overall']
    hamilton_scores = [
        hamilton_row['points_zscore'],
        hamilton_row['position_zscore'], 
        hamilton_row['speed_zscore'] if not pd.isna(hamilton_row['speed_zscore']) else 0,
        hamilton_row['grid_zscore'],
        hamilton_row['overall_zscore']
    ]
    
    bars = ax6.bar(categories, hamilton_scores, color=hamilton_color, alpha=0.8, edgecolor='black')
    ax6.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax6.axhline(y=1, color='red', linestyle='--', alpha=0.5, label='1 Std Above Mean')
    ax6.axhline(y=-1, color='red', linestyle='--', alpha=0.5, label='1 Std Below Mean')
    
    for bar, score in zip(bars, hamilton_scores):
        height = bar.get_height()
        ax6.text(bar.get_x() + bar.get_width()/2., height + 0.05 if height >= 0 else height - 0.15,
                f'{score:.2f}', ha='center', va='bottom' if height >= 0 else 'top', fontweight='bold')
    
    ax6.set_ylabel('Z-Score')
    ax6.set_title("Hamilton's Z-Score Profile")
    ax6.legend()
    ax6.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

def season_momentum_comparison(all_momentum):
    
    early_season = all_momentum[all_momentum['round'] <= 7]  # First 7 races
    late_season = all_momentum[all_momentum['round'] >= 13]  # Last 7 races
    early_avg = early_season.groupby('fullName')['position'].mean()
    late_avg = late_season.groupby('fullName')['position'].mean()
    
    momentum_comparison = pd.DataFrame({
        'early_season_avg': early_avg,
        'late_season_avg': late_avg
    }).dropna()
    momentum_comparison['improvement'] = momentum_comparison['early_season_avg'] - momentum_comparison['late_season_avg']
    top_drivers = all_momentum.groupby('fullName')['points'].sum().index
    momentum_subset = momentum_comparison[momentum_comparison.index.isin(top_drivers)]
    
    plt.figure(figsize=(12, 8))
    
    colors = ['green' if x > 0 else 'red' for x in momentum_subset['improvement']]
    plt.scatter(momentum_subset['early_season_avg'], momentum_subset['late_season_avg'], 
               c=colors, s=200, alpha=0.7, edgecolors='black', linewidth=2)
    min_pos = min(momentum_subset['early_season_avg'].min(), momentum_subset['late_season_avg'].min())
    max_pos = max(momentum_subset['early_season_avg'].max(), momentum_subset['late_season_avg'].max())
    plt.plot([min_pos, max_pos], [min_pos, max_pos], 'k--', alpha=0.5, linewidth=2)
    
    for fullName, row in momentum_subset.iterrows():
        plt.annotate(f'{fullName}', 
                    (row['early_season_avg'], row['late_season_avg']),
                    xytext=(5, 5), textcoords='offset points', fontsize=10)
    
    plt.xlabel('Early Season Average Position (Rounds 1-7)', fontsize=12)
    plt.ylabel('Late Season Average Position (Rounds 13-19)', fontsize=12)
    plt.title('Season Momentum: Early vs Late Season Performance', fontsize=16, fontweight='bold')
    plt.grid(True, alpha=0.3)
    plt.scatter([], [], c='green', s=100, label='Improved (green)', alpha=0.7)
    plt.scatter([], [], c='red', s=100, label='Declined (red)', alpha=0.7)
    plt.legend()
    plt.tight_layout()
    plt.show()
    
    return momentum_comparison

            

modeling_data = results_2015.merge(
    drivers[['driverId', 'forename', 'surname']], on='driverId', how='left').merge(
    constructors[['constructorId', 'name']], on='constructorId', how='left', suffixes=('', '_constructor')).merge(
    races_2015[['raceId', 'name', 'round']], on='raceId', how='left', suffixes=('', '_race'))

modeling_data['driver_name'] = modeling_data['forename'] + ' ' + modeling_data['surname']
modeling_data['race_name'] = modeling_data['name_race']
modeling_data['position'] = pd.to_numeric(modeling_data['position'], errors='coerce')
modeling_data['grid'] = pd.to_numeric(modeling_data['grid'], errors='coerce')
modeling_data['points'] = pd.to_numeric(modeling_data['points'], errors='coerce')
modeling_data = modeling_data.dropna(subset=['position', 'grid', 'constructor_name'])
modeling_data = modeling_data[modeling_data['position'] <= 20]
modeling_data = modeling_data[modeling_data['grid'] <= 24]

constructor_model_data = modeling_data.copy()
constructor_dummies = pd.get_dummies(
    constructor_model_data['constructor_name'], 
    prefix='constructor',
    drop_first=False)

if 'constructor_Mercedes' in constructor_dummies.columns:
    constructor_dummies = constructor_dummies.drop('constructor_Mercedes', axis=1)
    reference_constructor = 'Mercedes'
else:
    reference_constructor = sorted(constructor_model_data['constructor_name'].unique())[0]
    constructor_dummies = constructor_dummies.drop(f'constructor_{reference_constructor}', axis=1)

X_constructor = pd.concat([
    constructor_model_data[['grid']],
    constructor_dummies], axis=1)

y = constructor_model_data['position']

constructor_model = LinearRegression()
constructor_model.fit(X_constructor, y)
constructor_r2 = r2_score(y, constructor_model.predict(X_constructor))
grid_coef = constructor_model.coef_[0]
constructor_coefs = dict(zip(constructor_dummies.columns, constructor_model.coef_[1:]))

constructor_effects = {reference_constructor: 0.0}
for col, coef in constructor_coefs.items():
    constructor_name = col.replace('constructor_', '')
    constructor_effects[constructor_name] = coef

modeling_data['constructor_baseline'] = modeling_data['constructor_name'].map(constructor_effects)
modeling_data['expected_position_from_car'] = (
    constructor_model.intercept_ + 
    grid_coef * modeling_data['grid'] + 
    modeling_data['constructor_baseline'])

modeling_data['driver_performance'] = (
    modeling_data['position'] - modeling_data['expected_position_from_car'])

driver_skill = modeling_data.groupby('driver_name').agg({
    'driver_performance': ['mean', 'std', 'count'],
    'position': 'mean',
    'grid': 'mean',
    'points': 'sum'}).round(3)

driver_skill.columns = ['avg_performance_vs_car', 'performance_std', 'race_count', 
                       'avg_position', 'avg_grid', 'total_points']
driver_skill = driver_skill.reset_index()

min_races = 5
driver_skill_filtered = driver_skill[driver_skill['race_count'] >= min_races].copy()
driver_skill_filtered['skill_rank'] = driver_skill_filtered['avg_performance_vs_car'].rank()
driver_skill_filtered = driver_skill_filtered.sort_values('avg_performance_vs_car')

modeling_data['hierarchical_prediction'] = (
    modeling_data['expected_position_from_car'] + 
    modeling_data['driver_name'].map(
        driver_skill.set_index('driver_name')['avg_performance_vs_car']).fillna(0))

hierarchical_r2 = r2_score(modeling_data['position'], modeling_data['hierarchical_prediction'])
hierarchical_mae = np.mean(np.abs(modeling_data['position'] - modeling_data['hierarchical_prediction']))

baseline_model = LinearRegression()
baseline_model.fit(modeling_data[['grid']], modeling_data['position'])
baseline_r2 = r2_score(modeling_data['position'], baseline_model.predict(modeling_data[['grid']]))

total_variance = np.var(modeling_data['position'])
grid_effect_variance = np.var(baseline_model.predict(modeling_data[['grid']]))
constructor_effect_variance = np.var(modeling_data['constructor_baseline'])
driver_effect_variance = np.var(modeling_data['driver_name'].map(
    driver_skill.set_index('driver_name')['avg_performance_vs_car']
).fillna(0))

grid_pct = (grid_effect_variance / total_variance) * 100
constructor_pct = (constructor_effect_variance / total_variance) * 100  
driver_pct = (driver_effect_variance / total_variance) * 100
explained_pct = hierarchical_r2 * 100

if constructor_pct + driver_pct > 0:
    car_vs_driver_ratio = constructor_pct / driver_pct

constructor_ranking = []
for constructor, effect in constructor_effects.items():
    total_points = modeling_data[modeling_data['constructor_name'] == constructor]['points'].sum()
    constructor_ranking.append({
        'constructor': constructor,
        'car_effect': effect,
        'total_points': total_points
    })

constructor_ranking = sorted(constructor_ranking, key=lambda x: x['car_effect'])
expected_top_drivers = ['Lewis Hamilton', 'Sebastian Vettel', 'Nico Rosberg']
expected_top_constructors = ['Mercedes', 'Ferrari', 'Williams']
actual_top_drivers = driver_skill_filtered.head(3)['driver_name'].tolist()
actual_top_constructors = [item['constructor'] for item in constructor_ranking[:3]]
driver_matches = len(set(expected_top_drivers) & set(actual_top_drivers))
constructor_matches = len(set(expected_top_constructors) & set(actual_top_constructors))

fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('F1 2015: Driver Skill vs Car Performance Regression Analysis', 
             fontsize=18, fontweight='bold', y=0.96)

# 1. Model Performance Comparison
ax1 = axes[0, 0]
models = ['Baseline', 'Constructor', 'Driver', 'Full Model']
r2_scores = [r2_baseline, r2_constructor, r2_driver, r2_full]
colors = ['lightgray', 'orange', 'lightblue', 'green']

bars = ax1.bar(models, r2_scores, color=colors, edgecolor='black', alpha=0.8)
ax1.set_ylabel('R² Score', fontsize=12)
ax1.set_title('Model Performance Comparison', fontweight='bold', fontsize=14)
ax1.set_ylim(0, 1.0)
ax1.grid(True, alpha=0.3)

for bar, score in zip(bars, r2_scores):
    height = bar.get_height()
    ax1.text(bar.get_x() + bar.get_width()/2., height + 0.01,
             f'{score:.3f}', ha='center', va='bottom', fontweight='bold')

# 2. Top 15 Drivers by Skill (Car-Independent)
ax2 = axes[0, 1]
top_drivers = driver_analysis.head(15)
colors_drivers = ['red' if 'Hamilton' in name else 'purple' if 'Rosberg' in name 
                 else 'orange' if 'Verstappen' in name else 'lightblue' 
                 for name in top_drivers['driver_name']]

y_pos = np.arange(len(top_drivers))
bars = ax2.barh(y_pos, top_drivers['skill_coefficient'], color=colors_drivers, 
                edgecolor='black', alpha=0.8)
ax2.set_yticks(y_pos)
ax2.set_yticklabels([name.split()[-1] for name in top_drivers['driver_name']], fontsize=10)
ax2.set_xlabel('Skill Coefficient (Lower = Better)', fontsize=12)
ax2.set_title('Top 15 Drivers by Skill (Car-Independent)', fontweight='bold', fontsize=14)
ax2.grid(True, alpha=0.3, axis='x')
ax2.invert_yaxis()

# 3. Constructor Car Performance
ax3 = axes[0, 2]
constructor_colors = ['blue' if 'Mercedes' in name else 'red' if 'Ferrari' in name 
                     else 'orange' if 'Red Bull' in name else 'lightgreen' 
                     for name in constructor_analysis['constructor_name']]

y_pos_const = np.arange(len(constructor_analysis))
bars = ax3.barh(y_pos_const, constructor_analysis['car_coefficient'], 
                color=constructor_colors, edgecolor='black', alpha=0.8)
ax3.set_yticks(y_pos_const)
ax3.set_yticklabels([name[:15] for name in constructor_analysis['constructor_name']], fontsize=9)
ax3.set_xlabel('Car Performance Coefficient (Lower = Better)', fontsize=12)
ax3.set_title('Constructor Car Performance', fontweight='bold', fontsize=14)
ax3.grid(True, alpha=0.3, axis='x')
ax3.invert_yaxis()

# 4. Actual vs Predicted (Full Model)
ax4 = axes[1, 0]
scatter = ax4.scatter(y_test, y_pred_full, alpha=0.6, s=50, 
                     c=y_test, cmap='viridis', edgecolor='black', linewidth=0.5)
ax4.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 
         'r--', lw=2, alpha=0.8)
ax4.set_xlabel('Actual Position', fontsize=12)
ax4.set_ylabel('Predicted Position', fontsize=12)
ax4.set_title(f'Actual vs Predicted (R² = {r2_full:.3f})', fontweight='bold', fontsize=14)
ax4.grid(True, alpha=0.3)

# 5. Residuals vs Predicted
ax5 = axes[1, 1]
residuals = y_test - y_pred_full
ax5.scatter(y_pred_full, residuals, alpha=0.6, s=50, color='blue', edgecolor='black', linewidth=0.5)
ax5.axhline(y=0, color='red', linestyle='--', linewidth=2, alpha=0.8)
ax5.set_xlabel('Predicted Position', fontsize=12)
ax5.set_ylabel('Residuals', fontsize=12)
ax5.set_title('Residuals vs Predicted', fontweight='bold', fontsize=14)
ax5.grid(True, alpha=0.3)

# 6. Driver Skill vs Actual Performance
ax6 = axes[1, 2]
for _, row in driver_analysis.iterrows():
    if 'Hamilton' in row['driver_name']:
        ax6.scatter(row['skill_coefficient'], row['avg_position'], 
                   color='red', s=200, marker='*', edgecolor='black', linewidth=2,
                   label='Hamilton', zorder=5)
    elif 'Rosberg' in row['driver_name']:
        ax6.scatter(row['skill_coefficient'], row['avg_position'], 
                   color='purple', s=200, marker='*', edgecolor='black', linewidth=2,
                   label='Rosberg', zorder=5)
    elif 'Verstappen' in row['driver_name']:
        ax6.scatter(row['skill_coefficient'], row['avg_position'], 
                   color='orange', s=200, marker='*', edgecolor='black', linewidth=2,
                   label='Verstappen', zorder=5)
    else:
        ax6.scatter(row['skill_coefficient'], row['avg_position'], 
                   alpha=0.7, s=60, color='lightblue', edgecolor='black', linewidth=0.5)

ax6.set_xlabel('Skill Coefficient (Car-Independent)', fontsize=12)
ax6.set_ylabel('Average Finishing Position', fontsize=12)
ax6.set_title('Driver Skill vs Actual Performance', fontweight='bold', fontsize=14)
ax6.grid(True, alpha=0.3)
ax6.legend()
ax6.invert_yaxis()

plt.tight_layout(rect=[0, 0.02, 1, 0.94])
plt.show()
            

pit_analysis = pit_stops_2015.merge(
    results_2015[['raceId', 'driverId', 'constructorId']], 
    on=['raceId', 'driverId'], how='left').merge(
    constructors[['constructorId', 'name']], 
    on='constructorId', how='left')
pit_analysis.rename(columns={'name': 'constructor_name'}, inplace=True)

pit_stats = pit_analysis.groupby('constructor_name').agg({
    'duration_seconds': ['mean', 'std', 'count'],
    'stop': 'mean',
    'lap': 'mean'}).round(3)

pit_stats.columns = ['avg_duration', 'std_duration', 'total_stops', 'avg_stops_per_race', 'avg_pit_lap']
pit_stats['consistency_score'] = 1 / (pit_stats['std_duration'] + 0.1)
pit_stats = pit_stats[pit_stats['total_stops'] >= 10].sort_values('avg_duration')  # Filter teams with sufficient data

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('F1 2015 Pit Stop Performance Analysis', fontsize=16, fontweight='bold')

pit_stats_top = pit_stats.head(10)

# Subplot 1: Average pit stop time
bars1 = ax1.bar(range(len(pit_stats_top)), pit_stats_top['avg_duration'], 
                color='red', alpha=0.7, edgecolor='black')
ax1.set_xlabel('Constructor')
ax1.set_ylabel('Average Pit Stop Time (seconds)')
ax1.set_title('Average Pit Stop Duration by Constructor')
ax1.set_xticks(range(len(pit_stats_top)))
ax1.set_xticklabels([name[:8] for name in pit_stats_top.index], rotation=45)
ax1.grid(True, alpha=0.3, axis='y')

# Subplot 2: Pit stop consistency
bars2 = ax2.bar(range(len(pit_stats_top)), pit_stats_top['consistency_score'], 
                color='blue', alpha=0.7, edgecolor='black')
ax2.set_xlabel('Constructor')
ax2.set_ylabel('Consistency Score (Higher = More Consistent)')
ax2.set_title('Pit Stop Consistency by Constructor')
ax2.set_xticks(range(len(pit_stats_top)))
ax2.set_xticklabels([name[:8] for name in pit_stats_top.index], rotation=45)
ax2.grid(True, alpha=0.3, axis='y')

# Subplot 3: Pit stop time distribution (box plot)
all_constructors = pit_stats.index
pit_data_all = []
labels_all = []

for constructor in all_constructors:
    constructor_times = pit_analysis[pit_analysis['constructor_name'] == constructor]['duration_seconds'].dropna()
    if len(constructor_times) > 5:
        pit_data_all.append(constructor_times)
        labels_all.append(constructor[:8])

ax3.boxplot(pit_data_all, labels=labels_all)
ax3.set_ylabel('Pit Stop Duration (seconds)')
ax3.set_title('Pit Stop Time Distribution - All Teams')
ax3.tick_params(axis='x', rotation=45)
ax3.grid(True, alpha=0.3, axis='y')

# Subplot 4: Strategy patterns (scatter plot)
scatter = ax4.scatter(pit_stats['avg_stops_per_race'], pit_stats['avg_duration'], 
                     s=pit_stats['total_stops']/1000, alpha=0.6, c='purple')
ax4.set_xlabel('Average Stops per Race')
ax4.set_ylabel('Average Pit Duration (seconds)')
ax4.set_title('Strategy vs Speed (Size = Total Stops)')

plt.tight_layout()
plt.show()

            

def convert_duration(duration_str):
    try:
        if pd.isna(duration_str) or duration_str == '\\N':
            return np.nan
        duration_str = str(duration_str)
        if ':' in duration_str:
            parts = duration_str.split(':')
            return float(parts[0]) * 60 + float(parts[1])
        return float(duration_str)
    except:
        return np.nan

pit_stops_2015['duration_seconds'] = pit_stops_2015['duration'].apply(convert_duration)
pit_stops_2015 = pit_stops_2015[(pit_stops_2015['duration_seconds'] >= 1) & 
                                (pit_stops_2015['duration_seconds'] <= 60)]
pit_analysis = pit_stops_2015.merge(results_2015[['raceId', 'driverId', 'constructorId']], 
                                   on=['raceId', 'driverId'])
pit_analysis = pit_analysis.merge(constructors[['constructorId', 'name']], on='constructorId')

pit_efficiency = pit_analysis.groupby(['constructorId', 'name']).agg({
    'duration_seconds': ['mean', 'std', 'count'],
    'stop': 'mean'}).round(3)

pit_efficiency.columns = ['_'.join(col) for col in pit_efficiency.columns]
pit_efficiency = pit_efficiency.reset_index()
if len(pit_efficiency) > 0:
    min_time = pit_efficiency['duration_seconds_mean'].min()
    max_time = pit_efficiency['duration_seconds_mean'].max()
    pit_efficiency['pit_stop_efficiency'] = ((max_time - pit_efficiency['duration_seconds_mean']) / 
                                            (max_time - min_time)) * 100
else:
    pit_efficiency['pit_stop_efficiency'] = 0

if len(pit_analysis) > 0:
    strategy_analysis = pit_analysis.groupby(['name', 'raceId', 'driverId']).agg({
        'stop': 'max',
        'lap': ['mean', 'std']}).reset_index()

    strategy_analysis.columns = ['name', 'raceId', 'driverId', 'total_stops', 'avg_pit_lap', 'pit_timing_std']
    strategy_summary = strategy_analysis.groupby('name').agg({
        'total_stops': ['mean', 'std'],
        'avg_pit_lap': 'mean',
        'pit_timing_std': 'mean'}).round(3)
    strategy_summary.columns = ['total_stops_mean', 'total_stops_std', 'avg_pit_lap_mean', 'pit_timing_std_mean']
    strategy_summary = strategy_summary.reset_index()
    strategy_summary['strategy_consistency'] = 1 / (1 + strategy_summary['total_stops_std'].fillna(1))
    strategy_summary['timing_efficiency'] = 1 / (1 + strategy_summary['pit_timing_std_mean'].fillna(1))
    strategy_summary['strategy_efficiency'] = ((strategy_summary['strategy_consistency'] + 
                                             strategy_summary['timing_efficiency']) / 2) * 100
else:
    strategy_summary = pd.DataFrame(columns=['name', 'strategy_efficiency'])

comprehensive_efficiency = constructor_metrics.merge(
    pit_efficiency[['name', 'pit_stop_efficiency']], on='name', how='left').merge(strategy_summary[['name', 'strategy_efficiency']], on='name', how='left')
comprehensive_efficiency['pit_stop_efficiency'] = comprehensive_efficiency['pit_stop_efficiency'].fillna(0)
comprehensive_efficiency['strategy_efficiency'] = comprehensive_efficiency['strategy_efficiency'].fillna(0)
comprehensive_efficiency = comprehensive_efficiency.sort_values('overall_efficiency', ascending=False)

cluster_features = ['overall_efficiency', 'pit_stop_efficiency', 'strategy_efficiency', 'reliability_efficiency']
cluster_data = comprehensive_efficiency[cluster_features].fillna(0)

scaler = StandardScaler()
scaled_features = scaler.fit_transform(cluster_data)
kmeans = KMeans(n_clusters=3, random_state=42)
comprehensive_efficiency['cluster'] = kmeans.fit_predict(scaled_features)
cluster_names = {0: 'Development Teams', 1: 'Mid-tier Teams', 2: 'Operational Leaders'}
comprehensive_efficiency['cluster_name'] = comprehensive_efficiency['cluster'].map(cluster_names)

fig = plt.figure(figsize=(18, 12))
fig.suptitle('2015 F1 - Advanced Statistics & Efficiency Analysis', fontsize=20, fontweight='bold', y=0.98)

# 1. Overall Operational Efficiency
plt.subplot(2, 3, 1)
top_10 = comprehensive_efficiency.head(10)
bars = plt.bar(range(len(top_10)), top_10['overall_efficiency'])

plt.title('Overall Operational Efficiency', fontsize=14, fontweight='bold')
plt.ylabel('Efficiency Score')
plt.xlabel('Constructor Rank')
plt.xticks(range(len(top_10)), [name.replace(' ', '\n') for name in top_10['name']], 
          rotation=45, ha='right', fontsize=10)
plt.ylim(0, 110)
plt.grid(True, alpha = 0.3)

for bar, score in zip(bars, top_10['overall_efficiency']):
   plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1, 
            f'{score:.0f}', ha='center', va='bottom', fontweight='bold')

# 2. Efficiency Components - Top 6 Teams
plt.subplot(2, 3, 2)
top_6 = comprehensive_efficiency.head(6)
x = np.arange(len(top_6))
width = 0.25

bars1 = plt.bar(x - width, top_6['pit_stop_efficiency'], width, label='Pit Stop')
bars2 = plt.bar(x, top_6['reliability_efficiency'], width, label='Operations')
bars3 = plt.bar(x + width, top_6['strategy_efficiency'], width, label='Strategy')

plt.title('Efficiency Components - Top 6 Teams', fontsize=14, fontweight='bold')
plt.ylabel('Efficiency Score')
plt.xlabel('Constructor')
plt.xticks(x, [name.replace(' ', '\n') for name in top_6['name']], fontsize=10)
plt.legend()
plt.ylim(0, 110)
plt.grid(True, alpha = 0.3)

# 3. Operational Clusters
plt.subplot(2, 3, 3)
for cluster_name in comprehensive_efficiency['cluster_name'].unique():
   cluster_data = comprehensive_efficiency[comprehensive_efficiency['cluster_name'] == cluster_name]
   plt.scatter(cluster_data['overall_efficiency'], cluster_data['strategy_efficiency'], 
              s=100, label=cluster_name)

   for _, row in cluster_data.iterrows():
       plt.annotate(row['name'], 
                   (row['overall_efficiency'], row['strategy_efficiency']),
                   xytext=(5, 5), textcoords='offset points', fontsize=8)

plt.title('Operational Clusters', fontsize=14, fontweight='bold')
plt.xlabel('Operations Efficiency')
plt.ylabel('Strategy Efficiency')
plt.legend()
plt.grid(True, alpha = 0.3)

# 4. Operational Efficiency vs Championship Success
plt.subplot(2, 3, 4)
plt.scatter(comprehensive_efficiency['overall_efficiency'], comprehensive_efficiency['points_sum'], s=100)
z = np.polyfit(comprehensive_efficiency['overall_efficiency'], comprehensive_efficiency['points_sum'], 1)
p = np.poly1d(z)
x_trend = np.linspace(comprehensive_efficiency['overall_efficiency'].min(), 
                    comprehensive_efficiency['overall_efficiency'].max(), 100)
plt.plot(x_trend, p(x_trend), "r--", linewidth=2)

for _, row in comprehensive_efficiency.iterrows():
   plt.annotate(row['name'], (row['overall_efficiency'], row['points_sum']),
               xytext=(5, 5), textcoords='offset points', fontsize=9)

plt.title('Operational Efficiency vs Championship Success', fontsize=14, fontweight='bold')
plt.xlabel('Overall Efficiency Score')
plt.ylabel('Championship Points')
plt.grid(True, alpha = 0.3)

# 5. Efficiency Metrics Correlation Matrix
plt.subplot(2, 3, 5)
corr_metrics = ['pit_stop_efficiency', 'reliability_efficiency', 'strategy_efficiency', 'points_sum']
available_metrics = [col for col in corr_metrics if col in comprehensive_efficiency.columns]
correlation_matrix = comprehensive_efficiency[available_metrics].corr()

im = plt.imshow(correlation_matrix, cmap='RdBu_r', vmin=-1, vmax=1)
plt.colorbar(im, shrink=0.8)

for i in range(len(available_metrics)):
   for j in range(len(available_metrics)):
       plt.text(j, i, f'{correlation_matrix.iloc[i, j]:.2f}', 
               ha='center', va='center', fontweight='bold')

plt.title('Efficiency Metrics Correlation Matrix', fontsize=14, fontweight='bold')
metric_labels = ['Pit\nEfficiency', 'Operations\nEfficiency', 'Strategy\nEfficiency', 'Championship\nPoints']
plt.xticks(range(len(available_metrics)), metric_labels[:len(available_metrics)], fontsize=9)
plt.yticks(range(len(available_metrics)), metric_labels[:len(available_metrics)], fontsize=9)

# 6. Pit Excellence vs Success
plt.subplot(2, 3, 6)
sizes = comprehensive_efficiency['reliability_efficiency'] * 3 + 50

plt.scatter(comprehensive_efficiency['pit_stop_efficiency'], 
          comprehensive_efficiency['points_sum'], s=sizes)

significant_teams = comprehensive_efficiency[comprehensive_efficiency['points_sum'] > 30]
for _, row in significant_teams.iterrows():
   plt.annotate(row['name'], (row['pit_stop_efficiency'], row['points_sum']),
               xytext=(5, 5), textcoords='offset points', fontsize=9)

plt.title('Pit Excellence vs Success (Size = Operations)', fontsize=14, fontweight='bold')
plt.xlabel('Pit Stop Efficiency Score')
plt.ylabel('Championship Points')
plt.grid(True, alpha = 0.3)

plt.tight_layout()
plt.subplots_adjust(top=0.93, hspace=0.3, wspace=0.3)
plt.show()
            

comprehensive_analysis = constructor_performance.merge(
    strategy_summary, left_on='name_constructor', right_on='name', how='left'
).merge(
    quali_operations[['constructorId', 'qualifying_operational_efficiency']], 
    on='constructorId', how='left'
).merge(
    utilization_summary[['constructorId', 'operational_reliability']], 
    on='constructorId', how='left'
)

numeric_cols = comprehensive_analysis.select_dtypes(include=[np.number]).columns
comprehensive_analysis[numeric_cols] = comprehensive_analysis[numeric_cols].fillna(
    comprehensive_analysis[numeric_cols].median())
comprehensive_analysis['championship_points_threshold'] = comprehensive_analysis['points_sum'] >= 10  # Meaningful points threshold
comprehensive_analysis['top_6_finish_rate'] = np.where(
    comprehensive_analysis['position_numeric_mean'].notna(),
    (comprehensive_analysis['position_numeric_mean'] <= 6).astype(float),
    0)
leader_avg_position = comprehensive_analysis['position_numeric_mean'].min()
comprehensive_analysis['pace_relative_to_leader'] = np.where(
    comprehensive_analysis['position_numeric_mean'].notna(),
    1 / (1 + (comprehensive_analysis['position_numeric_mean'] - leader_avg_position)),
    0)
def assign_performance_tier(points_sum):
    if points_sum >= 100: 
        return 'Top 4'
    elif points_sum >= 20: 
        return 'midfield'
    else: 
        return 'backmarker'

comprehensive_analysis['performance_tier'] = comprehensive_analysis['points_sum'].apply(assign_performance_tier)
weights = {
    'points_per_race': 0.40, 
    'pace_relative_to_leader': 0.20,
    'championship_points_threshold': 0.10, 
    'top_6_finish_rate': 0.10, 
    'reliability_rate': 0.08,
    'performance_consistency': 0.05,
    'qualifying_operational_efficiency': 0.04, 
    'strategy_consistency': 0.02, 
    'operational_reliability': 0.01 
}
tier_sensitive_metrics = ['performance_consistency', 'strategy_consistency', 'reliability_rate']
championship_metrics = ['points_per_race', 'pace_relative_to_leader', 'top_6_finish_rate']

for metric in championship_metrics:
    if metric in comprehensive_analysis.columns:
        min_val = comprehensive_analysis[metric].min()
        max_val = comprehensive_analysis[metric].max()
        if max_val > min_val:
            comprehensive_analysis[f'{metric}_normalized'] = (
                comprehensive_analysis[metric] - min_val
            ) / (max_val - min_val)
        else:
            comprehensive_analysis[f'{metric}_normalized'] = 0

for metric in tier_sensitive_metrics:
    if metric in comprehensive_analysis.columns:
        comprehensive_analysis[f'{metric}_normalized'] = 0
        for tier in ['top 4', 'midfield', 'backmarker']:
            tier_data = comprehensive_analysis[comprehensive_analysis['performance_tier'] == tier]
            if len(tier_data) > 0:
                min_val = tier_data[metric].min()
                max_val = tier_data[metric].max()
                if max_val > min_val:
                    mask = comprehensive_analysis['performance_tier'] == tier
                    comprehensive_analysis.loc[mask, f'{metric}_normalized'] = (
                        tier_data[metric] - min_val
                    ) / (max_val - min_val)
                else:
                    mask = comprehensive_analysis['performance_tier'] == tier
                    comprehensive_analysis.loc[mask, f'{metric}_normalized'] = 0.5
for metric in ['championship_points_threshold']:
    if metric in comprehensive_analysis.columns:
        comprehensive_analysis[f'{metric}_normalized'] = comprehensive_analysis[metric].astype(float)

comprehensive_analysis['team_efficiency_score'] = sum([
    comprehensive_analysis.get(f'{metric}_normalized', 0) * weight 
    for metric, weight in weights.items()])

comprehensive_analysis['championship_weighted_score'] = (
    comprehensive_analysis['points_per_race'] * 0.6 +
    comprehensive_analysis['pace_relative_to_leader'] * 0.4
)
correlation_metrics = [
    'points_per_race', 'reliability_rate', 'performance_consistency',
    'qualifying_efficiency', 'race_efficiency', 'position_improvement',
    'strategy_consistency', 'pit_efficiency', 'operational_reliability'
]

available_metrics = [col for col in correlation_metrics if col in comprehensive_analysis.columns]
correlation_matrix = comprehensive_analysis[available_metrics].corr()

active_constructors = comprehensive_analysis[comprehensive_analysis['raceId_count'] >= 10]
if len(active_constructors) >= 3:
    scaler = StandardScaler()
    scaled_features = scaler.fit_transform(active_constructors[available_metrics])
    pca = PCA(n_components=min(3, len(available_metrics)))
    pca_features = pca.fit_transform(scaled_features)
    
    for i in range(pca.n_components_):
        active_constructors[f'PC{i+1}'] = pca_features[:, i]
    n_clusters = min(3, len(active_constructors))
    kmeans = KMeans(n_clusters=n_clusters, random_state=42)
    active_constructors['cluster'] = kmeans.fit_predict(scaled_features)

print(f"Creating feature importance analysis...")
print(f"Comprehensive analysis shape: {comprehensive_analysis.shape}")
print(f"Available columns: {list(comprehensive_analysis.columns)}")

core_features = ['reliability_rate', 'position_numeric_mean', 'grid_mean', 'points_count']
feature_cols = [col for col in core_features if col in comprehensive_analysis.columns]

print(f"Selected features: {feature_cols}")

if len(feature_cols) >= 2 and len(comprehensive_analysis) >= 3:
    try:
        X = comprehensive_analysis[feature_cols].fillna(0)
        y = comprehensive_analysis['points_per_race'].fillna(0)
        
        print(f"X shape: {X.shape}, y shape: {y.shape}")
        rf_model = RandomForestRegressor(n_estimators=50, random_state=42)
        rf_model.fit(X, y)
        feature_importance = pd.DataFrame({
            'feature': feature_cols,
            'importance': rf_model.feature_importances_
        }).sort_values('importance', ascending=False)
        
        print(feature_importance)
        
    except Exception as e:
        feature_importance = None
else:
    feature_importance = None
fig = plt.figure(figsize=(18, 12))

# 1. Tier-Based Performance Analysis
plt.subplot(2, 3, 1)
tier_performance = comprehensive_analysis.groupby('performance_tier').agg({
 'points_per_race': 'mean',
 'reliability_rate': 'mean',
 'championship_weighted_score': 'mean'
}).round(3)

x = np.arange(len(tier_performance.index))
width = 0.25

bars1 = plt.bar(x - width, tier_performance['points_per_race'], width, label='Points/Race', color='#2E8B57')
bars2 = plt.bar(x, tier_performance['reliability_rate'], width, label='Reliability', color='#4682B4')
bars3 = plt.bar(x + width, tier_performance['championship_weighted_score'], width, label='Championship Score', color='#CD5C5C')

plt.xlabel('Performance Tier', fontsize=12)
plt.ylabel('Average Score', fontsize=12)
plt.title('Performance Metrics by Championship Tier', fontsize=16, fontweight='bold')
plt.xticks(x, tier_performance.index.str.title(), fontsize=11)
plt.yticks(fontsize=11)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)

# 2. Feature Importance Analysis
plt.subplot(2, 3, 2)
backup_metrics = ['reliability_rate', 'position_numeric_mean', 'grid_mean', 'performance_consistency']
metric_display_names = ['Reliability', 'Avg Position', 'Avg Grid', 'Consistency']
correlations = []

for i, metric in enumerate(backup_metrics):
    if metric in comprehensive_analysis.columns:
        corr = abs(comprehensive_analysis[metric].corr(comprehensive_analysis['points_per_race']))
    if not pd.isna(corr):
        correlations.append({'metric': metric_display_names[i], 'correlation': corr})

if correlations:
    corr_df = pd.DataFrame(correlations).sort_values('correlation', ascending=False)
bars = plt.barh(range(len(corr_df)), corr_df['correlation'], color='#FF6B35')
plt.yticks(range(len(corr_df)), corr_df['metric'], fontsize=11)
plt.xlabel('Absolute Correlation with Points/Race', fontsize=12)
plt.title('Key Performance Predictors', fontsize=16, fontweight='bold')
plt.grid(True, alpha=0.3)

for i, (bar, corr) in enumerate(zip(bars, corr_df['correlation'])):
    plt.text(bar.get_width() + 0.01, bar.get_y() + bar.get_height()/2, 
             f'{corr:.3f}', ha='left', va='center', fontsize=10, fontweight='bold')

# 3. Operational Efficiency Breakdown
plt.subplot(2, 3, 3)
if len(active_constructors) > 0:
 metrics_for_radar = ['reliability_rate', 'performance_consistency', 
                     'qualifying_efficiency', 'race_efficiency']
 available_radar_metrics = [col for col in metrics_for_radar if col in active_constructors.columns]
 
 if len(available_radar_metrics) >= 3:
     top_3_teams = active_constructors.head(3)
     angles = np.linspace(0, 2*np.pi, len(available_radar_metrics), endpoint=False).tolist()
     angles += angles[:1]
     
     ax = plt.subplot(2, 3, 3, projection='polar')
     colors = ['#FF6B6B', '#4ECDC4', '#45B7D1']
     for i, (_, team) in enumerate(top_3_teams.iterrows()):
         values = [team[metric] for metric in available_radar_metrics]
         values += values[:1]
         
         ax.plot(angles, values, linewidth=1, label=team['name_constructor'], color=colors[i], markersize=8)
         ax.fill(angles, values, alpha=0.25, color=colors[i])
     
     clean_radar_labels = [label.replace('_', ' ').title() for label in available_radar_metrics]
     ax.set_thetagrids(np.degrees(angles[:-1]), clean_radar_labels, fontsize=10)
     plt.title('Top 3 Teams - Operational Efficiency', fontsize=16, fontweight='bold', pad=20)
     plt.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0), fontsize=11)

# 4. Pit Stop Performance Distribution
plt.subplot(2, 3, 4)
if not strategy_summary.empty:
 pit_data_clean = strategy_summary.dropna(subset=['duration_seconds_mean'])
 plt.hist(pit_data_clean['duration_seconds_mean'], bins=10, color='#87CEEB', edgecolor='black', alpha=0.7)
 plt.axvline(pit_data_clean['duration_seconds_mean'].mean(), color='red', linestyle='--', linewidth=3,
             label=f'Mean: {pit_data_clean["duration_seconds_mean"].mean():.2f}s')
 plt.xlabel('Average Pit Stop Duration (seconds)', fontsize=12)
 plt.ylabel('Number of Constructors', fontsize=12)
 plt.title('Pit Stop Performance', fontsize=16, fontweight='bold')
 plt.xticks(fontsize=11)
 plt.yticks(fontsize=11)
 plt.legend(fontsize=11)
 plt.grid(True, alpha=0.3)

# 5. Principal Component Analysis
plt.subplot(2, 3, 5)
if len(active_constructors) >= 3 and 'PC1' in active_constructors.columns:
 cluster_colors = ['#E74C3C', '#3498DB', '#2ECC71']
 for cluster in range(n_clusters):
     cluster_data = active_constructors[active_constructors['cluster'] == cluster]
     plt.scatter(cluster_data['PC1'], cluster_data['PC2'], 
                label=f'Cluster {cluster+1}', s=120, color=cluster_colors[cluster], alpha=0.7)
 
 for _, row in active_constructors.iterrows():
     plt.annotate(row['name_constructor'][:8], (row['PC1'], row['PC2']),
                 xytext=(5, 5), textcoords='offset points', fontsize=10)
 
 plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%} variance)', fontsize=12)
 plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%} variance)', fontsize=12)
 plt.title('Constructor Performance Clustering (PCA)', fontsize=16, fontweight='bold')
 plt.xticks(fontsize=11)
 plt.yticks(fontsize=11)
 plt.legend(fontsize=11)
 plt.grid(True, alpha=0.3)

# 6. Position Improvement Analysis
plt.subplot(2, 3, 6)
improvement_data = comprehensive_analysis.dropna(subset=['position_improvement'])
colors = ['#E74C3C' if x < 0 else '#27AE60' for x in improvement_data['position_improvement']]
bars = plt.barh(range(len(improvement_data)), improvement_data['position_improvement'], color=colors)
clean_constructor_names = [name.replace('_', ' ') for name in improvement_data['name_constructor']]
plt.yticks(range(len(improvement_data)), clean_constructor_names, fontsize=10)
plt.xlabel('Average Position Improvement (Grid to Finish)', fontsize=12)
plt.title('Race Performance vs Qualifying', fontsize=16, fontweight='bold')
plt.xticks(fontsize=11)
plt.axvline(x=0, color='black', linestyle='-', linewidth=1)
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
            

distributions = {}
clean_rookies = rookie_df[
    (rookie_df['total_points'] >= 0) & 
    (rookie_df['avg_position'] <= 20) &
    (rookie_df['races_started'] >= 10) & 
    (rookie_df['rookie_year'] <= 2015)].copy(
metrics_to_model = [
    'total_points', 'avg_position', 'points_per_race', 
    'top_10_finishes', 'finish_rate', 'best_position']

for metric in metrics_to_model:
    data = clean_rookies[metric].dropna()
    if len(data) > 5:
        dist_fits = []
        try:
            params = stats.norm.fit(data)
            ks_stat, p_val = stats.kstest(data, lambda x: stats.norm.cdf(x, *params))
            dist_fits.append(('norm', params, p_val))
        except:
            pass
        if (data > 0).all():
            try:
                params = stats.lognorm.fit(data)
                ks_stat, p_val = stats.kstest(data, lambda x: stats.lognorm.cdf(x, *params))
                dist_fits.append(('lognorm', params, p_val))
            except:
                pass
        if (data > 0).all():
            try:
                params = stats.gamma.fit(data)
                ks_stat, p_val = stats.kstest(data, lambda x: stats.gamma.cdf(x, *params))
                dist_fits.append(('gamma', params, p_val))
            except:
                pass
        if dist_fits:
            best_dist = max(dist_fits, key=lambda x: x[2])  # Highest p-value
            distributions[metric] = {
                'distribution': best_dist[0],
                'params': best_dist[1],
                'p_value': best_dist[2]}
        else:
            distributions[metric] = {
                'distribution': 'empirical',
                'data': data.values}
n_simulations = 10000
simulated_seasons = []

for sim in range(n_simulations):
    season = {}
    for metric, dist_info in distributions.items():
        if dist_info['distribution'] == 'norm':
            value = np.random.normal(*dist_info['params'])
        elif dist_info['distribution'] == 'lognorm':
            value = stats.lognorm.rvs(*dist_info['params'])
        elif dist_info['distribution'] == 'gamma':
            value = stats.gamma.rvs(*dist_info['params'])
        elif dist_info['distribution'] == 'empirical':
            value = np.random.choice(dist_info['data'])
        else:
            value = np.nan
        if metric == 'avg_position':
            value = max(1, min(20, value))
        elif metric == 'best_position':
            value = max(1, min(20, value))
        elif metric == 'finish_rate':
            value = max(0, min(1, value))
        elif metric in ['total_points', 'points_per_race', 'top_10_finishes']:
            value = max(0, value)
        season[metric] = value
    simulated_seasons.append(season)
simulated_rookies = pd.DataFrame(simulated_seasons)

verstappen_percentiles = {}
key_metrics = ['total_points', 'avg_position', 'points_per_race', 'top_10_finishes', 'finish_rate']

for metric in key_metrics:
    if metric in max_data and metric in simulated_rookies.columns:
        max_value = max_data[metric]
        if metric == 'avg_position': 
            percentile = (simulated_rookies[metric] >= max_value).mean() * 100
        else: 
            percentile = (simulated_rookies[metric] <= max_value).mean() * 100
        verstappen_percentiles[metric] = percentile
for metric, percentile in verstappen_percentiles.items():
    metric_name = metric.replace('_', ' ').title()
    print(f"   {metric_name:<20}: {percentile:5.1f}th percentile")
historical_percentiles = {}
for metric in key_metrics:
    if metric in max_data:
        max_value = max_data[metric]
        historical_data = rookie_df[metric].dropna()
        if len(historical_data) > 0:
            if metric == 'avg_position':  
                percentile = (historical_data >= max_value).mean() * 100
            else: 
                percentile = (historical_data <= max_value).mean() * 100
            historical_percentiles[metric] = percentile

fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('F1 Rookie Performance: Monte Carlo Simulation vs Verstappen vs Historical Rookies', 
             fontsize=18, fontweight='bold', y=0.96)
verstappen_color = '#FF4500' 
simulation_color = '#4682B4' 
historical_color = '#228B22'
historical_rookies_filtered = rookie_df[rookie_df['rookie_year'] <= 2015]

# 1. Total Points Distribution
ax1 = axes[0, 0]
ax1.hist(simulated_rookies['total_points'], bins=25, alpha=0.5, color=simulation_color, 
         label='Monte Carlo Simulation', density=True, edgecolor='black', linewidth=0.5)
ax1.hist(historical_rookies_filtered['total_points'].dropna(), bins=20, alpha=0.5, color=historical_color, 
         label='Historical Rookies', density=True, edgecolor='black', linewidth=0.5)
if 'total_points' in max_data:
    ax1.axvline(max_data['total_points'], color=verstappen_color, linewidth=4, 
                label=f'Verstappen ({max_data["total_points"]:.0f} points)')
ax1.set_xlabel('Total Championship Points', fontsize=12)
ax1.set_ylabel('Probability Density', fontsize=12)
ax1.set_title('Total Points Distribution', fontsize=14, fontweight='bold', pad=15)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.tick_params(axis='both', which='major', labelsize=10)

# 2. Average Position Distribution
ax2 = axes[0, 1]
ax2.hist(simulated_rookies['avg_position'], bins=25, alpha=0.5, color=simulation_color, 
         label='Monte Carlo Simulation', density=True, edgecolor='black', linewidth=0.5)
ax2.hist(historical_rookies_filtered['avg_position'].dropna(), bins=20, alpha=0.5, color=historical_color, 
         label='Historical Rookies', density=True, edgecolor='black', linewidth=0.5)
if 'avg_position' in max_data:
    ax2.axvline(max_data['avg_position'], color=verstappen_color, linewidth=4, 
                label=f'Verstappen (P{max_data["avg_position"]:.1f} avg)')
ax2.set_xlabel('Average Finishing Position', fontsize=12)
ax2.set_ylabel('Probability Density', fontsize=12)
ax2.set_title('Average Position Distribution', fontsize=14, fontweight='bold', pad=15)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.tick_params(axis='both', which='major', labelsize=10)

# 3. Points Per Race
ax3 = axes[0, 2]
ax3.hist(simulated_rookies['points_per_race'], bins=25, alpha=0.5, color=simulation_color, 
         label='Monte Carlo Simulation', density=True, edgecolor='black', linewidth=0.5)
ax3.hist(historical_rookies_filtered['points_per_race'].dropna(), bins=20, alpha=0.5, color=historical_color, 
         label='Historical Rookies', density=True, edgecolor='black', linewidth=0.5)
if 'points_per_race' in max_data:
    ax3.axvline(max_data['points_per_race'], color=verstappen_color, linewidth=4, 
                label=f'Verstappen ({max_data["points_per_race"]:.2f} pts/race)')
ax3.set_xlabel('Points Per Race', fontsize=12)
ax3.set_ylabel('Probability Density', fontsize=12)
ax3.set_title('Points Per Race Distribution', fontsize=14, fontweight='bold', pad=15)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.tick_params(axis='both', which='major', labelsize=10)

# 4. Top 10 Finishes
ax4 = axes[1, 0]
ax4.hist(simulated_rookies['top_10_finishes'], bins=range(0, 21), alpha=0.5, 
         color=simulation_color, label='Monte Carlo Simulation', density=True, edgecolor='black', linewidth=0.5)
ax4.hist(historical_rookies_filtered['top_10_finishes'].dropna(), bins=range(0, 21), alpha=0.5, 
         color=historical_color, label='Historical Rookies', density=True, edgecolor='black', linewidth=0.5)
if 'top_10_finishes' in max_data:
    ax4.axvline(max_data['top_10_finishes'], color=verstappen_color, linewidth=4, 
                label=f'Verstappen ({max_data["top_10_finishes"]:.0f} top-10s)')
ax4.set_xlabel('Number of Top 10 Finishes', fontsize=12)
ax4.set_ylabel('Probability Density', fontsize=12)
ax4.set_title('Top 10 Finishes Distribution', fontsize=14, fontweight='bold', pad=15)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.tick_params(axis='both', which='major', labelsize=10)

# 5. Verstappen Percentile Rankings
ax5 = axes[1, 1]
metrics = list(verstappen_percentiles.keys())
percentiles = list(verstappen_percentiles.values())
metric_labels = [m.replace('_', ' ').title() for m in metrics]

bars = ax5.bar(range(len(metrics)), percentiles, color=verstappen_color, alpha=0.8, 
               edgecolor='black', linewidth=1.5)
ax5.set_xticks(range(len(metrics)))
ax5.set_xticklabels(metric_labels, rotation=45, ha='right', fontsize=11)
ax5.set_ylabel('Percentile Ranking', fontsize=12)
ax5.set_title('Verstappen vs Monte Carlo Simulation\n(Percentile Rankings)', 
              fontsize=14, fontweight='bold', pad=15)
ax5.set_ylim(0, 100)
ax5.axhline(y=90, color='red', linestyle='--', alpha=0.8, linewidth=2, label='90th Percentile')
ax5.axhline(y=95, color='darkred', linestyle='-', alpha=0.8, linewidth=2, label='95th Percentile')
for bar, pct in zip(bars, percentiles):
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height + 1.5,
             f'{pct:.1f}%', ha='center', va='bottom', fontweight='bold', fontsize=11)

ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3)
ax5.tick_params(axis='both', which='major', labelsize=10)

# 6. Simulation Validation
ax6 = axes[1, 2]
historical_points = historical_rookies_filtered['total_points'].dropna().sort_values()
simulated_points = simulated_rookies['total_points'].sort_values()
if len(historical_points) != len(simulated_points):
    min_len = min(len(historical_points), len(simulated_points))
    historical_points = np.interp(np.linspace(0, 1, min_len), 
                                 np.linspace(0, 1, len(historical_points)), historical_points)
    simulated_points = np.interp(np.linspace(0, 1, min_len), 
                                np.linspace(0, 1, len(simulated_points)), simulated_points)

ax6.scatter(historical_points, simulated_points, alpha=0.7, s=30, 
           color=simulation_color, edgecolor='black', linewidth=0.5)
min_val = min(historical_points.min(), simulated_points.min())
max_val = max(historical_points.max(), simulated_points.max())
ax6.plot([min_val, max_val], [min_val, max_val], 'r--', lw=3, alpha=0.8, 
         label='Perfect Correlation')
correlation = np.corrcoef(historical_points, simulated_points)[0, 1]
ax6.text(0.05, 0.95, f'Correlation: {correlation:.3f}', 
         transform=ax6.transAxes, fontsize=11, fontweight='bold',
         bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

ax6.set_xlabel('Historical Points', fontsize=12)
ax6.set_ylabel('Simulated Points', fontsize=12)
ax6.set_title('Simulation Validation\n(Q-Q Plot)', fontsize=14, fontweight='bold', pad=15)
ax6.legend(fontsize=10)
ax6.grid(True, alpha=0.3)
ax6.tick_params(axis='both', which='major', labelsize=10)

plt.tight_layout(rect=[0, 0.02, 1, 0.94])
plt.subplots_adjust(hspace=0.35, wspace=0.25)
plt.show()