Conjugate Priors

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

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
 
def beta_binomial_sensitivity(
    successes: int, 
    trials: int,
    prior_list: list
) -> dict:
    """
    Conduct sensitivity analysis for Beta-Binomial model.
    
    Args:
        successes: Observed successes
        trials: Total trials
        prior_list: List of (name, alpha, beta) tuples
    
    Returns:
        Dictionary with posterior summaries for each prior
    """
    results = {}
    
    for name, alpha_prior, beta_prior in prior_list:
        # Posterior parameters
        alpha_post = alpha_prior + successes
        beta_post = beta_prior + trials - successes
        
        # Posterior summaries
        mean = alpha_post / (alpha_post + beta_post)
        var = (alpha_post * beta_post) / ((alpha_post + beta_post)**2 * (alpha_post + beta_post + 1))
        ci_lower = stats.beta.ppf(0.025, alpha_post, beta_post)
        ci_upper = stats.beta.ppf(0.975, alpha_post, beta_post)
        
        # Probability > 0.1 (example decision threshold)
        prob_above_threshold = 1 - stats.beta.cdf(0.1, alpha_post, beta_post)
        
        results[name] = {
            'prior': (alpha_prior, beta_prior),
            'posterior': (alpha_post, beta_post),
            'mean': mean,
            'std': np.sqrt(var),
            'ci': (ci_lower, ci_upper),
            'P(θ > 0.1)': prob_above_threshold
        }
    
    return results
 
 
def visualize_sensitivity(results: dict, title: str = "Sensitivity Analysis"):
    """
    Create visualization comparing posteriors under different priors.
    """
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    theta = np.linspace(0, 0.3, 500)
    colors = plt.cm.tab10(np.linspace(0, 1, len(results)))
    
    # Left plot: Posterior densities
    ax = axes[0]
    for i, (name, res) in enumerate(results.items()):
        alpha, beta = res['posterior']
        pdf = stats.beta.pdf(theta, alpha, beta)
        ax.plot(theta, pdf, color=colors[i], lw=2, label=name)
        ax.fill_between(theta, pdf, alpha=0.1, color=colors[i])
    
    ax.set_xlabel('θ', fontsize=12)
    ax.set_ylabel('Posterior Density', fontsize=12)
    ax.set_title('Posterior Distributions Under Different Priors', fontsize=12)
    ax.legend(fontsize=10)
    ax.grid(alpha=0.3)
    
    # Right plot: Summary comparison
    ax = axes[1]
    names = list(results.keys())
    means = [res['mean'] for res in results.values()]
    ci_lowers = [res['ci'][0] for res in results.values()]
    ci_uppers = [res['ci'][1] for res in results.values()]
    
    y_pos = np.arange(len(names))
    ax.errorbar(means, y_pos, 
                xerr=[np.array(means) - ci_lowers, np.array(ci_uppers) - np.array(means)],
                fmt='o', capsize=5, markersize=8)
    ax.axvline(0.1, color='red', linestyle='--', label='Threshold (0.1)')
    ax.set_yticks(y_pos)
    ax.set_yticklabels(names)
    ax.set_xlabel('θ', fontsize=12)
    ax.set_title('Posterior Means with 95% CIs', fontsize=12)
    ax.legend()
    ax.grid(alpha=0.3)
    
    plt.tight_layout()
    return fig
 
 
# Example: Testing a new drug for side effect rate
# Observed: 8 side effects in 100 patients
 
successes = 8
trials = 100
 
# Define range of reasonable priors
priors = [
    ("Uniform Beta(1,1)", 1, 1),
    ("Jeffreys Beta(0.5,0.5)", 0.5, 0.5),
    ("Skeptical Beta(1,9)", 1, 9),           # Prior mean 0.1
    ("Optimistic Beta(1,19)", 1, 19),        # Prior mean 0.05
    ("Historical Beta(10,90)", 10, 90),      # Strong prior at 0.1
]
 
results = beta_binomial_sensitivity(successes, trials, priors)
 
print("SENSITIVITY ANALYSIS")
print("=" * 60)
print(f"Data: {successes} side effects in {trials} patients")
print("-" * 60)
 
