# 1. Foundations ## 1.1 What is Bayesian inference Bayesian inference means updating one's belief about something as new information becomes available. Instead of returning a single fixed estimate, it returns a **degree of belief** — a full probability distribution over possible values — rather than a fixed proportion of outcomes the way the frequentist approach does. Everything starts from conditional probability: $P(A|B) = \frac{P(A \cap B)}{P(B)}$ Combining the definition of conditional probability for $P(A|B)$ and $P(B|A)$ gives **Bayes' theorem**: $P(A|B) = \frac{P(B|A)*P(A)}{P(B)}$ Where $P(B)$ (the normalizing constant, or "evidence") can be expanded using the law of total probability: $P(B) = P(B|A)*P(A) + P(B|\lnot A)*P(\lnot A)$ Intuition: $P(A)$ is what you believed before seeing evidence $B$ (the **prior**). $P(B|A)$ tells you how likely the evidence is if $A$ were true (the **likelihood**). Multiplying the two and normalizing by $P(B)$ tells you what to believe about $A$ *after* seeing the evidence (the **posterior**). Bayesian inference is just this update rule applied repeatedly as data arrives — today's posterior can become tomorrow's prior. ## 1.2 Bayesian vs frequentist | | Frequentist | Bayesian | | ----------------- | ----------------------------------------------- | -------------------------------------------------- | | Parameters | fixed, unknown constants | random variables with distributions | | Output | point estimate + p-value / CI | full posterior distribution | | Interpretation | long-run frequency over repeated sampling | degree of belief given the data observed | | Prior knowledge | not formally incorporated | explicitly incorporated via the prior | | Small samples | asymptotic assumptions can break down | still valid, though more sensitive to prior choice | | Question answered | "how surprising is this data if $H_0$ is true?" | "how probable is each hypothesis given the data?" | ## 1.3 Why use Bayesian inference - Natural handling of **uncertainty** — because parameters have distributions rather than single values, uncertainty propagates automatically through every downstream calculation (predictions, differences between groups, business metrics, etc.) - Possibility to **incorporate expert opinion or domain knowledge** into the model through the prior - No need to rely on fixed constants such as p-values, and no need to fix parameters as constants — they can be modeled as random variables - Statistically coherent even with **little data**, because the prior regularizes the estimate when the data alone is not very informative - **More flexible** than the frequentist approach for building custom, hierarchical, or multi-level models # 2. The Three Ingredients: Prior, Likelihood, Posterior Applying Bayes' theorem to a modeling problem, where $\theta$ represents the unknown parameters and $y$ represents observed data: $P(\theta|y) = \frac{P(y|\theta)*P(\theta)}{P(y)}$ | Term | Name | Meaning | | ------------------- | ------------------------------ | -------------------------------------------------------------------------------------------------------------------------- | | $P(\theta \vert y)$ | posterior | what we know about the parameters after seeing the data | | $P(\theta)$ | prior | what we know about the parameters before seeing any data | | $P(y \vert \theta)$ | likelihood | how probable the observed data is, for a given value of the parameters (the data-generating process) | | $P(y)$ | evidence / marginal likelihood | scaling constant that makes the posterior integrate to 1; obtained by integrating the numerator over all possible $\theta$ | ## 2.1 Prior distribution The prior distribution reflects what we know about the parameter **before observing any data**. It can be chosen based on expert opinion, common knowledge, previous data/research, or subjective belief. Broadly, priors fall on a spectrum: - **Informative prior**: strongly encodes existing knowledge (e.g. a previous experiment's posterior) - **Weakly informative prior**: rules out implausible values but otherwise lets the data dominate - **Non-informative / flat prior**: intended to have minimal influence on the posterior (e.g. uniform) Prior choice can meaningfully impact the posterior, **especially with little data** — as sample size grows, the likelihood dominates and the influence of the prior shrinks. To avoid cherry-picking, prior choices should be clearly stated and explainable based on previous research or sensible domain assumptions (a "prior sensitivity analysis" — refitting with a few different reasonable priors — is good practice to confirm conclusions aren't an artifact of the prior). ## 2.2 Likelihood The likelihood is not chosen *because* you are doing Bayesian inference — it is chosen because it matches the type of data-generating process you are modeling. The same likelihood functions are used in frequentist maximum-likelihood estimation; Bayesian inference just multiplies the likelihood by a prior instead of maximizing it alone. | outcome type | likelihood | | --------------------------------------------------- | --------------------- | | continuous | normal | | continuous with outliers | student t | | binary | bernoulli | | conversion rate (counts of successes out of trials) | binomial | | counts | poisson | | overdispersed counts | negative binomial | | positive skewed continuous | gamma / log-normal | | proportions in (0,1) | beta | | time to event | weibull / exponential | ## 2.3 Posterior distribution The posterior is the actual output of Bayesian inference: a full probability distribution over the parameter(s), combining what the prior said with what the data says. Because it's a distribution rather than a single number, it directly encodes uncertainty — narrow posteriors mean high confidence, wide posteriors mean the data (and prior) haven't pinned the parameter down tightly. A key practical advantage: instead of only testing whether $\mu_1 = \mu_2$ (as a frequentist hypothesis test would), Bayesian inference directly estimates the distribution of the **difference** $\mu_2 - \mu_1$, which can then answer richer questions like "what's the probability $\mu_2$ is at least 5% larger than $\mu_1$?" # 3. Computing the Posterior There are three general strategies for getting from prior + likelihood to posterior, in increasing order of practicality for real-world models: ## 3.1 Grid approximation Discretize the parameter space into a fine grid of candidate values, compute the (unnormalized) posterior density $P(y|\theta)*P(\theta)$ at each grid point, then normalize so the values sum to 1. This is intuitive and exact for simple problems, but **impractical with many parameters** — the number of grid points needed grows exponentially with the number of dimensions (the curse of dimensionality). ## 3.2 Conjugate priors If the prior and likelihood belong to a matching "conjugate" pair, the posterior has a closed-form solution in the same distribution family as the prior — no simulation needed. The classic example is **Beta-Binomial conjugacy**: if the prior is $Beta(\alpha, \beta)$ and the data is a binomial count of successes/failures, the posterior is simply $Beta(\alpha + successes,\ \beta + failures)$. This is fast and elegant but requires the prior to be an "unintuitive" mathematical shape that happens to be conjugate to the likelihood, which limits flexibility for custom models. ## 3.3 Markov Chain Monte Carlo (MCMC) The standard approach used in production Bayesian data analysis, because it works for **any prior and any likelihood**, including non-conjugate combinations and models with many parameters. **Markov chain**: a stochastic simulation process where the probability of moving to the next state depends only on the current state (not the full history), defined by a transition matrix of movement probabilities between states. Markov chains are widely used in retention modeling, attribution analysis, recommendation systems, and — critically here — as the foundation of MCMC sampling for Bayesian inference. Example — transition matrix for "what will the bear do next": | | hunt | eat | sleep | | ----- | ---- | ---- | ----- | | hunt | 0.1 | 0.8 | 0.1 | | eat | 0.05 | 0.4 | 0.55 | | sleep | 0.8 | 0.15 | 0.05 | After many time steps, the chain forgets its starting point and the transition probabilities converge to a **stationary distribution** — the same regardless of which state you started in: | | hunt | eat | sleep | | ----- | ---- | ---- | ----- | | hunt | 0.28 | 0.44 | 0.28 | | eat | 0.28 | 0.44 | 0.28 | | sleep | 0.28 | 0.44 | 0.28 | This is the core idea MCMC exploits: construct a Markov chain whose stationary distribution *is* the posterior you want, **then simulate the chain long enough** that its samples approximate draws from that posterior — even when you can't write the posterior down in closed form. **Monte Carlo estimation** (the "MC" half of MCMC): approximate a quantity by generating random numbers. Classic example — estimating $\pi$: draw a 10x10 square around a circle of radius $r=5$, sample 25 random points uniformly inside the square, count how many land inside the circle (say `xx` of them), and the circle's area is approximately `xx/25 * 100`, which lets you back out $\pi r^2 \approx 78.5$. MCMC applies this same "estimate by random sampling" logic, but the random walk is guided by the Markov chain so that it spends more time in high-probability regions of the posterior. ``` self.idata = pm.sample( draws=draws, tune=tune, chains=chains, cores=cores, target_accept=target_accept, random_seed=self.random_state, return_inferencedata=True, progressbar=progressbar, idata_kwargs={ "log_likelihood": True, }, ) ``` This draws samples from posterior distribution: $P(\theta|y)$. `draws` = number of posterior samples kept per chain after tuning. `tune` = number of warm-up samples per chain. during tuning, the sampler learns suitable step size, posterior geometry, mass matrix, how efficiently to explore the posterior then tuning samples are normally discarded. `chains` = number of independent MCMC chain. multiple chains help evaluate convergence. If they all reach and explore the same posterior distribution, $\hat{R}$ should be close to 1. `cores` = number of CPU cores run chains in parallel. `target_accept` = acceptance probability for NUTS, a high value cause the sampler to use smaller more careful steps. (0.8 is default and faster, 0.9 is common safe, 0.95 is useful when divergences occur, 0.99 is very cautious and slower) ### Convergence diagnostics Because MCMC samples are drawn sequentially from a chain, consecutive samples are **not independent** — each one is correlated with its neighbors. This means raw draw count overstates how much information you actually have. - **Effective Sample Size (ESS)** tells you how many *independent* samples the chain's draws are actually worth. - **ess_bulk**: effective sample size for the **center** of the posterior distribution. E.g. `ess_bulk=901` means that although MCMC generated thousands of raw draws, the information content is roughly equivalent to 901 **independent samples** for estimating central tendency (mean/median). - **ess_tail**: effective sample size for the **tails** of the distribution — relevant for how reliable your **uncertainty** intervals are. E.g. `ess_tail=992` means the tail-based uncertainty intervals are based on roughly 992 independent samples. - Rules of thumb: ESS > 1000 (excellent), ESS > 400 (good), ESS < 100 (concerning), ESS < 50 (problematic). When ESS > 1000: chains have converged, there are enough effective samples, posterior estimates are stable, and credible intervals are reliable. - **Low ESS — causes**: too little data, an overly complex model, strong multicollinearity between predictors, poorly chosen priors, non-identifiability (e.g. a treatment/test indicator that only appears within specific subgroups, so its effect can't be cleanly separated from the group effect), or outliers. - **Low ESS — fixes**: increase the number of draws (`pm.sample(draws, tune)`), use a stronger/more regularizing prior (e.g. `pm.Normal('beta', 0, 1)`), remove collinear variables, or simplify the model. - **Trace plots** (`pm.traceplot()`) visualize each chain's sampled values over iterations for every parameter — multiple chains are run in parallel, and when their traces overlap closely and look like "fuzzy caterpillars" with no trend, that's evidence of good convergence. Divergent-looking or trending traces suggest the chains haven't converged and results shouldn't be trusted yet. ![[Screenshot 2026-07-16 at 11.28.57 AM.png]] # 4. The Beta Distribution Often used in experimentation, Bayesian statistics, and conversion-rate modeling. **The Beta distribution is defined only between 0 and 1**, which makes it perfect for modeling conversion rate, CTR, retention rate, or any other proportion. $X \sim Beta(\alpha,\beta)$ - $\alpha$ = number of successes + prior successes - $\beta$ = number of failures + prior failures Because $\alpha$ and $\beta$ have this "successes/failures" interpretation, the Beta distribution is the **conjugate prior for the Binomial/Bernoulli likelihood** (see §3.2) — starting from $Beta(\alpha_0, \beta_0)$ and observing $s$ successes and $f$ failures gives a posterior of exactly $Beta(\alpha_0+s, \beta_0+f)$, with no simulation required. The shape of the Beta distribution is very flexible depending on $\alpha$ and $\beta$: - $Beta(10,10)$ → symmetric bell shape, centered at 0.5 - $Beta(2,10)$ → right-skewed, mass concentrated near low values - $Beta(1,1)$ → uniform shape (a flat, non-informative prior over [0,1]) Use `np.random.beta(a, b, size)` to generate sample draws from a Beta distribution. # 5. Interpreting Bayesian Results Even though Bayesian inference outputs a full distribution rather than a single number, you still typically want to summarize it: - **Point summaries**: `np.mean()`, `np.median()`, `np.percentile(75)` on the posterior samples - **Credible interval**: the Bayesian analogue of a confidence interval, but with a much more intuitive interpretation — a 95% credible interval means there is a 95% probability the true parameter lies within that range, *given the observed data and prior*. (A frequentist 95% CI instead means: if you repeated the experiment infinitely many times, 95% of the intervals constructed this way would contain the true parameter — it says nothing about the probability for this specific interval.) - **HDI (Highest Density Interval)**: `arviz.hdi(data, prob=0.95)` computes the narrowest interval containing 95% of the posterior mass — unlike a simple equal-tailed interval (e.g. the 2.5th–97.5th percentile), the HDI is the *shortest* interval with that coverage, which matters for skewed or multimodal posteriors. - **Density plots**: `sns.kdeplot(data, shade=True, label='prior')` visualizes the shape of a distribution (prior, likelihood, or posterior) directly. - **ROPE (Region of Practical Equivalence)**: a range around zero (or around "no meaningful difference") defined by domain judgment — if the credible interval falls entirely inside the ROPE, the effect is practically negligible even if it's technically nonzero; if it falls entirely outside, the effect is practically meaningful. # 6. Bayesian A/B Testing ## 6.1 Frequentist approach Based on hypothesis testing: checks whether A and B perform the same or not. Its key limitation is that it does **not** say how much better A is than B — only whether the observed difference is unlikely to have occurred by chance under the null hypothesis. ## 6.2 Bayesian approach Calculate the posterior distribution for each group separately, then compare them directly. This directly answers "what's the probability A is better than B?" and quantifies *how much* better it is — not just whether a difference exists. ``` diff = a_posterior - b_posterior prob_b_better = (diff > 0).mean() ``` Example output: ``` probability of B is better than A: 0.9875 95% CI of B is better than A: [0.00185995 0.0218329] ``` Further steps in a Bayesian A/B test: - **Expected loss**: estimate the cost of making the wrong decision — `loss = diff[diff < 0]` then `loss.mean()` gives the expected magnitude of loss in the scenario where the "winning" variant actually underperforms. - **Business translation**: apply the Bayesian result to business metrics (e.g. conversion lift × revenue per conversion) to run scenario analysis and translate statistical uncertainty into dollar-value uncertainty. - **Visualization**: `pm.forestplot(ads_costs, hdi_prob=0.99)` plots credible intervals for multiple groups/parameters side by side for easy comparison. ![[Pasted image 20260608075221.png|263]] ### Example - **step1. Prior** suppose compare two recovery score, $\mu_c = normal(50,20^2)$ and $\mu_t = normal(60,2^2)$. $sd_c = 20$ so this is uncertain about central point. for the variance, suppose $\sigma = halfnormal(10)$. if using student-t, then $v=Exponential(\frac{1}{29})+2$ - **step2. Likelihood** suppose we have observed data: ``` control = [52,49,50,53,47], mean_c = 50.2 treatment = [54,57,52,55,58], mean_t = 55.2 ``` the likelihood is $y_i = studentT(\mu,\sigma,v)$ or $y_i=normal(\mu,\sigma)$. to simplify, assume observed SD is $\sigma=3$, for each individual points likelihood: $likelihood_{ci} = y_{ci} | \mu_c = normal(\mu_c,3^2)$ $likelihood_{ti} = y_{ti} | \mu_t = normal(\mu_t,3^2) = \frac{1}{\sqrt{2\pi(3^2)}}exp(-\frac{(y_{ti}-\mu_{t})^2}{2(3^2)})$ $likelihood_{all} = \prod{p(y_{ti}|\mu_t)} = n(54|\mu_t,9)\times n(57|\mu_t,9)\times n(52|\mu_t,9)\times n(55|\mu_t,9)\times n(58|\mu_t,9)$ - **step3. Posterior** $Posterior \approx Likelihood \times Prior$, now PyMC gives thousands of draws and posterior draw - $precision_{priort} = \frac{1}{variance} = \frac{1}{4}$ - $precision_{data} = \frac{1}{3^2}\times 5$, 5 observations. then $\sigma = {1}/{precision_{data}}$ - $\mu_t = normal(60,2^2)$ and treatment likelihood is $\bar{y_t}=55.2$ and $\bar{y_t} | \mu_t = normal(\mu_t,\frac{1}{precision_{data}}) = nomal(\mu_t,1.8)$ - $precision_{postt} = precision_{priort} + precision_{data} = 0.8056$ and $\sigma_{post} = 1/0.8056=1.2414$ - $\mu_{postt} = \frac{precision_{priort} \times \mu_{prior} + precision_{data} \times \mu_{data}}{precision_{post}} = \frac{0.25\times 60 + 0.5556\times 55.2}{0.8056} = 56.69$ - $\mu_{postt} | y_t = normal(56.69,1.2414)$ - precision works as weight, $\mu_{postt}=\frac{0.25}{0.8056}\times 60 + \frac{0.5556}{0.8056}\times 55.2=56.69$ - from the calculation, the broader the prior distribution, the less weight it is in post distribution which means post distribution still rely on observed data - then we got posterior distribution, `pymc` draws follow posterior distribution to generate new data ``` mu_control = [49.8,50.3,50.1,49.9] mu_treatment = [53.7,54.1,53.5,53.9] difference = [3.9,3.8,3.4,4.0] ``` - **step4. Credible interval** suppose there are 8000 posterior draws and 8000 difference then get 95% hdi, `P(treat>control) = np.mean(difference>0)` - **step5. Practical significance** suppose business need difference > 2 then `P(practical difference) = np.mean(difference>2)` - **overall**. Bayesian don't fully reply on prior nor observed data, but balance them and calculate a posterior distribution, then generate possible data to do significant test # 7. Bayesian Regression A **Bayesian regression model** treats regression parameters (like slopes and intercepts) as probability distributions rather than fixed, single values. It combines prior beliefs with observed data using Bayes' theorem to output a full distribution of predictions, actively quantifying prediction uncertainty (rather than reporting a single point estimate with a separately-computed standard error). ## 7.1 Estimating regression parameters - Grid approximation → impractical with many parameters (§3.1) - Choose conjugate priors and simulate from a known posterior → fast, but forces unintuitive prior choices (§3.2) - Simulate from the posterior via MCMC, works even with non-conjugate priors → standard approach in practice (§3.3) - Assume the parameter draws are given (e.g. reusing posterior draws from a previous fit) ## 7.2 Tooling - `scikit-learn`'s `linear_model` module offers `BayesianRidge` for a quick, out-of-the-box approach with limited flexibility - For full probabilistic modeling with custom priors/likelihoods, **PyMC** or **Stan** are the standard choices ## 7.3 Reading the output Bayesian regression outputs a distribution for every coefficient, e.g. $\beta_0 \sim N(5,2)$ — visualize with `pm.plot_posterior()` or compare coefficients across a model with `pm.forestplot()`. The regression's expected point estimate for a new observation is still computed the familiar way, just using the posterior mean (or full posterior) of each coefficient: $\hat{y} = intercept + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_n x_n$ $sd_{coef}$ (the posterior standard deviation of a coefficient) is used to randomly sample plausible coefficient values around the mean, which is how prediction uncertainty gets propagated into forecasts. # 8. Hierarchical Bayesian Models Hierarchical Bayesian modeling estimates the effect of a strategy (or treatment) while accounting for the fact that observations belong to **different groups** — such as tickers and market conditions — that naturally differ in their baseline performance. Ignoring this grouping risks **Simpson's paradox**: an aggregate result that contradicts what every subgroup individually shows, often because groups were imbalanced in sample size or which treatment they received. **Worked example** — comparing Strategy A vs Strategy B across stocks: | Ticker | strategy | return lift | | ------ | -------- | ----------- | | AAPL | A | 1% | | AAPL | A | 2% | | AAPL | B | 3% | | AAPL | B | 4% | | TSLA | A | -2% | | TSLA | B | -1% | - **Simple A/B test (ignores ticker)**: $A=(1+2-2)/3=0.33,\ B=(3+4-1)/3=2,\ diff=B-A=1.67$ → "Strategy B is better by 1.67%" - **Looking within each group**: AAPL → $A=(1+2)/2=1.5,\ B=(3+4)/2=3.5,\ diff=2$. TSLA → $diff=B-A=(-1)-(-2)=1$. A naive average across groups (ignoring sample size) would give AAPL $diff=2.5$ and TSLA $diff=-1.5$ if computed the wrong way — ordinary regression treats AAPL and TSLA equally, ignoring the fact that AAPL has more observations and a different baseline. - **Hierarchical model**: uses a prior distribution like $\beta_{ticker} \sim Normal(0, \sigma_{ticker})$, so each ticker's effect is "shrunk" toward the overall average by an amount that depends on how much data that ticker has (small-sample groups get pulled more strongly toward the population mean — this is called **partial pooling**, and it prevents small groups from being overfit). As a result, the model might land on AAPL = 1.9, TSLA = -0.7 — informed by, but not purely equal to, either the raw subgroup averages or the pooled overall average. When we say the model "controls for ticker and market," we mean it isolates the strategy effect after **accounting for differences** attributable to the stock itself or the market environment. A positive `beta_strategy` coefficient therefore means Strategy B still improves return lift even after adjusting for the fact that some tickers and market regimes naturally have higher or lower returns — making the estimated strategy effect much more credible, and helping explain apparent Simpson's paradoxes where the aggregate result differs from subgroup results. # 9. Model Diagnostics and Comparison - `pm.traceplot()` / `pm.forestplot()` visualize each parameter's posterior distribution, showing multiple Markov chains' sampled draws — when the chains' lines overlap closely, that indicates good convergence (see §3.3 for what to do if they don't). - To compare different candidate models, use `pm.compare()` or `pm.summary()` to compare **WAIC** (Widely Applicable Information Criterion — an estimate of out-of-sample predictive accuracy that penalizes model complexity, analogous to AIC but valid for Bayesian posterior distributions), or visualize with `pm.compareplot()`. Lower WAIC generally indicates better expected predictive performance. # 10. Making Predictions from a Bayesian Regression To sample predictive draws: build the model with `pm.GLM.from_formula()`, then draw from the posterior predictive distribution with `pm.fast_sample_posterior_predictive()`, and inspect the result with `pm.plot_posterior()`. To evaluate how good the model is, plot the distribution of test errors and check whether it's centered around 0 — `az.hdi()` gives the credible interval for that error distribution. Core PyMC building blocks used throughout: `pm.Model()`, `pm.sample()`, `pm.Normal()`. # 11. Projects and Resources - [Bayesian project](https://colab.research.google.com/drive/1S9zOL6NNaNnaCZaaohwyKF4hc1z9R31N#scrollTo=jWNolr-JhTC5) - [Bayesian framework](https://app.hex.tech/019eb406-8c5e-755b-986f-e0693fc2a83c/hex/Bayesian-Test-Playground-033pzz7HLh5xNZYyKaBqdv/draft/logic)