Computation and Software

Bayesian data analysis depends fundamentally on computational tools for fitting models, generating simulations, and conducting inference.
R is mainly used for data handling, visualization, regression modeling, optimization, and programming, while Stan is used for Bayesian inference and posterior sampling.
Even when Stan is employed, R remains central for data preprocessing, plotting, and post-analysis model checking.

The general computational philosophy is to start simple and increase model complexity gradually.
Instead of running a large, complex model overnight, one should fit smaller models quickly, use earlier inferences as starting points, and visually assess fit before proceeding further.

---

Key Computational Tasks in Bayesian Analysis

1. Vector and Matrix Manipulations
   • Core to expressing statistical models in algebraic form.
   • Used in regression, linear algebra operations, and computing covariance matrices.
   • Examples include computing $X^T X$, $(X^T X)^{-1} X^T y$, or variance matrices for regression parameters.
2. Computing Probability Densities
   • Evaluating prior $p(\theta)$ and likelihood $p(y|\theta)$ functions.
   • Requires understanding standard distributions (Normal, Gamma, Beta, Bernoulli, etc.).
   • In R, these are easily computed via built-in density functions like dnorm, dbeta, etc.
3. Drawing Random Samples (Simulation)
   • Generating random draws from probability distributions to approximate integrals and posterior distributions.
   • Core to Monte Carlo simulation methods used in Bayesian inference.
4. Structured Programming
   • Includes writing loops, defining user functions, and modular coding for reproducibility and scalability.
   • Used for automating iterative model fitting and diagnostic procedures.
5. Regression Estimation and Variance Calculation
   • Computing estimates and uncertainty of regression coefficients using analytical or simulated methods.
   • Example: $\hat{\beta} = (X^T X)^{-1} X^T y$, with covariance $\sigma^2 (X^T X)^{-1}$.
6. Graphics and Visualization
   • Scatterplots, histograms, overlayed regression lines, and multi-panel graphs are essential for exploratory analysis and posterior checking.

---

Computational Strategy

• Incremental model building: Begin with simple models, evaluate their fit, and progressively increase complexity.
• Iterative validation: After fitting each model, perform posterior predictive checks and visualization before moving on.
• Avoid blind computation: Rather than maximizing precision in one large model, iteratively refine and diagnose.
• Simulation as a diagnostic: Outlier simulated values can signal model or parameterization problems.

---

Summarizing Inferences by Simulation

Simulation is the central mechanism of Bayesian inference.
Even when analytic integration of the posterior is impossible, simulations can approximate it effectively.

A fundamental idea is the duality between a probability density function and a histogram:
given enough simulated samples, the histogram approximates the underlying density $p(\theta)$.

Practical Example

To estimate the 95th percentile of $p(\theta)$:

1. Draw $S$ samples from $p(\theta)$.
2. Sort them and take the $0.95S^{th}$ ordered value.
   Typically, $S = 1000$ is sufficient.

Simulations also reveal modeling issues — for instance, if extremely large or small values occur, this may indicate poor prior choices or model mis-specification.

---

Generating Random Samples

Modern computation uses pseudorandom number generators (PRNGs) that generate deterministic but statistically uniform sequences in $[0,1]$.
These serve as a base for transforming uniform samples into draws from any distribution through mathematical transformation rules.

---

Sampling Using the Inverse Cumulative Distribution Function (CDF)

The inverse CDF (or inverse transform) method allows sampling from both discrete and continuous distributions.

Given a CDF $F(v)$:

$F(v^*) = P(v \le v^*) = \begin{cases} \sum_{v \le v^*} p(v), & \text{if discrete}, \\ \int_{-\infty}^{v^*} p(v)\, dv, & \text{if continuous.} \end{cases}$

Algorithm

1. Draw $U \sim \text{Uniform}(0,1)$.
2. Compute $v = F^{-1}(U)$.

This guarantees that $v$ is distributed according to $p(v)$.

Example — Exponential Distribution

For $v \sim \text{Exponential}(\lambda)$:

$F(v) = 1 – e^{-\lambda v}.$

Setting $U = F(v)$ gives:

$v = -\frac{\log(1 – U)}{\lambda}.$

Since $1 – U$ is also uniformly distributed:

$v = -\frac{\log U}{\lambda}.$

Hence, exponential random variables can be generated by transforming uniform random draws.

This inverse transform sampling method is simple, general, and forms the foundation for more advanced sampling techniques (e.g., Gibbs sampling, Metropolis–Hastings).

---

Simulation of Posterior and Posterior Predictive Quantities

In Bayesian inference, two main types of simulation are performed:

1. Posterior Simulation — draws from $p(\theta | y)$, representing parameter uncertainty.
2. Posterior Predictive Simulation — draws from $p(\tilde{y} | y)$, representing predictive uncertainty for new or unobserved data.

Assume $S$ draws from the joint posterior:

$(\theta^s, \tilde{y}^s), \quad s = 1, 2, \ldots, S.$

These are stored in a structured table:

Simulation Draw	Parameters (θ₁…θₖ)	Predictive Quantities (ỹ₁…ỹₙ)
1	θ₁¹ … θₖ¹	ỹ₁¹ … ỹₙ¹
2	θ₁² … θₖ²	ỹ₁² … ỹₙ²
…	…	…
S	θ₁ˢ … θₖˢ	ỹ₁ˢ … ỹₙˢ

---

Using Simulated Draws

From these simulated draws, we can compute any derived quantity:

• Transformed parameters: e.g., $\theta_1 / \theta_3$​.
• Posterior probabilities: e.g., $P(\tilde{y}_1 + \tilde{y}_2 > e^{\theta_1})$.
• Credible intervals:
  • For instance, a 95% credible interval for $\theta_jθj​ satisfies
    • $P(\theta_j < a) = 0.025, \quad P(\theta_j > b) = 0.025.$
    • These can be estimated from simulation order statistics, e.g., the 25th and 975th values when $S = 1000$.

Typically, summaries are presented using 50% and 95% posterior intervals, providing central and wide uncertainty bands.

---

Advantages of Simulation-Based Inference

1. General Applicability – works even when no closed-form solution exists.
2. Transparency – simulated samples provide intuitive visual understanding of uncertainty.
3. Model Diagnostics – unrealistic simulated values indicate potential model problems.
4. Ease of Implementation – modern tools like R and Stan simplify posterior simulation.
5. Scalability – easily extends from small models to large hierarchical systems.

---

Accuracy of Simulation

The precision of simulation-based inference depends on $S$, the number of posterior draws.

• $S = 1000$: sufficient for typical summaries (mean, median, 95th percentile).
• $S = 10{,}000+$: recommended for estimating tail probabilities or joint posteriors.

The accuracy of simulation results is further examined, after readers gain experience applying posterior simulation to practical examples.

---

Overall Perspective

The computational foundation of Bayesian data analysis integrates:

• R for preprocessing, visualization, and model diagnostics,
• Stan for efficient sampling and model fitting, and
• Simulation as the bridge between model, data, and inference.

By combining algebraic computation, probability modeling, and simulation, Bayesian computation provides a flexible, interpretable, and empirically grounded framework for handling uncertainty in statistical analysis.