Bayesian Mixture Priors & System Reliability: ISYE 6420 Tutorial

Understanding Bayesian Inference with Mixture Priors

In modern data science, Bayesian methods are widely used to incorporate prior knowledge into statistical models. A common scenario is when the prior distribution itself is a mixture of simpler distributions, reflecting multiple possible states of nature. For example, in A/B testing for a new AI-driven app feature, you might believe the effect size is either small negative or small positive, leading to a mixture normal prior. This tutorial walks through the key steps to compute posterior distributions given a mixture prior, using concepts from ISYE 6420 spring 2025 homework 2.

The Setup: Normal Likelihood with Mixture Normal Prior

Suppose data yi are i.i.d. from N(μ, 1) for i=1,...,n. The prior on μ is a 50-50 mixture of two normal distributions: μ ∼ 0.5 N(−1,1) + 0.5 N(1,1). This prior might represent a belief that the true mean is either near -1 or near 1, with equal probability. After observing a single data point y=1, we want the posterior distribution p(μ|y).

Step 1: Applying Bayes' Theorem

Bayes' theorem gives p(μ|y) ∝ p(y|μ) p(μ). The likelihood p(y|μ) is φ(y; μ, 1). The prior is a sum: 0.5 φ(μ; -1,1) + 0.5 φ(μ; 1,1). Thus, the product becomes a sum of two terms: 0.5 φ(y; μ,1) φ(μ; -1,1) + 0.5 φ(y; μ,1) φ(μ; 1,1).

To simplify each term, we use the identity:

φ(x; μ₁, σ₁²) φ(x; μ₂, σ₂²) = φ(x; (μ₁/σ₁² + μ₂/σ₂²)/(1/σ₁²+1/σ₂²), 1/(1/σ₁²+1/σ₂²)) × φ(μ₁-μ₂; 0, σ₁²+σ₂²)

This identity is a powerful tool in Bayesian conjugate analysis. Applying it to the first term with x=μ, μ₁=y, σ₁²=1, μ₂=-1, σ₂²=1 gives:

φ(μ; y,1) φ(μ; -1,1) = φ(μ; (y/1 + (-1)/1)/(1/1+1/1), 1/(1/1+1/1)) × φ(y-(-1); 0, 1+1) = φ(μ; (y-1)/2, 1/2) × φ(y+1; 0, 2)

Similarly, the second term becomes φ(μ; (y+1)/2, 1/2) × φ(y-1; 0, 2). With y=1, the first term's factor φ(2;0,2) and second term's factor φ(0;0,2).

Step 2: Identifying the Posterior Mixture

The posterior is proportional to a weighted sum of two normal densities. After normalizing, the posterior is a mixture of two normals with means (y-1)/2 = 0 and (y+1)/2 = 1, both with variance 1/2. The weights are proportional to the normalizing constants: w₁ ∝ 0.5 φ(2;0,2) and w₂ ∝ 0.5 φ(0;0,2). Since φ(0;0,2) = 1/(√(4π)) and φ(2;0,2) = 1/(√(4π)) exp(-1), the weights are w₁ = exp(-1)/(1+exp(-1)) ≈ 0.2689 and w₂ = 1/(1+exp(-1)) ≈ 0.7311. Thus, the posterior is:

p(μ|y=1) = 0.2689 N(0, 0.5) + 0.7311 N(1, 0.5)

This result shows how the observed data shifts the prior mixture toward the component that better fits the data.

System Reliability: k-out-of-n Systems with Weibull Components

Reliability engineering often deals with systems that function if at least k out of n components work. Such systems are common in cloud computing (e.g., at least 7 out of 10 servers must be operational) or in autonomous vehicle sensor arrays. Here, we consider a k-out-of-n system with n=10, k=7, and component lifetimes following a Weibull distribution with shape r=1.3 and scale λ=1/20 per month.

Probability the System is Operational at Time t

The lifetime T of a component satisfies P(T ≥ t) = exp(-λ t^r). For t in months, the probability a component works at time t is p(t) = exp(-t^r/20). The number of working components X ~ Binomial(n, p(t)). The system is operational if X ≥ k. Thus,

P(system operational at t) = P(X ≥ k) = 1 - F(k-1; n, p(t))

where F is the binomial CDF. For t=6 months, p(6) = exp(-6^1.3/20). Compute 6^1.3 ≈ 10.27, so p(6)=exp(-10.27/20)=exp(-0.5135)≈0.598. Then P(X≥7) = 1 - binocdf(6,10,0.598). Using Python: 1 - scipy.stats.binom.cdf(6,10,0.598) gives about 0.382.

Conditional Probability Given System Operational at t=6

Given the system was operational at t=6, what is the probability exactly 7 components were working? This is P(X=7 | X≥7). Using Bayes:

