import numpy as np
import matplotlib.pyplot as plt

# Parameters of the normal distribution
mu = 0
sigma = 4
n = 1000

np.random.seed(42)

# Generate random sample
samples = np.random.normal(mu, sigma, n)

# Sample statistics
sample_mean = np.mean(samples)
sample_std = np.std(samples)

print(f"Sample mean: {sample_mean:.4f}")
print(f"Sample standard deviation: {sample_std:.4f}")

# X-axis for theoretical curve
x = np.linspace(mu - 4*sigma, mu + 4*sigma, 1000)

# Theoretical probability density function
pdf = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((x - mu) / sigma)**2)

# Figure setup
plt.figure(figsize=(8, 5))

# Histogram of random sample
plt.hist(samples, bins=30, density=True,
         alpha=0.6, color='lightgray', edgecolor='black',
         label='Sample')

# Theoretical curve
plt.plot(x, pdf, color='black', linewidth=2, label=r'Theoretical')

# Sample mean
plt.axvline(sample_mean, color='blue', linestyle='--', linewidth=2,
            label='Mean (sample)')

# Theoretical mean
plt.axvline(mu, color='red', linestyle=':', linewidth=2,
            label='Mean (theory)')

# Theoretical ±1σ (shaded area)
#plt.fill_between(x, 0, pdf, where=(x >= mu - sigma) & (x <= mu + sigma),
#                 color='red', alpha=0.1, label=r'±1σ')

# Font sizes
plt.title('Normal Distribution – Random Sample vs. Theory', fontsize=20)
plt.xlabel('x', fontsize=18)
plt.ylabel('Density', fontsize=18)

# Legend at upper right
plt.legend(fontsize=16, loc='upper right')

# Customize ticks
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.locator_params(axis='x', nbins=8)
plt.locator_params(axis='y', nbins=6)

# Tight layout
plt.tight_layout()

# Save as SVG for Wikipedia
plt.savefig("normal_distribution_wikipedia_tics_legend.svg", format="svg")

plt.show()