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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

import numpy as np
from scipy import stats
import warnings
 
def demonstrate_correlation_blindness():
    """
    Show how mean-field completely misses posterior correlations.
    
    Example: Bayesian linear regression with correlated predictors.
    """
    
    print("Correlation Blindness in Mean-Field VI")
    print("=" * 60)
    
    # Generate data with correlated predictors
    np.random.seed(42)
    n = 100
    
    # Correlated predictors: x1 and x2 have correlation 0.9
    rho = 0.9
    cov_x = np.array([[1, rho], [rho, 1]])
    X = np.random.multivariate_normal([0, 0], cov_x, n)
    
    # True coefficients
    true_beta = np.array([1.0, 1.0])
    sigma = 1.0
    
    # Generate response
    y = X @ true_beta + np.random.randn(n) * sigma
    
    print(f"
Data: n={n}, p=2")
    print(f"Predictor correlation: {rho}")
    print(f"True coefficients: {true_beta}")
    
    # =====================================================
    # Compute true posterior (analytically for linear regression)
    # =====================================================
    
    # Prior: β ~ N(0, τ²I) with τ² = 10
    tau_sq = 10.0
    prior_precision = np.eye(2) / tau_sq
    
    # Posterior: β | y ~ N(μ_post, Σ_post)
    # Σ_post^{-1} = X'X / σ² + prior_precision
    # μ_post = Σ_post @ X'y / σ²
    
    XtX = X.T @ X
    Xty = X.T @ y
    
    post_precision = XtX / sigma**2 + prior_precision
    post_cov = np.linalg.inv(post_precision)
    post_mean = post_cov @ Xty / sigma**2
    
    post_corr = post_cov[0, 1] / np.sqrt(post_cov[0, 0] * post_cov[1, 1])
    
    print(f"
True Posterior:")
    print(f"  Means: [{post_mean[0]:.3f}, {post_mean[1]:.3f}]")
    print(f"  Variances: [{post_cov[0,0]:.4f}, {post_cov[1,1]:.4f}]")
    print(f"  Correlation: {post_corr:.3f}")
    print(f"  Covariance: {post_cov[0,1]:.4f}")
    
    # =====================================================
    # Mean-field approximation
    # =====================================================
    
    # Mean-field: q(β) = q(β1) q(β2)
    # Each q(βj) is Gaussian, but independent
    
    # The mean-field optimal marginals have the same means
    # but diagonal covariance (ignoring off-diagonal)
    
    # For Gaussian, mean-field marginals match true marginal means
    mf_means = post_mean.copy()
    
    # But variances come from diagonal of precision matrix
    # (not diagonal of covariance!)
    # For mean-field: var_j = 1 / precision_jj
    mf_vars = 1.0 / np.diag(post_precision)
    
    print(f"
Mean-Field Approximation:")
    print(f"  Means: [{mf_means[0]:.3f}, {mf_means[1]:.3f}]")
    print(f"  Variances: [{mf_vars[0]:.4f}, {mf_vars[1]:.4f}]")
    print(f"  Correlation: 0.000 (by construction)")
    print(f"  Covariance: 0.000 (by construction)")
    
    # =====================================================
    # Compare consequences
    # =====================================================
    
    print(f"
" + "=" * 60)
    print("Consequences of Correlation Blindness:")
    print("=" * 60)
    
    # 1. Variance comparison
    print(f"
1. Marginal Variance Comparison:")
    print(f"   True variances: [{post_cov[0,0]:.4f}, {post_cov[1,1]:.4f}]")
    print(f"   MF variances:   [{mf_vars[0]:.4f}, {mf_vars[1]:.4f}]")
    print(f"   MF underestimates variance by factor: "
          f"[{post_cov[0,0]/mf_vars[0]:.2f}x, {post_cov[1,1]/mf_vars[1]:.2f}x]")
    
    # 2. Sum of coefficients
    print(f"
2. Variance of β₁ + β₂:")
    true_var_sum = post_cov[0,0] + post_cov[1,1] + 2*post_cov[0,1]
    mf_var_sum = mf_vars[0] + mf_vars[1]  # No covariance term!
    print(f"   True: Var(β₁ + β₂) = {true_var_sum:.4f}")
    print(f"   MF:   Var(β₁ + β₂) = {mf_var_sum:.4f}")
    print(f"   Error: {abs(true_var_sum - mf_var_sum)/true_var_sum * 100:.1f}%")
    
    # 3. Credible regions
    print(f"
3. Joint 95% Credible Region:")
    print(f"   True: ellipse (captures {post_corr:.2f} correlation)")
    print(f"   MF:   axis-aligned rectangle (wrong shape entirely)")
    
    # 4. Predictive at new point
    x_new = np.array([1.0, 1.0])
    
    true_pred_var = x_new @ post_cov @ x_new + sigma**2
    mf_pred_var = x_new @ np.diag(mf_vars) @ x_new + sigma**2
    
    print(f"
4. Predictive Variance at x=[1, 1]:")
    print(f"   True: {true_pred_var:.4f}")
    print(f"   MF:   {mf_pred_var:.4f}")
    print(f"   MF overestimates by {(mf_pred_var/true_pred_var - 1)*100:.1f}%")
    
    return {
        'true_cov': post_cov,
        'mf_vars': mf_vars,
        'post_corr': post_corr
    }
 
 
if __name__ == "__main__":
    demonstrate_correlation_blindness()

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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

import numpy as np
from scipy import stats
 
def analyze_variance_underestimation():
    """
    Quantify how mean-field underestimates variance.
    """
    
    print("Variance Underestimation in Mean-Field VI")
    print("=" * 60)
    
    # True bivariate Gaussian posterior
    # Strong positive correlation
    rho = 0.8
    true_cov = np.array([[1.0, rho], [rho, 1.0]])
    
    # Mean-field: diagonal covariance
    # Best MF approximation to match marginals
    mf_cov = np.diag([1.0, 1.0])
    
    print(f"
True posterior: bivariate Gaussian, ρ = {rho}")
    print(f"Mean-field: product of independent Gaussians")
    
    # Compare variances of various quantities
    quantities = [
        ("z₁", np.array([1, 0])),
        ("z₂", np.array([0, 1])),
        ("z₁ + z₂", np.array([1, 1])),
        ("z₁ - z₂", np.array([1, -1])),
        ("2z₁ + z₂", np.array([2, 1])),
    ]
    
    print(f"
{'Quantity':<12} {'True Var':<12} {'MF Var':<12} {'Ratio':<10} {'Error'}")
    print("-" * 60)
    
    for name, w in quantities:
        true_var = w @ true_cov @ w
        mf_var = w @ mf_cov @ w
        ratio = true_var / mf_var
        error = (mf_var - true_var) / true_var * 100
        
        print(f"{name:<12} {true_var:<12.4f} {mf_var:<12.4f} {ratio:<10.2f} {error:+.1f}%")
    
    print()
    print("Key observations:")
    print("- For z₁ + z₂: True var includes positive cov term, MF misses it")
    print("- For z₁ - z₂: True var is reduced by positive cov, MF overestimates")
    print("- Mean-field errors compound for linear combinations")
    
    # Coverage probability experiment
    print(f"
" + "=" * 60)
    print("Credible Interval Coverage Analysis")
    print("=" * 60)
    
    n_samples = 10000
    
    # Generate samples from true posterior
    true_samples = np.random.multivariate_normal([0, 0], true_cov, n_samples)
    
    # Mean-field "posterior" 
    mf_sample_z1 = np.random.randn(n_samples)
    mf_sample_z2 = np.random.randn(n_samples)
    
    # For z1 + z2, compute 95% intervals
    true_sum = true_samples[:, 0] + true_samples[:, 1]
    mf_sum = mf_sample_z1 + mf_sample_z2
    
    # True 95% interval
    true_lower = np.percentile(true_sum, 2.5)
    true_upper = np.percentile(true_sum, 97.5)
    
    # MF 95% interval
    mf_lower = np.percentile(mf_sum, 2.5)
    mf_upper = np.percentile(mf_sum, 97.5)
    
    print(f"
95% Credible Intervals for z₁ + z₂:")
    print(f"  True: [{true_lower:.3f}, {true_upper:.3f}] (width: {true_upper-true_lower:.3f})")
    print(f"  MF:   [{mf_lower:.3f}, {mf_upper:.3f}] (width: {mf_upper-mf_lower:.3f})")
    
    # What coverage does MF interval achieve under true posterior?
    mf_coverage = np.mean((true_sum >= mf_lower) & (true_sum <= mf_upper))
    print(f"
MF '95%' interval actual coverage: {mf_coverage*100:.1f}%")
    print(f"MF interval is too narrow!")
 
 
def scaling_with_dimensions():
    """
    Show how variance underestimation scales with dimensionality.
    """
    
    print("
" + "=" * 60)
    print("Scaling of Variance Underestimation with Dimension")
    print("=" * 60)
    print()
    
    rho = 0.3  # Moderate correlation
    
    print(f"Equal correlation ρ = {rho} between all pairs")
    print(f"Quantity: sum of all variables")
    print()
    print(f"{'Dimension':<12} {'True Var(sum)':<15} {'MF Var(sum)':<15} {'Underest. Factor'}")
    print("-" * 60)
    
    for d in [2, 5, 10, 20, 50]:
        # True covariance: compound symmetry
        true_cov = rho * np.ones((d, d)) + (1 - rho) * np.eye(d)
        
        # MF: diagonal
        mf_cov = np.eye(d)
        
        # Variance of sum = 1ᵀ Σ 1
        ones = np.ones(d)
        true_var_sum = ones @ true_cov @ ones
        mf_var_sum = ones @ mf_cov @ ones
        
        factor = true_var_sum / mf_var_sum
        
        print(f"{d:<12} {true_var_sum:<15.2f} {mf_var_sum:<15.2f} {factor:.2f}x")
    
    print()
    print("As dimensionality increases, the underestimation of")
    print("variance for sums/averages gets WORSE, not better!")
 
 
if __name__ == "__main__":
    analyze_variance_underestimation()
    scaling_with_dimensions()

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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

import numpy as np
from scipy import stats
 
def demonstrate_mode_seeking():
    """
    Show how mean-field focuses on single modes.
    """
    
    print("Mode-Seeking Behavior in Mean-Field VI")
    print("=" * 60)
    
    # True posterior: mixture of two Gaussians
    # Represents a bimodal posterior over a single parameter
    
    mu1, sigma1 = -2.0, 0.5
    mu2, sigma2 = 2.0, 0.5
    weight1 = 0.6  # Slightly more mass in mode 1
    
    print(f"
True posterior: mixture of Gaussians")
    print(f"  Mode 1: N({mu1}, {sigma1}²) with weight {weight1}")
    print(f"  Mode 2: N({mu2}, {sigma2}²) with weight {1-weight1}")
    
    # True posterior statistics
    true_mean = weight1 * mu1 + (1 - weight1) * mu2
    true_var = (weight1 * (sigma1**2 + mu1**2) + 
                (1 - weight1) * (sigma2**2 + mu2**2) - 
                true_mean**2)
    
    print(f"
True posterior statistics:")
    print(f"  Mean: {true_mean:.3f}")
    print(f"  Variance: {true_var:.3f}")
    print(f"  Std Dev: {np.sqrt(true_var):.3f}")
    
    # Mean-field: single Gaussian
    # KL minimization will focus on ONE mode
    
    print(f"
Mean-field approximation (Gaussian):")
    
    # Mean-field focused on mode 1
    mf_mode1_mean = mu1
    mf_mode1_var = sigma1**2
    
    # Mean-field focused on mode 2
    mf_mode2_mean = mu2
    mf_mode2_var = sigma2**2
    
    print(f"
  If initialized near mode 1:")
    print(f"    Mean: {mf_mode1_mean:.3f}")
    print(f"    Variance: {mf_mode1_var:.3f}")
    
    print(f"
  If initialized near mode 2:")
    print(f"    Mean: {mf_mode2_mean:.3f}")
    print(f"    Variance: {mf_mode2_var:.3f}")
    
    # What if we naively averaged across modes?
    print(f"
  Naive Gaussian matching true mean/var:")
    print(f"    Mean: {true_mean:.3f}")
    print(f"    Variance: {true_var:.3f}")
    print(f"    (This is NOT what mean-field produces!)")
    
    # Coverage analysis
    print(f"
" + "=" * 60)
    print("Coverage Analysis")
    print("=" * 60)
    
    # Generate samples from true bimodal posterior
    n = 10000
    assignments = np.random.binomial(1, weight1, n)
    true_samples = np.where(
        assignments,
        np.random.normal(mu1, sigma1, n),
        np.random.normal(mu2, sigma2, n)
    )
    
    # Mean-field 95% CI (using mode 1)
    mf1_ci_low = mu1 - 1.96 * sigma1
    mf1_ci_high = mu1 + 1.96 * sigma1
    
    coverage_mf1 = np.mean((true_samples >= mf1_ci_low) & 
                           (true_samples <= mf1_ci_high))
    
    print(f"
Mean-field (mode 1) '95%' CI: [{mf1_ci_low:.2f}, {mf1_ci_high:.2f}]")
    print(f"  Actual coverage: {coverage_mf1*100:.1f}%")
    print(f"  COMPLETELY misses mode 2!")
    
    # True 95% CI (empirical)
    true_ci_low = np.percentile(true_samples, 2.5)
    true_ci_high = np.percentile(true_samples, 97.5)
    
    print(f"
True 95% CI: [{true_ci_low:.2f}, {true_ci_high:.2f}]")
    print(f"  (Much wider, covers both modes)")
 
 
def mixture_label_switching():
    """
    Illustrate label switching issue in mixture models.
    """
    
    print("
" + "=" * 60)
    print("Label Switching in Mixture Models")
    print("=" * 60)
    
    print("""
    Consider a 2-component Gaussian mixture:
    
      p(data | μ₁, μ₂, π) = Σᵢ [π N(xᵢ | μ₁, σ²) + (1-π) N(xᵢ | μ₂, σ²)]
    
    The likelihood is SYMMETRIC in (μ₁, μ₂):
      Swapping μ₁ ↔ μ₂ and π ↔ (1-π) gives identical likelihood.
    
    True posterior has (at least) 2 equivalent modes:
      Mode A: μ₁ ≈ 0, μ₂ ≈ 5 
      Mode B: μ₁ ≈ 5, μ₂ ≈ 0  (swapped)
    
    Mean-field will pick ONE mode based on initialization.
    It cannot represent both simultaneously.
    
    Consequences:
    • Posterior mean of μ₁ is NOT meaningful
    • Must apply 'label processing' to interpret results
    • Uncertainty about cluster means is underestimated
    """)
 
 
if __name__ == "__main__":
    demonstrate_mode_seeking()
    mixture_label_switching()

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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

import numpy as np
 
def funnel_posterior_failure():
    """
    Demonstrate mean-field failure on funnel posteriors.
    
    The Neal's funnel is a canonical example from hierarchical models.
    """
    
    print("Funnel Posterior Failure Case")
    print("=" * 60)
    
    print("""
    Neal's Funnel:
    
      v ~ N(0, 3²)              # Wide prior on log-variance
      x | v ~ N(0, exp(v))      # Conditional Gaussian
    
    The marginal posterior p(x, v | data) has a funnel shape:
    - When v is large: x has high variance (wide)
    - When v is small: x has low variance (narrow)
    
    Mean-field q(x)q(v) uses FIXED variance for x
    - Cannot adapt to v-dependent width
    - Either too wide at narrow end OR too narrow at wide end
    """)
    
    # Simulate the issue
    n_samples = 10000
    
    # True funnel samples
    v_samples = np.random.randn(n_samples) * 3
    x_samples = np.random.randn(n_samples) * np.exp(v_samples / 2)
    
    # Analyze x variance conditional on v
    v_bins = [(-np.inf, -2), (-2, 0), (0, 2), (2, np.inf)]
    
    print("
True x variance by v region:")
    for low, high in v_bins:
        mask = (v_samples >= low) & (v_samples < high)
        if np.sum(mask) > 100:
            x_var = np.var(x_samples[mask])
            print(f"  v ∈ [{low:+.0f}, {high:+.0f}): Var(x) = {x_var:.2f}")
    
    # Mean-field uses single variance for all v
    mf_x_var = np.var(x_samples)
    print(f"
Mean-field uses single Var(x) = {mf_x_var:.2f}")
    print("This is WAY too large for small v, too small for large v!")
 
 
def constraint_failure():
    """
    Demonstrate failure when posterior has constraints.
    """
    
    print("
" + "=" * 60)
    print("Constraint Satisfaction Failure")
    print("=" * 60)
    
    print("""
    Scenario: Simplex constraint
    
    True posterior: (z₁, z₂, z₃) must sum to 1
      Lives on a 2D triangle in 3D space
    
    Mean-field: q(z₁)q(z₂)q(z₃)
      Each marginal is independent
      Product is a 3D region
    
    Problem:
      Mean-field places probability OUTSIDE the valid simplex!
      Probability leaks to invalid regions.
    """)
    
    # Simulate
    n_samples = 10000
    
    # True samples: on simplex
    true_samples = np.random.dirichlet([2, 2, 2], n_samples)
    
    # Mean-field: match marginals independently
    mf_means = np.mean(true_samples, axis=0)
    mf_vars = np.var(true_samples, axis=0)
    
    print(f"True marginal means: {mf_means.round(3)}")
    print(f"True marginal vars: {mf_vars.round(4)}")
    
    # Sample from mean-field (independent marginals)
    # Using truncated normal as crude approximation
    mf_samples = np.column_stack([
        np.random.normal(mf_means[i], np.sqrt(mf_vars[i]), n_samples)
        for i in range(3)
    ])
    
    # Check how many violate simplex
    sums = mf_samples.sum(axis=1)
    valid = (mf_samples >= 0).all(axis=1) & np.isclose(sums, 1, atol=0.3)
    
    print(f"
Mean-field samples summing to ~1 (±0.3): {np.mean(valid)*100:.1f}%")
    print(f"Mean-field samples with all positive: {np.mean((mf_samples >= 0).all(axis=1))*100:.1f}%")
    print("
Mean-field puts substantial probability in invalid regions!")
 
 
def summarize_failures():
    """Summary of when to avoid mean-field."""
    
    print("
" + "=" * 60)
    print("When to Avoid Mean-Field VI")
    print("=" * 60)
    
    print("""
    Avoid mean-field when:
    
    1. Strong posterior correlations matter for your analysis
       → Use structured VI or full-covariance VI
    
    2. Accurate uncertainty quantification is critical
       → Use MCMC or importance-weighted VI
    
    3. Posterior is known to be multi-modal
       → Use mixture VI or tempering methods
    
    4. Parameters have constraints (simplex, positive definite, etc.)
       → Use constrained parameterizations
    
    5. Hierarchical models with funnel geometries
       → Use non-centered parameterizations
    
    Use mean-field when:
    
    1. Point estimates are sufficient
    2. Problem scale requires fast inference
    3. Posterior correlations are expected to be weak
    4. Initial exploration before refined analysis
    """)
 
 
if __name__ == "__main__":
    funnel_posterior_failure()
    constraint_failure()
    summarize_failures()