import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

# Function to simulate binomial distribution and normal approximation
def binomial_vs_normal(N, p=0.5):
    # Simulating the binomial distribution
    binomial_dist = np.random.binomial(N, p, 10000)
    
    # Calculating mean and standard deviation for the normal approximation
    mean = N * p
    std_dev = np.sqrt(N * p * (1 - p))
    
    # Generating values for the normal approximation
    normal_dist = np.random.normal(mean, std_dev, 10000)
    
    # Plotting both distributions
    plt.figure(figsize=(10, 6))
    
    # Histogram for binomial distribution
    plt.hist(binomial_dist, bins=30, density=True, alpha=0.6, color='blue', label='Binomial Distribution')
    
    # Histogram for normal approximation
    plt.hist(normal_dist, bins=30, density=True, alpha=0.6, color='green', label='Normal Approximation')
    
    # Overlay normal distribution curve for better comparison
    x = np.linspace(min(binomial_dist), max(binomial_dist), 100)
    plt.plot(x, stats.norm.pdf(x, mean, std_dev), 'k', linewidth=2, label='Normal PDF')
    
    plt.title(f"Binomial(N={N}, p={p}) vs. Normal Approximation")
    plt.legend(loc='best')
    plt.xlabel('Number of Successes')
    plt.ylabel('Density')
    plt.show()

# Simulate with N=20, N=50, and N=100
binomial_vs_normal(20)
binomial_vs_normal(50)
binomial_vs_normal(10000000)