P(X=7 | X≥7) = P(X=7) / P(X≥7) = [dbinom(7,10,0.598)] / 0.382

dbinom(7,10,0.598) = C(10,7) * 0.598^7 * (0.402)^3 ≈ 120 * 0.032 * 0.065 ≈ 0.250. So the conditional probability ≈ 0.250/0.382 ≈ 0.654.

Future Reliability Given Past Operation

If the system was operational at t=6, what is the probability it will still be operational at t=9? We need P(system operational at t=9 | X≥7 at t=6). This is more complex because component states are dependent over time. However, due to the memoryless property of the exponential distribution? No, Weibull is not memoryless. We must use the joint distribution. A simpler approach: given that a component was working at t=6, its remaining lifetime distribution is Weibull with a conditional survival function. For a component that survived to t=6, the probability it works at t=9 is P(T≥9 | T≥6) = exp(- (9^1.3 - 6^1.3)/20). Compute 9^1.3 ≈ 18.72, so difference = 18.72-10.27=8.45, then exp(-8.45/20)=exp(-0.4225)=0.655. So the conditional probability a component works at 9 given it worked at 6 is 0.655. Now, given at least 7 worked at 6, the number working at 9 is not simply binomial because the survival probabilities are independent across components? Actually, given each component's status at 6, the future is independent. But we only know that at least 7 worked; we don't know which ones. The exact calculation requires summing over all possible numbers at 6. However, a common approximation is to assume the system is as good as new? Not valid here. For this homework, you might compute by conditioning on the exact number at 6: P(system at 9 | X≥7 at 6) = Σ_{x=7}^{10} P(X_9≥7 | X_6=x) P(X_6=x)/P(X_6≥7). Given X_6=x, the number working at 9 is Binomial(x, p_cond) where p_cond=0.655. So P(X_9≥7 | X_6=x) = 1 - binocdf(6, x, 0.655). Then compute weighted average. This is a typical Bayesian reliability analysis.

Bayesian Inference for a Uniform Distribution with Pareto Prior

Consider i.i.d. data U_i ~ Uniform(0,θ). The maximum M = max(U_i) is sufficient. The likelihood is p(M|θ) = n M^{n-1}/θ^n for θ ≥ M. Using a noninformative Jeffreys prior π(θ)=1/θ, the posterior is Pareto: p(θ|M) = n M^n / θ^{n+1} for θ>M. This is Pa(c=M, α=n). For a dataset with n=351 and M=0.7 (the true θ), the posterior is Pa(0.7, 351).

Point Estimate and Credible Set

A common Bayesian estimate is the posterior mean: E[θ|M] = α c/(α-1) = 351*0.7/350 = 0.702. The 95% equi-tailed credible set is (c, θ_upper) where the lower bound is c (since Pareto is decreasing) and upper bound such that P(θ ≤ θ_upper)=0.95. Solve 1 - (c/θ_upper)^α = 0.95 → (c/θ_upper)^α = 0.05 → θ_upper = c / (0.05)^{1/α} = 0.7 / (0.05)^{1/351}. Since 0.05^{1/351} ≈ exp(ln(0.05)/351) ≈ exp(-2.9957/351) ≈ exp(-0.00853) ≈ 0.9915, so θ_upper ≈ 0.7/0.9915 ≈ 0.706. Thus the 95% credible set is (0.7, 0.706). The true value 0.7 is at the boundary; technically it is included because θ>c, so 0.7 is not in the open interval, but often we consider [c, θ_upper]. In practice, 0.7 is the minimum possible, so the credible set includes it.

Sensitivity to Prior

If we replace Jeffreys prior with a Pareto prior Pa(c0, α0) where c0 < M and α0 small, the posterior becomes Pa(max(c0, M), α0+n). For example, with c0=0.1, α0=0.1, posterior is Pa(0.7, 351.1). The mean becomes 351.1*0.7/350.1 ≈ 0.702, similar. As α0 increases, the prior becomes more informative, pulling the posterior mean toward c0. If c0 is far from M, the posterior may shift. This experiment highlights the robustness of Bayesian inference with large sample sizes.

Visualizing the Posterior

A plot of the posterior PDF for the Pareto case shows a steep decline from θ=0.7. The credible interval marks the region from 0.7 to 0.706, and the posterior mean at 0.702. Such plots are essential for communicating uncertainty.

Conclusion

This tutorial covered three key Bayesian topics: mixture priors, system reliability, and Pareto posteriors. These concepts are directly applicable to modern AI systems (e.g., reliability of neural network ensembles) and financial modeling (e.g., mixture models for asset returns). By mastering these techniques, you'll be well-equipped for advanced Bayesian coursework and real-world data analysis.