This is a coin flip simulator. It compare theoretical binomial distribution with experimental results
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()
Top comments (0)