for name, res in results.items():
    print(f"\n{name}:")
    print(f"  Prior: Beta{res['prior']}")
    print(f"  Posterior mean: {res['mean']:.4f} ± {res['std']:.4f}")
    print(f"  95% CI: [{res['ci'][0]:.4f}, {res['ci'][1]:.4f}]")
    print(f"  P(rate > 10%): {res['P(θ > 0.1)']:.4f}")
 
# Check if conclusions are robust
probs = [res['P(θ > 0.1)'] for res in results.values()]
print("\n" + "=" * 60)
print("ROBUSTNESS CHECK:")
if all(p > 0.95 for p in probs):
    print("✓ All priors conclude P(θ > 0.1) > 95% - ROBUST CONCLUSION")
elif all(p < 0.05 for p in probs):
    print("✓ All priors conclude P(θ > 0.1) < 5% - ROBUST CONCLUSION")
else:
    print("⚠ Conclusions depend on prior choice - SENSITIVE")
    print(f"  Range of P(θ > 0.1): [{min(probs):.4f}, {max(probs):.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

import numpy as np
from scipy.special import gammaln, betaln
 
def beta_binomial_log_marginal(k: int, n: int, alpha: float, beta: float) -> float:
    """
    Log marginal likelihood for Beta-Binomial model.
    
    This is log P(k | n, alpha, beta) = log ∫ P(k | θ, n) P(θ | α, β) dθ
    
    With conjugacy, this has closed form.
    """
    # Log binomial coefficient
    log_binom = gammaln(n + 1) - gammaln(k + 1) - gammaln(n - k + 1)
    
    # Log ratio of Beta functions
    log_beta_ratio = (betaln(alpha + k, beta + n - k) - betaln(alpha, beta))
    
    return log_binom + log_beta_ratio
 
 
def compare_priors_bayes_factor(
    k: int, n: int,
    prior1: tuple, prior1_name: str,
    prior2: tuple, prior2_name: str
) -> dict:
    """
    Compare two prior specifications using Bayes factor.
    """
    log_ml1 = beta_binomial_log_marginal(k, n, prior1[0], prior1[1])
    log_ml2 = beta_binomial_log_marginal(k, n, prior2[0], prior2[1])
    
    log_bf = log_ml1 - log_ml2
    bf = np.exp(log_bf)
    
    # Interpretation
    if bf > 100:
        interpretation = f"Decisive evidence for {prior1_name}"
    elif bf > 10:
        interpretation = f"Strong evidence for {prior1_name}"
    elif bf > 3:
        interpretation = f"Moderate evidence for {prior1_name}"
    elif bf > 1/3:
        interpretation = "Inconclusive"
    elif bf > 1/10:
        interpretation = f"Moderate evidence for {prior2_name}"
    elif bf > 1/100:
        interpretation = f"Strong evidence for {prior2_name}"
    else:
        interpretation = f"Decisive evidence for {prior2_name}"
    
    return {
        'log_marginal_1': log_ml1,
        'log_marginal_2': log_ml2,
        'bayes_factor': bf,
        'log_bayes_factor': log_bf,
        'interpretation': interpretation
    }
 
 
# Example: Does historical data support an informative prior?
 
# Observed: 12 conversions in 150 visitors
k, n = 12, 150
 
# Model 1: Informative prior based on historical data (8% conversion)
prior_informative = (8, 92)  # Beta(8, 92), mean = 0.08
 
# Model 2: Non-informative prior
prior_uniform = (1, 1)  # Beta(1, 1)
 
result = compare_priors_bayes_factor(
    k, n,
    prior_informative, "Informative (historical)",
    prior_uniform, "Uniform"
)
 
print("BAYESIAN MODEL COMPARISON")
print("=" * 60)
print(f"Data: {k} conversions in {n} visitors ({100*k/n:.1f}%)")
print(f"\nModel 1: Beta{prior_informative} (prior mean = {prior_informative[0]/sum(prior_informative):.3f})")
print(f"Model 2: Beta{prior_uniform} (prior mean = {prior_uniform[0]/sum(prior_uniform):.3f})")
print(f"\nLog marginal likelihood (Model 1): {result['log_marginal_1']:.4f}")
print(f"Log marginal likelihood (Model 2): {result['log_marginal_2']:.4f}")
print(f"\nBayes Factor (M1 vs M2): {result['bayes_factor']:.4f}")
print(f"Interpretation: {result['interpretation']}")
 
# The informative prior, if well-calibrated, should have higher marginal likelihood
# because it doesn't waste probability on implausible parameter values

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

import json
from dataclasses import dataclass, asdict
from typing import Optional
import numpy as np
from scipy import stats
 
@dataclass
class BayesianTracker:
    """
    Production-ready Bayesian tracker for binary outcomes.
    
    Features:
    - O(1) update complexity
    - Serializable state (JSON-compatible)
    - Thread-safe design (immutable updates)
    - Full posterior access
    """
    alpha: float
    beta: float
    n_observations: int = 0
    
    @classmethod
    def create(cls, alpha_prior: float = 1.0, beta_prior: float = 1.0):
        """Factory method with default uniform prior."""
        return cls(alpha=alpha_prior, beta=beta_prior, n_observations=0)
    
    def update(self, successes: int, failures: int) -> 'BayesianTracker':
        """
        Return new tracker with updated posterior.
        Original tracker is unchanged (immutable pattern).
        """
        return BayesianTracker(
            alpha=self.alpha + successes,
            beta=self.beta + failures,
            n_observations=self.n_observations + successes + failures
        )
    
    @property
    def mean(self) -> float:
        """Posterior mean."""
        return self.alpha / (self.alpha + self.beta)
    
    @property
    def variance(self) -> float:
        """Posterior variance."""
        a, b = self.alpha, self.beta
        return (a * b) / ((a + b)**2 * (a + b + 1))
    
    def credible_interval(self, confidence: float = 0.95) -> tuple:
        """Equal-tailed credible interval."""
        tail = (1 - confidence) / 2
        return (
            stats.beta.ppf(tail, self.alpha, self.beta),
            stats.beta.ppf(1 - tail, self.alpha, self.beta)
        )
    
    def prob_greater_than(self, threshold: float) -> float:
        """P(θ > threshold | data)."""
        return 1 - stats.beta.cdf(threshold, self.alpha, self.beta)
    
    def to_json(self) -> str:
        """Serialize to JSON for storage/transmission."""
        return json.dumps(asdict(self))
    
    @classmethod
    def from_json(cls, json_str: str) -> 'BayesianTracker':
        """Deserialize from JSON."""
        return cls(**json.loads(json_str))
    
    @classmethod
    def combine(cls, trackers: list) -> 'BayesianTracker':
        """
        Combine multiple trackers (e.g., from parallel workers).
        
        For Beta-Binomial, this is exact: sum the count contributions.
        """
        if not trackers:
            return cls.create()
        
        # Subtract 1 from each alpha/beta to get counts, sum, add back
        # Actually: just use the formula
        # If trackers share the same prior, combine as:
        # α_combined = α_prior + Σ(α_i - α_prior) = α_prior + Σ successes_i
        
        # For correctness, assume all trackers started with same prior
        prior_alpha = trackers[0].alpha - trackers[0].n_observations * trackers[0].mean
        prior_beta = trackers[0].beta - trackers[0].n_observations * (1 - trackers[0].mean)
        
        # Simple approach: sum increments over prior
        total_alpha = sum(t.alpha for t in trackers) - (len(trackers) - 1) * prior_alpha
        total_beta = sum(t.beta for t in trackers) - (len(trackers) - 1) * prior_beta
        total_n = sum(t.n_observations for t in trackers)
        
        return cls(alpha=total_alpha, beta=total_beta, n_observations=total_n)
 
 
# Example: Distributed A/B Test Processing
# Simulate parallel workers each processing a shard of data
 
np.random.seed(42)
true_rate_A = 0.10
true_rate_B = 0.12
 
# Create trackers with shared prior
prior_A = BayesianTracker.create(alpha_prior=1, beta_prior=1)
prior_B = BayesianTracker.create(alpha_prior=1, beta_prior=1)
 
# Simulate 3 workers each processing 1000 users per variant
n_workers = 3
n_per_worker = 1000
 
worker_trackers_A = []
worker_trackers_B = []
 
for w in range(n_workers):
    # Worker w's data
    successes_A = np.random.binomial(n_per_worker, true_rate_A)
    successes_B = np.random.binomial(n_per_worker, true_rate_B)
    
    tracker_A = prior_A.update(successes_A, n_per_worker - successes_A)
    tracker_B = prior_B.update(successes_B, n_per_worker - successes_B)
    
    worker_trackers_A.append(tracker_A)
    worker_trackers_B.append(tracker_B)
    
    print(f"Worker {w}: A={successes_A}/{n_per_worker}, B={successes_B}/{n_per_worker}")
 
# Combine results from all workers
combined_A = BayesianTracker.combine(worker_trackers_A)
combined_B = BayesianTracker.combine(worker_trackers_B)
 
print(f"\nCombined Results:")
print(f"Variant A: mean={combined_A.mean:.4f}, 95% CI={combined_A.credible_interval()}")
print(f"Variant B: mean={combined_B.mean:.4f}, 95% CI={combined_B.credible_interval()}")
 
# Serialize for storage
state_A = combined_A.to_json()
print(f"\nSerialized state: {state_A}")
 
# Restore later
restored_A = BayesianTracker.from_json(state_A)
print(f"Restored mean: {restored_A.mean:.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

import numpy as np
from scipy import stats
 
def simulation_based_calibration_check(
    alpha_prior: float, 
    beta_prior: float,
    n_trials: int,
    n_simulations: int = 10000
) -> dict:
    """
    Check calibration of Beta-Binomial model via simulation.
    
    For a well-calibrated model, the true parameter should fall
    in the X% credible interval about X% of the time.
    """
    results = {
        'in_50_ci': 0,
        'in_80_ci': 0,
        'in_95_ci': 0,
        'in_99_ci': 0,
        'rank_statistics': []  # For detailed calibration check
    }
    
    for _ in range(n_simulations):
        # 1. Sample true theta from prior
        true_theta = np.random.beta(alpha_prior, beta_prior)
        
        # 2. Generate data from true theta
        k = np.random.binomial(n_trials, true_theta)
        
        # 3. Compute posterior
        alpha_post = alpha_prior + k
        beta_post = beta_prior + n_trials - k
        
        # 4. Check if true_theta falls in various CIs
        for ci_level, key in [(0.50, 'in_50_ci'), (0.80, 'in_80_ci'),
                               (0.95, 'in_95_ci'), (0.99, 'in_99_ci')]:
            tail = (1 - ci_level) / 2
            lower = stats.beta.ppf(tail, alpha_post, beta_post)
            upper = stats.beta.ppf(1 - tail, alpha_post, beta_post)
            if lower <= true_theta <= upper:
                results[key] += 1
        
        # 5. Rank statistic: what quantile is true_theta in posterior?
        rank = stats.beta.cdf(true_theta, alpha_post, beta_post)
        results['rank_statistics'].append(rank)
    
    # Compute coverage rates
    for key in ['in_50_ci', 'in_80_ci', 'in_95_ci', 'in_99_ci']:
        results[key] = results[key] / n_simulations
    
    return results
 
 
# Check calibration with different sample sizes
print("CALIBRATION CHECK: Beta-Binomial Model")
print("=" * 60)
 
for n_trials in [10, 50, 200]:
    results = simulation_based_calibration_check(
        alpha_prior=1, beta_prior=1, n_trials=n_trials, n_simulations=10000
    )
    
    print(f"\nn_trials = {n_trials}")
    print(f"  50% CI coverage: {results['in_50_ci']:.3f} (target: 0.500)")
    print(f"  80% CI coverage: {results['in_80_ci']:.3f} (target: 0.800)")
    print(f"  95% CI coverage: {results['in_95_ci']:.3f} (target: 0.950)")
    print(f"  99% CI coverage: {results['in_99_ci']:.3f} (target: 0.990)")
    
    # Rank statistics should be uniform if calibrated
    ranks = np.array(results['rank_statistics'])
    ks_stat, ks_pvalue = stats.kstest(ranks, 'uniform')
    print(f"  Rank uniformity (KS test p-value): {ks_pvalue:.4f}")
    if ks_pvalue > 0.05:
        print("  ✓ Well-calibrated (ranks are uniform)")
    else:
        print("  ⚠ Potential calibration issue")
 
# The Beta-Binomial model should be perfectly calibrated since we're 
# simulating from the exact same model we're fitting