Common Probability Distributions

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

import numpy as np
from scipy import stats
 
# Standard normal distribution
standard_normal = stats.norm(loc=0, scale=1)
 
# Key quantiles
print("Standard Normal Distribution Z ~ N(0,1)")
print("=" * 50)
 
# Probability within k standard deviations
for k in [1, 2, 3, 4]:
    prob = standard_normal.cdf(k) - standard_normal.cdf(-k)
    print(f"P(-{k} < Z < {k}) = {prob:.6f} ({prob*100:.3f}%)")
 
print()
 
# Critical values for confidence intervals
confidence_levels = [0.90, 0.95, 0.99]
print("Critical z-values for two-tailed confidence intervals:")
for cl in confidence_levels:
    alpha = 1 - cl
    z_critical = standard_normal.ppf(1 - alpha/2)
    print(f"  {cl*100:.0f}% CI: z = ±{z_critical:.4f}")
 
print()
 
# Converting between X ~ N(μ,σ²) and Z ~ N(0,1)
mu, sigma = 100, 15  # Example: IQ scores
x = 130  # A specific IQ score
 
z = (x - mu) / sigma  # Standardize
prob_below = standard_normal.cdf(z)
print(f"IQ scores: X ~ N({mu}, {sigma}²)")
print(f"IQ = {x} → z-score = {z:.2f}")
print(f"P(IQ < {x}) = P(Z < {z:.2f}) = {prob_below:.4f}")
print(f"Percentile: {prob_below*100:.1f}th")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
 
def demonstrate_clt(original_dist, dist_name, n_samples=10000):
    """
    Demonstrate CLT by showing how sample means converge to Gaussian.
    """
    sample_sizes = [1, 2, 5, 10, 30, 100]
    fig, axes = plt.subplots(2, 3, figsize=(14, 8))
    
    for ax, n in zip(axes.flatten(), sample_sizes):
        # Generate n_samples sets of n observations each
        data = original_dist.rvs(size=(n_samples, n))
        sample_means = data.mean(axis=1)
        
        # True parameters of sample mean distribution
        true_mean = original_dist.mean()
        true_std = original_dist.std() / np.sqrt(n)
        
        # Plot histogram of sample means
        ax.hist(sample_means, bins=50, density=True, alpha=0.7, 
                color='steelblue', edgecolor='navy')
        
        # Overlay the normal approximation
        x = np.linspace(sample_means.min(), sample_means.max(), 200)
        normal_pdf = stats.norm.pdf(x, true_mean, true_std)
        ax.plot(x, normal_pdf, 'r-', linewidth=2, label='Normal approx.')
        
        ax.set_title(f'n = {n}')
        ax.set_xlabel('Sample Mean')
        ax.legend()
    
    plt.suptitle(f'CLT Demonstration: {dist_name}', fontsize=14)
    plt.tight_layout()
    plt.show()
 
# Demonstrate with highly non-normal distributions
 
# 1. Exponential (skewed right)
print("Exponential Distribution (λ=1):")
exponential = stats.expon(scale=1)
print(f"  Mean: {exponential.mean():.2f}, Variance: {exponential.var():.2f}")
print(f"  Skewness: {exponential.stats(moments='s')[0]:.2f} (highly skewed)")
# demonstrate_clt(exponential, "Exponential(λ=1)")
 
