distribution_analysis.py

"""
looking at how the disaster and GDP data is distributed
using histograms, boxplots, scatter plots etc
"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

# make it look nice
sns.set_style("whitegrid")
plt.rcParams['figure.figsize'] = (16, 10)

print("="*80)
print("DISTRIBUTION ANALYSIS: NATURAL DISASTERS & GDP DATA")
print("="*80)
print()

# load everything
print("Loading datasets...")
disasters = pd.read_excel('naturalDisasters.xlsx')
gdp_growth = pd.read_csv('gdp_growth.csv')
gdp = pd.read_csv('gdp.csv')
gdp_per_capita = pd.read_csv('gdp_per_capita.csv')

print(f"✓ Disasters: {disasters.shape}")
print(f"✓ GDP Growth: {gdp_growth.shape}")
print(f"✓ GDP: {gdp.shape}")
print(f"✓ GDP Per Capita: {gdp_per_capita.shape}")
print()

# clean up the disaster data
disasters['Start Year'] = pd.to_numeric(disasters['Start Year'], errors='coerce')
disasters['Total Deaths'] = pd.to_numeric(disasters['Total Deaths'], errors='coerce')
disasters['Total Affected'] = pd.to_numeric(disasters['Total Affected'], errors='coerce')
disasters['Total Damage USD'] = pd.to_numeric(disasters['Total Damage (\'000 US$)'], errors='coerce') / 1000

# only keep years we care about
disasters_filtered = disasters[(disasters['Start Year'] >= 1960) & (disasters['Start Year'] <= 2020)]

# get GDP ready
year_cols = [col for col in gdp_growth.columns if col.isdigit()]

# group by year
disasters_by_year = disasters_filtered.groupby('Start Year').agg({
    'Disaster Type': 'count',
    'Total Deaths': 'sum',
    'Total Affected': 'sum',
    'Total Damage USD': 'sum'
}).reset_index()
disasters_by_year.columns = ['Year', 'Total_Disasters', 'Total_Deaths', 'Total_Affected', 'Economic_Damage_Million_USD']

# calculate global GDP stuff
global_gdp_growth = []
global_gdp_total = []

for year_col in year_cols:
    year = int(year_col)
    if 1960 <= year <= 2020:
        gdp_growth_values = gdp_growth[year_col].dropna()
        gdp_values = gdp[year_col].dropna()
        
        global_gdp_growth.append({
            'Year': year,
            'Global_GDP_Growth': gdp_growth_values.mean(),
            'GDP_Growth_Std': gdp_growth_values.std(),
            'GDP_Growth_Median': gdp_growth_values.median()
        })
        
        global_gdp_total.append({
            'Year': year,
            'Global_GDP_Billion': gdp_values.sum() / 1e9
        })

gdp_growth_df = pd.DataFrame(global_gdp_growth)
gdp_total_df = pd.DataFrame(global_gdp_total)

# combine everything together
merged_data = disasters_by_year.merge(gdp_growth_df, on='Year', how='outer')
merged_data = merged_data.merge(gdp_total_df, on='Year', how='outer')
merged_data = merged_data.sort_values('Year').reset_index(drop=True)

# if there were no disasters that year, just use 0
merged_data['Total_Disasters'] = merged_data['Total_Disasters'].fillna(0)
merged_data['Total_Deaths'] = merged_data['Total_Deaths'].fillna(0)
merged_data['Total_Affected'] = merged_data['Total_Affected'].fillna(0)
merged_data['Economic_Damage_Million_USD'] = merged_data['Economic_Damage_Million_USD'].fillna(0)

print("="*80)
print("SUMMARY STATISTICS")
print("="*80)
print()

# variables we want to look at
variables = {
    'Total_Disasters': 'Number of Natural Disasters',
    'Total_Deaths': 'Deaths from Disasters',
    'Total_Affected': 'People Affected',
    'Economic_Damage_Million_USD': 'Economic Damage (Million USD)',
    'Global_GDP_Growth': 'Global GDP Growth Rate (%)',
    'Global_GDP_Billion': 'Global GDP (Billion USD)'
}

summary_stats = []

for var, description in variables.items():
    if var in merged_data.columns:
        data = merged_data[var].dropna()
        
        stats_dict = {
            'Variable': description,
            'Count': len(data),
            'Mean': data.mean(),
            'Median': data.median(),
            'Std Dev': data.std(),
            'Min': data.min(),
            'Max': data.max(),
            'Q1 (25%)': data.quantile(0.25),
            'Q3 (75%)': data.quantile(0.75),
            'IQR': data.quantile(0.75) - data.quantile(0.25),
            'Skewness': stats.skew(data),
            'Kurtosis': stats.kurtosis(data)
        }
        summary_stats.append(stats_dict)

summary_df = pd.DataFrame(summary_stats)
print(summary_df.to_string(index=False))
print()

# what do the distributions look like?
print("="*80)
print("DISTRIBUTION SHAPE ANALYSIS")
print("="*80)
print()

for var, description in variables.items():
    if var in merged_data.columns:
        data = merged_data[var].dropna()
        skew = stats.skew(data)
        kurt = stats.kurtosis(data)
        
        print(f"{description}:")
        print(f"  Skewness: {skew:.3f} - ", end="")
        if abs(skew) < 0.5:
            print("Approximately symmetric")
        elif skew > 0:
            print("Right-skewed (positive skew) - tail extends to right")
        else:
            print("Left-skewed (negative skew) - tail extends to left")
        
        print(f"  Kurtosis: {kurt:.3f} - ", end="")
        if abs(kurt) < 0.5:
            print("Normal-like distribution")
        elif kurt > 0:
            print("Heavy tails (leptokurtic) - more outliers than normal")
        else:
            print("Light tails (platykurtic) - fewer outliers than normal")
        print()

# check if anything looks normal
print("="*80)
print("NORMALITY TESTS (Shapiro-Wilk)")
print("="*80)
print()

for var, description in variables.items():
    if var in merged_data.columns:
        data = merged_data[var].dropna()
        if len(data) >= 3:
            stat, p_value = stats.shapiro(data)
            is_normal = "Yes" if p_value > 0.05 else "No"
            print(f"{description:40s}: p={p_value:.4f} - Normal? {is_normal}")
print()

print("="*80)
print("GENERATING DISTRIBUTION VISUALIZATIONS")
print("="*80)
print()

# make the big distribution figure
fig = plt.figure(figsize=(20, 24))
gs = fig.add_gridspec(6, 3, hspace=0.4, wspace=0.3)

plot_configs = [
    ('Total_Disasters', 'Number of Disasters', 'tab:red'),
    ('Total_Deaths', 'Deaths', 'tab:orange'),
    ('Total_Affected', 'People Affected', 'tab:purple'),
    ('Economic_Damage_Million_USD', 'Economic Damage (Million USD)', 'tab:olive'),
    ('Global_GDP_Growth', 'GDP Growth Rate (%)', 'tab:green'),
    ('Global_GDP_Billion', 'Global GDP (Billion USD)', 'tab:blue')
]

for idx, (var, label, color) in enumerate(plot_configs):
    data = merged_data[var].dropna()
    
    # histogram on the left
    ax_hist = fig.add_subplot(gs[idx, 0])
    ax_hist.hist(data, bins=20, color=color, alpha=0.7, edgecolor='black', linewidth=1.2)
    ax_hist.axvline(data.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {data.mean():.2f}')
    ax_hist.axvline(data.median(), color='blue', linestyle='--', linewidth=2, label=f'Median: {data.median():.2f}')
    ax_hist.set_xlabel(label, fontsize=10, fontweight='bold')
    ax_hist.set_ylabel('Frequency', fontsize=10)
    ax_hist.set_title(f'Histogram: {label}', fontsize=11, fontweight='bold')
    ax_hist.legend(fontsize=9)
    ax_hist.grid(True, alpha=0.3)
    
    # boxplot in the middle
    ax_box = fig.add_subplot(gs[idx, 1])
    bp = ax_box.boxplot(data, vert=True, patch_artist=True, widths=0.6,
                         boxprops=dict(facecolor=color, alpha=0.7),
                         medianprops=dict(color='red', linewidth=2),
                         whiskerprops=dict(linewidth=1.5),
                         capprops=dict(linewidth=1.5))
    ax_box.set_ylabel(label, fontsize=10, fontweight='bold')
    ax_box.set_title(f'Boxplot: {label}', fontsize=11, fontweight='bold')
    ax_box.grid(True, alpha=0.3, axis='y')
    
    # count how many outliers
    q1, q3 = data.quantile(0.25), data.quantile(0.75)
    iqr = q3 - q1
    lower_bound = q1 - 1.5 * iqr
    upper_bound = q3 + 1.5 * iqr
    outliers = data[(data < lower_bound) | (data > upper_bound)]
    ax_box.text(0.5, 0.98, f'Outliers: {len(outliers)}', transform=ax_box.transAxes,
                fontsize=9, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    # Q-Q plot on the right
    ax_qq = fig.add_subplot(gs[idx, 2])
    stats.probplot(data, dist="norm", plot=ax_qq)
    ax_qq.set_title(f'Q-Q Plot: {label}', fontsize=11, fontweight='bold')
    ax_qq.grid(True, alpha=0.3)

plt.suptitle('Distribution Analysis: All Key Variables', fontsize=16, fontweight='bold', y=0.995)
plt.savefig('distribution_analysis.png', dpi=300, bbox_inches='tight')
print("✓ Distribution analysis saved: distribution_analysis.png")

# now make a scatter plot matrix
print()
print("Generating scatter plot matrix...")
fig2, axes = plt.subplots(3, 3, figsize=(18, 16))

scatter_vars = [
    ('Total_Disasters', 'Disasters'),
    ('Total_Deaths', 'Deaths'),
    ('Economic_Damage_Million_USD', 'Damage ($M)')
]

for i in range(3):
    for j in range(3):
        ax = axes[i, j]
        
        if i == j:
            # histograms down the diagonal
            var, label = scatter_vars[i]
            data = merged_data[var].dropna()
            ax.hist(data, bins=20, color='steelblue', alpha=0.7, edgecolor='black')
            ax.set_ylabel('Frequency', fontsize=10)
            ax.set_title(label, fontsize=12, fontweight='bold')
            ax.grid(True, alpha=0.3)
        else:
            # scatter plots everywhere else
            var_x, label_x = scatter_vars[j]
            var_y, label_y = scatter_vars[i]
            
            x_data = merged_data[var_x].dropna()
            y_data = merged_data[var_y].dropna()
            
            # line up the data
            common_idx = merged_data[[var_x, var_y]].dropna().index
            x = merged_data.loc[common_idx, var_x]
            y = merged_data.loc[common_idx, var_y]
            
            ax.scatter(x, y, alpha=0.6, s=80, c=merged_data.loc[common_idx, 'Year'], 
                      cmap='viridis', edgecolors='black', linewidth=0.5)
            
            # show the correlation
            if len(x) > 2:
                corr = x.corr(y)
                ax.text(0.05, 0.95, f'r = {corr:.3f}', transform=ax.transAxes,
                       fontsize=10, verticalalignment='top',
                       bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
            
            if i == 2:
                ax.set_xlabel(label_x, fontsize=10, fontweight='bold')
            if j == 0:
                ax.set_ylabel(label_y, fontsize=10, fontweight='bold')
            
            ax.grid(True, alpha=0.3)

plt.suptitle('Scatter Plot Matrix: Disaster Variables', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('scatter_matrix.png', dpi=300, bbox_inches='tight')
print("✓ Scatter plot matrix saved: scatter_matrix.png")

# more scatter plots for disaster vs GDP
print()
print("Generating disaster-GDP relationship plots...")
fig3, axes = plt.subplots(2, 2, figsize=(16, 12))

# first plot
ax1 = axes[0, 0]
scatter1 = ax1.scatter(merged_data['Total_Disasters'], merged_data['Global_GDP_Growth'],
                       c=merged_data['Year'], cmap='plasma', s=120, alpha=0.7, 
                       edgecolors='black', linewidth=1)
ax1.set_xlabel('Number of Disasters', fontsize=12, fontweight='bold')
ax1.set_ylabel('GDP Growth Rate (%)', fontsize=12, fontweight='bold')
ax1.set_title('Disasters vs GDP Growth', fontsize=13, fontweight='bold')
ax1.grid(True, alpha=0.3)
plt.colorbar(scatter1, ax=ax1, label='Year')

# add the correlation number
corr1 = merged_data[['Total_Disasters', 'Global_GDP_Growth']].corr().iloc[0, 1]
ax1.text(0.05, 0.95, f'Correlation: {corr1:.4f}', transform=ax1.transAxes,
         fontsize=11, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.7))

# second plot
ax2 = axes[0, 1]
scatter2 = ax2.scatter(merged_data['Total_Deaths'], merged_data['Global_GDP_Growth'],
                       c=merged_data['Year'], cmap='plasma', s=120, alpha=0.7,
                       edgecolors='black', linewidth=1)
ax2.set_xlabel('Total Deaths', fontsize=12, fontweight='bold')
ax2.set_ylabel('GDP Growth Rate (%)', fontsize=12, fontweight='bold')
ax2.set_title('Deaths vs GDP Growth', fontsize=13, fontweight='bold')
ax2.grid(True, alpha=0.3)
plt.colorbar(scatter2, ax=ax2, label='Year')

corr2 = merged_data[['Total_Deaths', 'Global_GDP_Growth']].corr().iloc[0, 1]
ax2.text(0.05, 0.95, f'Correlation: {corr2:.4f}', transform=ax2.transAxes,
         fontsize=11, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.7))

# third plot
ax3 = axes[1, 0]
scatter3 = ax3.scatter(merged_data['Economic_Damage_Million_USD'], merged_data['Global_GDP_Growth'],
                       c=merged_data['Year'], cmap='plasma', s=120, alpha=0.7,
                       edgecolors='black', linewidth=1)
ax3.set_xlabel('Economic Damage (Million USD)', fontsize=12, fontweight='bold')
ax3.set_ylabel('GDP Growth Rate (%)', fontsize=12, fontweight='bold')
ax3.set_title('Economic Damage vs GDP Growth', fontsize=13, fontweight='bold')
ax3.grid(True, alpha=0.3)
plt.colorbar(scatter3, ax=ax3, label='Year')

corr3 = merged_data[['Economic_Damage_Million_USD', 'Global_GDP_Growth']].corr().iloc[0, 1]
ax3.text(0.05, 0.95, f'Correlation: {corr3:.4f}', transform=ax3.transAxes,
         fontsize=11, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.7))

# fourth plot
ax4 = axes[1, 1]
scatter4 = ax4.scatter(merged_data['Total_Disasters'], merged_data['Global_GDP_Billion'],
                       c=merged_data['Year'], cmap='plasma', s=120, alpha=0.7,
                       edgecolors='black', linewidth=1)
ax4.set_xlabel('Number of Disasters', fontsize=12, fontweight='bold')
ax4.set_ylabel('Global GDP (Billion USD)', fontsize=12, fontweight='bold')
ax4.set_title('Disasters vs Global GDP', fontsize=13, fontweight='bold')
ax4.grid(True, alpha=0.3)
plt.colorbar(scatter4, ax=ax4, label='Year')

corr4 = merged_data[['Total_Disasters', 'Global_GDP_Billion']].corr().iloc[0, 1]
ax4.text(0.05, 0.95, f'Correlation: {corr4:.4f}', transform=ax4.transAxes,
         fontsize=11, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.7))

plt.suptitle('Disaster-GDP Relationship Analysis', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('disaster_gdp_scatter.png', dpi=300, bbox_inches='tight')
print("✓ Disaster-GDP scatter plots saved: disaster_gdp_scatter.png")

# make a correlation heatmap
print()
print("Generating correlation heatmap...")
fig4, ax = plt.subplots(figsize=(12, 10))

corr_vars = ['Total_Disasters', 'Total_Deaths', 'Total_Affected', 
             'Economic_Damage_Million_USD', 'Global_GDP_Growth', 'Global_GDP_Billion']
corr_matrix = merged_data[corr_vars].corr()

mask = np.triu(np.ones_like(corr_matrix, dtype=bool), k=1)
sns.heatmap(corr_matrix, annot=True, fmt='.3f', cmap='RdYlGn', center=0,
            square=True, linewidths=2, cbar_kws={"shrink": 0.8},
            annot_kws={'size': 11, 'weight': 'bold'}, mask=mask, ax=ax)

ax.set_title('Correlation Matrix: All Variables', fontsize=15, fontweight='bold', pad=20)
labels = ['Disasters', 'Deaths', 'Affected', 'Damage ($M)', 'GDP Growth (%)', 'GDP ($B)']
ax.set_xticklabels(labels, rotation=45, ha='right', fontsize=11)
ax.set_yticklabels(labels, rotation=0, fontsize=11)

plt.tight_layout()
plt.savefig('correlation_matrix_full.png', dpi=300, bbox_inches='tight')
print("✓ Correlation matrix saved: correlation_matrix_full.png")