Coin Flip Simulator

import random
import math
import matplotlib.pyplot as plt
import numpy as np
from collections import Counter

def Nchoosek(N, k):
    return math.factorial(N) / (math.factorial(k) * math.factorial(N - k))

def binomial(n,k, p):
    return Nchoosek(n, k) * math.pow(p, k) * math.pow(1-p, n-k)


def run_experiment(n_flips, n_experiments, p_heads=0.5):
    """
    Simulate coin flips and count number of heads.

    :param n_flips: Number of coin flips per experiment.
    :param n_experiments: Number of experiments to run.
    :param p_heads: Probability of getting heads (default 0.5 for a fair coin)

    :return: List of head counts for each experiment.
    """
    results = []

    for experiment in range(n_experiments):
        flips = [random.random() < p_heads for _ in range(n_flips)]
        heads_count = sum(flips)
        results.append(heads_count)

    return results


def compare_theoretical_vs_experimental():
    """Compare theoretical binomial distribution with experimental results."""
    # Parameters
    n_flips = 10          # Number of coin flips per experiment
    n_experiments = 1000  # Number of experiments to run
    p_heads = 0.5         # Probability of heads (fair coin)

    # Run the experiment
    print(f"Running {n_experiments} experiments with {n_flips} coin flips each...")
    results = run_experiment(n_flips, n_experiments, p_heads)

    # Count occurrences of each number of heads
    result_counts = Counter(results)

    # Convert to probabilities
    experimental_probs = {k: v/n_experiments for k, v in result_counts.items()}

    # Calculate theoretical probabilities using your binomial function
    theoretical_probs = {k: binomial(n_flips, k, p_heads) for k in range(n_flips + 1)}

    # Prepare data for plotting
    k_values = list(range(n_flips + 1))
    theoretical_values = [theoretical_probs[k] for k in k_values]
    experimental_values = [experimental_probs.get(k, 0) for k in k_values]

    # Plot the results
    plt.figure(figsize=(12, 6))

    # Bar width and positions
    width = 0.35
    x = np.array(k_values)

    # Create bars
    plt.bar(x - width/2, theoretical_values, width, label='Theoretical Probability')
    plt.bar(x + width/2, experimental_values, width, label='Experimental Probability')

    # Add labels and title
    plt.xlabel('Number of Heads')
    plt.ylabel('Probability')
    plt.title(f'Binomial Distribution: {n_flips} coin flips, p(heads)={p_heads}')
    plt.xticks(k_values)
    plt.legend()
    plt.grid(axis='y', linestyle='--', alpha=0.7)

    # Add a text box with statistics
    mean_theoretical = n_flips * p_heads
    var_theoretical = n_flips * p_heads * (1 - p_heads)
    mean_experimental = sum(k * count for k, count in result_counts.items()) / n_experiments

    stats_text = f"Theoretical Mean: {mean_theoretical}\n"
    stats_text += f"Theoretical Variance: {var_theoretical}\n"
    stats_text += f"Experimental Mean: {mean_experimental:.2f}"

    plt.figtext(0.15, 0.8, stats_text, bbox=dict(facecolor='white', alpha=0.8))

    # Show plot
    plt.tight_layout()
    #plt.show()
    plt.savefig('coinflip.png')

    # Print summary
    print("\nSummary:")
    print(f"Theoretical Mean: {mean_theoretical}")
    print(f"Theoretical Variance: {var_theoretical}")
    print(f"Experimental Mean: {mean_experimental:.2f}")

    # Print probabilities table
    print("\nProbabilities Table:")
    print("Heads | Theoretical | Experimental | Difference")
    print("-" * 50)
    for k in k_values:
        theo = theoretical_probs[k]
        exp = experimental_probs.get(k, 0)
        diff = exp - theo
        print(f"{k:5d} | {theo:.6f} | {exp:.6f} | {diff:.6f}")

compare_theoretical_vs_experimental()