Expectation and Variance

import numpy as np
from scipy import stats

# === Theoretical Values ===
# Normal distribution: mu=5, sigma=2
mu, sigma = 5, 2
print(f"E[X] = {mu}")
print(f"Var(X) = {sigma**2}")

# === Empirical Verification via Sampling ===
np.random.seed(42)
samples = np.random.normal(mu, sigma, size=100_000)
print(f"Sample mean: {samples.mean():.4f}")
print(f"Sample variance: {samples.var():.4f}")
print(f"Sample std dev: {samples.std():.4f}")

# === Discrete Distribution Moments ===
# Fair die
die = np.arange(1, 7)
prob = np.ones(6) / 6
E_die = np.sum(die * prob)
E_die2 = np.sum(die**2 * prob)
Var_die = E_die2 - E_die**2
print(f"E[die] = {E_die:.4f}, Var(die) = {Var_die:.4f}")

# === Custom Random Variable ===
def compute_expectation(values, probs):
    """Compute E[X] for a discrete random variable."""
    return np.sum(values * probs)

def compute_variance(values, probs):
    """Compute Var(X) for a discrete random variable."""
    mean = compute_expectation(values, probs)
    return compute_expectation((values - mean)**2, probs)

values = np.array([0, 1, 2, 3, 4])
probs = np.array([0.1, 0.2, 0.3, 0.25, 0.15])
print(f"E[X] = {compute_expectation(values, probs):.4f}")
print(f"Var(X) = {compute_variance(values, probs):.4f}")

# === Moment Generating Function ===
def mgf_normal(t, mu, sigma):
    """MGF of N(mu, sigma^2)."""
    return np.exp(mu * t + 0.5 * sigma**2 * t**2)

t = 0.1
print(f"M_X({t}) = {mgf_normal(t, mu, sigma):.6f}")

# === Empirical MGF ===
empirical_mgf = np.mean(np.exp(samples * t))
print(f"Empirical M_X({t}) = {empirical_mgf:.6f}")

# === Higher Moments with Scipy ===
from scipy.stats import skew, kurtosis
print(f"Skewness: {skew(samples):.4f}")
print(f"Excess Kurtosis: {kurtosis(samples):.4f}")