# 2. Uniform (symmetric, non-normal)
print("
Uniform[0,1] Distribution:")
uniform = stats.uniform(0, 1)
print(f"  Mean: {uniform.mean():.2f}, Variance: {uniform.var():.4f}")
# demonstrate_clt(uniform, "Uniform[0,1]")
 
# 3. Bernoulli (discrete)
print("
Bernoulli(p=0.3) Distribution:")
bernoulli = stats.bernoulli(0.3)
print(f"  Mean: {bernoulli.mean():.2f}, Variance: {bernoulli.var():.4f}")
 
# Manual CLT demonstration
np.random.seed(42)
n = 50  # Sample size
n_experiments = 10000
 
# Draw n_experiments sample means, each from n observations
bernoulli_means = np.random.binomial(n, 0.3, n_experiments) / n
 
# Compare to normal
mu = 0.3
sigma = np.sqrt(0.3 * 0.7 / n)
print(f"
For n=50 Bernoulli samples:")
print(f"  μ_X̄ = {mu:.3f}, σ_X̄ = {sigma:.4f}")
print(f"  Observed: mean = {bernoulli_means.mean():.4f}, std = {bernoulli_means.std():.4f}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

import numpy as np
from scipy import stats
 
def gaussian_mle(samples: np.ndarray) -> dict:
    """
    Maximum Likelihood Estimation for Gaussian parameters.
    
    Returns both MLE estimates and unbiased estimates with confidence intervals.
    """
    n = len(samples)
    
    # MLE estimates
    mu_mle = np.mean(samples)
    var_mle = np.var(samples, ddof=0)  # ddof=0 for MLE (divide by n)
    
    # Unbiased estimates  
    var_unbiased = np.var(samples, ddof=1)  # ddof=1 for unbiased (divide by n-1)
    sigma_unbiased = np.std(samples, ddof=1)
    
    # Standard error of the mean
    se_mean = sigma_unbiased / np.sqrt(n)
    
    # Confidence intervals for μ using t-distribution
    t_critical = stats.t.ppf(0.975, df=n-1)  # 95% CI
    ci_mean = (mu_mle - t_critical * se_mean, 
               mu_mle + t_critical * se_mean)
    
    # Confidence interval for σ² using chi-squared distribution
    # (n-1)s²/σ² ~ χ²(n-1)
    chi2_lower = stats.chi2.ppf(0.025, df=n-1)
    chi2_upper = stats.chi2.ppf(0.975, df=n-1)
    ci_var = ((n-1) * var_unbiased / chi2_upper,
              (n-1) * var_unbiased / chi2_lower)
    
    return {
        'n': n,
        'mu_mle': mu_mle,
        'sigma2_mle': var_mle,
        'sigma2_unbiased': var_unbiased,
        'sigma_unbiased': sigma_unbiased,
        'se_mean': se_mean,
        'ci_mean_95': ci_mean,
        'ci_var_95': ci_var
    }
 
# Example: Estimate parameters from data
np.random.seed(42)
true_mu, true_sigma = 50, 10
n = 100
data = np.random.normal(true_mu, true_sigma, n)
 
results = gaussian_mle(data)
print(f"True parameters: μ = {true_mu}, σ = {true_sigma}, σ² = {true_sigma**2}")
print(f"
MLE estimates:")
print(f"  μ̂ = {results['mu_mle']:.4f}")
print(f"  σ̂² (MLE) = {results['sigma2_mle']:.4f}")
print(f"  σ̂² (unbiased) = {results['sigma2_unbiased']:.4f}")
print(f"
95% CI for μ: [{results['ci_mean_95'][0]:.4f}, {results['ci_mean_95'][1]:.4f}]")
print(f"95% CI for σ²: [{results['ci_var_95'][0]:.4f}, {results['ci_var_95'][1]:.4f}]")
 
# Demonstrate bias in variance estimation
print("
--- Bias Demonstration ---")
n_simulations = 10000
mle_vars = []
unbiased_vars = []
 
for _ in range(n_simulations):
    sample = np.random.normal(true_mu, true_sigma, n)
    mle_vars.append(np.var(sample, ddof=0))
    unbiased_vars.append(np.var(sample, ddof=1))
 
print(f"True σ² = {true_sigma**2}")
print(f"E[σ̂²_MLE] ≈ {np.mean(mle_vars):.4f} (biased)")
print(f"E[s²] ≈ {np.mean(unbiased_vars):.4f} (unbiased)")
print(f"Expected MLE bias: {(n-1)/n * true_sigma**2:.4f}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

import numpy as np
from scipy import stats
 
# 1. LINEAR REGRESSION FROM GAUSSIAN LIKELIHOOD
def linear_regression_mle(X, y):
    """
    Derive linear regression from Gaussian noise assumption.
    y = X @ w + ε, where ε ~ N(0, σ²I)
    
    MLE for w: minimize ||y - Xw||² 
    Solution: w = (X'X)^{-1} X'y
    """
    # Add bias term
    X_bias = np.column_stack([np.ones(len(X)), X])
    
    # MLE (equivalent to OLS)
    w_mle = np.linalg.lstsq(X_bias, y, rcond=None)[0]
    
    # Predicted values
    y_pred = X_bias @ w_mle
    
    # MLE for σ²
    residuals = y - y_pred
    sigma2_mle = np.mean(residuals ** 2)
    
    # Log-likelihood
    n = len(y)
    log_likelihood = -n/2 * np.log(2 * np.pi * sigma2_mle) - n/2
    
    return {
        'weights': w_mle,
        'sigma2': sigma2_mle,
        'log_likelihood': log_likelihood,
        'predictions': y_pred
    }
 
# 2. GAUSSIAN NAIVE BAYES
class GaussianNaiveBayes:
    """
    Gaussian Naive Bayes classifier.
    Assumes features are Gaussian given class: P(x_j | y=c) = N(μ_jc, σ²_jc)
    """
    def fit(self, X, y):
        self.classes = np.unique(y)
        self.n_features = X.shape[1]
        
        # Store class priors and parameters
        self.priors = {}
        self.means = {}
        self.variances = {}
        
        for c in self.classes:
            X_c = X[y == c]
            self.priors[c] = len(X_c) / len(X)
            self.means[c] = X_c.mean(axis=0)
            self.variances[c] = X_c.var(axis=0) + 1e-9  # Smoothing
    
    def _log_likelihood(self, x, c):
        """Compute log P(x | y=c) assuming Gaussian features."""
        log_prob = 0
        for j in range(self.n_features):
            log_prob += stats.norm.logpdf(
                x[j], 
                self.means[c][j], 
                np.sqrt(self.variances[c][j])
            )
        return log_prob
    
    def predict_proba(self, X):
        """Compute posterior probabilities P(y=c | x)."""
        log_posteriors = np.zeros((len(X), len(self.classes)))
        
        for i, x in enumerate(X):
            for j, c in enumerate(self.classes):
                log_posteriors[i, j] = (np.log(self.priors[c]) + 
                                        self._log_likelihood(x, c))
        
        # Normalize using log-sum-exp
        log_norm = np.log(np.exp(log_posteriors).sum(axis=1, keepdims=True))
        posteriors = np.exp(log_posteriors - log_norm)
        return posteriors
    
    def predict(self, X):
        return self.classes[self.predict_proba(X).argmax(axis=1)]
 
# Example usage
np.random.seed(42)
X = np.vstack([
    np.random.normal([0, 0], [1, 1], (50, 2)),
    np.random.normal([3, 3], [1, 1], (50, 2))
])
y = np.array([0]*50 + [1]*50)
 
gnb = GaussianNaiveBayes()
gnb.fit(X, y)
accuracy = (gnb.predict(X) == y).mean()
print(f"Gaussian Naive Bayes accuracy: {accuracy:.2%}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
 
def comprehensive_normality_test(data, alpha=0.05):
    """
    Perform comprehensive normality testing with visual and statistical tests.
    """
    n = len(data)
    results = {'n': n, 'alpha': alpha}
    
    # Descriptive statistics
    results['mean'] = np.mean(data)
    results['std'] = np.std(data, ddof=1)
    results['skewness'] = stats.skew(data)
    results['kurtosis'] = stats.kurtosis(data)  # Excess kurtosis (0 for normal)
    
    # Statistical tests
    # 1. Shapiro-Wilk (best for n < 5000)
    if n <= 5000:
        stat, p = stats.shapiro(data)
        results['shapiro_wilk'] = {'statistic': stat, 'p_value': p, 
                                    'normal': p > alpha}
    
    # 2. Kolmogorov-Smirnov (against standard normal after standardization)
    standardized = (data - np.mean(data)) / np.std(data, ddof=1)
    stat, p = stats.kstest(standardized, 'norm')
    results['kolmogorov_smirnov'] = {'statistic': stat, 'p_value': p,
                                      'normal': p > alpha}
    
    # 3. Anderson-Darling
    result = stats.anderson(data, dist='norm')
    # Use 5% significance level (index 2)
    results['anderson_darling'] = {
        'statistic': result.statistic,
        'critical_value_5%': result.critical_values[2],
        'normal': result.statistic < result.critical_values[2]
    }
    
    # 4. D'Agostino-Pearson (requires n >= 20)
    if n >= 20:
        stat, p = stats.normaltest(data)
        results['dagostino_pearson'] = {'statistic': stat, 'p_value': p,
                                         'normal': p > alpha}
    
    return results
 
def print_normality_results(results):
    print(f"
Normality Test Results (n = {results['n']}, α = {results['alpha']})")
    print("=" * 60)
    print(f"Skewness: {results['skewness']:.4f} (0 for normal)")
    print(f"Excess Kurtosis: {results['kurtosis']:.4f} (0 for normal)")
    print("-" * 60)
    
    if 'shapiro_wilk' in results:
        r = results['shapiro_wilk']
        print(f"Shapiro-Wilk: W={r['statistic']:.4f}, p={r['p_value']:.4f} "
              f"→ {'Normal' if r['normal'] else 'NOT Normal'}")
    
    r = results['kolmogorov_smirnov']
    print(f"Kolmogorov-Smirnov: D={r['statistic']:.4f}, p={r['p_value']:.4f} "
          f"→ {'Normal' if r['normal'] else 'NOT Normal'}")
    
    r = results['anderson_darling']
    print(f"Anderson-Darling: A²={r['statistic']:.4f}, crit={r['critical_value_5%']:.4f} "
          f"→ {'Normal' if r['normal'] else 'NOT Normal'}")
    
    if 'dagostino_pearson' in results:
        r = results['dagostino_pearson']
        print(f"D'Agostino-Pearson: K²={r['statistic']:.4f}, p={r['p_value']:.4f} "
              f"→ {'Normal' if r['normal'] else 'NOT Normal'}")
 
# Test with different distributions
np.random.seed(42)
 
print("
" + "="*60)
print("Testing NORMAL data:")
normal_data = np.random.normal(50, 10, 200)
print_normality_results(comprehensive_normality_test(normal_data))
 
print("
" + "="*60)
print("Testing EXPONENTIAL data (should reject normality):")
exp_data = np.random.exponential(10, 200)
print_normality_results(comprehensive_normality_test(exp_data))