DEV Community

Senthil Kumaran
Senthil Kumaran

Posted on • Originally published at senthil.learntosolveit.com

Coin Flip Simulator

This is a coin flip simulator. It compare theoretical binomial distribution with experimental results

Image description


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()
Enter fullscreen mode Exit fullscreen mode

Top comments (0)