Bayesian forecasting treats predictions as probability distributions that update as new evidence arrives. Rather than producing a single point estimate, the Bayesian approach maintains a full belief distribution — a prior — and systematically refines it with each new piece of data using Bayes' theorem. This produces calibrated uncertainty estimates, naturally handles sparse data, and provides a principled framework for sequential updating that is essential for real-time forecasting.

## Key Points

- **Prior**: Your initial forecast before new data
- **Likelihood**: How probable the observed data is under different scenarios
- **Posterior**: Your updated forecast after incorporating new data
1. Bayesian forecasting maintains full probability distributions, not just point estimates, giving calibrated uncertainty
2. Conjugate priors (Beta-Binomial, Normal-Normal, Gamma-Poisson) provide analytical updates without computation
3. Prior selection encodes domain knowledge; uninformative priors let data dominate, informative priors regularize with sparse data
4. Sequential updating via Kalman-filter-style methods enables real-time forecast refinement as evidence arrives
5. Bayesian networks model causal relationships, enabling "what if" reasoning for prediction
6. Calibration tracking is essential: regularly check whether your 80% confidence intervals contain the truth 80% of the time
7. PyMC and similar tools make complex Bayesian models accessible without deriving posteriors analytically
8. The information value calculation helps decide which data to collect next — prioritize observations that maximally reduce uncertainty

Bayesian Forecasting

Overview

Bayesian forecasting treats predictions as probability distributions that update as new evidence arrives. Rather than producing a single point estimate, the Bayesian approach maintains a full belief distribution — a prior — and systematically refines it with each new piece of data using Bayes' theorem. This produces calibrated uncertainty estimates, naturally handles sparse data, and provides a principled framework for sequential updating that is essential for real-time forecasting.

Bayes' Theorem

The Foundation

P(H|E) = P(E|H) * P(H) / P(E)

Where:
  P(H|E)  = Posterior: probability of hypothesis H given evidence E
  P(E|H)  = Likelihood: probability of seeing evidence E if H is true
  P(H)    = Prior: probability of H before seeing evidence
  P(E)    = Evidence: total probability of seeing this evidence

In forecasting terms:

• Prior: Your initial forecast before new data
• Likelihood: How probable the observed data is under different scenarios
• Posterior: Your updated forecast after incorporating new data

Simple Example: Election Forecasting

def bayesian_election_update(
    prior_prob_win: float,     # Prior probability candidate wins
    poll_result: float,         # Poll shows X% support
    poll_accuracy: float,       # Historical poll accuracy (std dev)
    true_support_if_win: float  # Expected poll result if candidate will win
) -> float:
    """
    Update election win probability based on a new poll.

    Example:
        prior = 0.60 (60% chance of winning)
        poll shows 52% support
        polls have ~3% std dev
        winning candidates typically poll at 54%
    """
    from scipy.stats import norm

    # Likelihood: P(poll_result | candidate wins)
    likelihood_win = norm.pdf(poll_result, loc=true_support_if_win, scale=poll_accuracy)

    # Likelihood: P(poll_result | candidate loses)
    true_support_if_lose = 1.0 - true_support_if_win  # Simplified
    likelihood_lose = norm.pdf(poll_result, loc=true_support_if_lose, scale=poll_accuracy)

    # Posterior via Bayes
    numerator = likelihood_win * prior_prob_win
    denominator = (likelihood_win * prior_prob_win +
                   likelihood_lose * (1 - prior_prob_win))

    posterior = numerator / denominator
    return posterior

# Example usage
posterior = bayesian_election_update(
    prior_prob_win=0.60,
    poll_result=0.52,
    poll_accuracy=0.03,
    true_support_if_win=0.54
)
print(f"Updated probability: {posterior:.1%}")

Prior Selection

Types of Priors

class PriorSelector:
    """Guide for selecting appropriate prior distributions."""

    @staticmethod
    def uninformative_prior(param_type: str) -> dict:
        """
        Minimally informative priors when you have little prior knowledge.
        """
        priors = {
            'proportion': {
                'distribution': 'Beta(1, 1)',  # Uniform on [0,1]
                'description': 'Equal weight to all proportions',
                'scipy': 'stats.beta(1, 1)'
            },
            'positive_real': {
                'distribution': 'HalfCauchy(0, 10)',
                'description': 'Weakly informative for scale parameters',
                'scipy': 'stats.halfcauchy(scale=10)'
            },
            'real_valued': {
                'distribution': 'Normal(0, 100)',
                'description': 'Vague normal centered at zero',
                'scipy': 'stats.norm(0, 100)'
            },
            'count': {
                'distribution': 'Gamma(0.01, 0.01)',
                'description': 'Vague prior for rate parameters',
                'scipy': 'stats.gamma(0.01, scale=100)'
            }
        }
        return priors.get(param_type, priors['real_valued'])

    @staticmethod
    def informative_prior(domain: str, parameter: str) -> dict:
        """
        Domain-informed priors encoding real knowledge.
        """
        priors = {
            ('marketing', 'conversion_rate'): {
                'distribution': 'Beta(2, 50)',
                'description': 'Typical web conversion rate ~2-5%',
                'mean': 2/52, 'reasoning': 'Industry benchmarks'
            },
            ('finance', 'annual_return'): {
                'distribution': 'Normal(0.07, 0.15)',
                'description': 'Historical equity market return',
                'mean': 0.07, 'reasoning': 'Long-run market data'
            },
            ('software', 'defect_rate'): {
                'distribution': 'Gamma(2, 5)',
                'description': 'Defects per KLOC',
                'mean': 10, 'reasoning': 'Industry CMMI data'
            }
        }
        return priors.get((domain, parameter), {})

    @staticmethod
    def prior_from_expert_elicitation(
        median: float, p25: float, p75: float
    ) -> dict:
        """
        Construct a prior from expert-elicited quantiles.
        Ask expert: "What is the median? What range captures the middle 50%?"
        """
        import numpy as np
        from scipy.optimize import minimize

        def loss(params):
            mu, sigma = params
            from scipy.stats import norm
            predicted_median = norm.ppf(0.5, mu, sigma)
            predicted_p25 = norm.ppf(0.25, mu, sigma)
            predicted_p75 = norm.ppf(0.75, mu, sigma)
            return (
                (predicted_median - median)**2 +
                (predicted_p25 - p25)**2 +
                (predicted_p75 - p75)**2
            )

        result = minimize(loss, x0=[median, (p75 - p25) / 1.35])
        mu, sigma = result.x

        return {
            'distribution': f'Normal({mu:.3f}, {sigma:.3f})',
            'mean': mu,
            'std': sigma,
            'description': f'Fitted to expert quantiles: median={median}, IQR=[{p25}, {p75}]'
        }

Conjugate Priors

When Prior and Posterior Have the Same Form

Conjugate priors make Bayesian updating analytically tractable:

class ConjugatePriors:
    """Conjugate prior-likelihood pairs for common scenarios."""

    @staticmethod
    def beta_binomial(
        prior_alpha: float, prior_beta: float,
        successes: int, failures: int
    ) -> dict:
        """
        Beta prior + Binomial likelihood = Beta posterior.
        Perfect for: conversion rates, click-through rates, win rates.

        Example: Estimating a website's conversion rate
            Prior: Beta(2, 50) — expect ~4% conversion
            Data: 15 conversions out of 200 visitors
            Posterior: Beta(2+15, 50+185) = Beta(17, 235)
        """
        post_alpha = prior_alpha + successes
        post_beta = prior_beta + failures

        return {
            'posterior': f'Beta({post_alpha}, {post_beta})',
            'mean': post_alpha / (post_alpha + post_beta),
            'mode': (post_alpha - 1) / (post_alpha + post_beta - 2) if (post_alpha > 1 and post_beta > 1) else None,
            'variance': (post_alpha * post_beta) / ((post_alpha + post_beta)**2 * (post_alpha + post_beta + 1)),
            'ci_95': _beta_ci(post_alpha, post_beta, 0.95)
        }

    @staticmethod
    def normal_normal(
        prior_mean: float, prior_var: float,
        data_mean: float, data_var: float, n: int
    ) -> dict:
        """
        Normal prior + Normal likelihood = Normal posterior.
        Perfect for: continuous measurements, forecasting values.

        The posterior precision = prior precision + data precision.
        """
        prior_precision = 1 / prior_var
        data_precision = n / data_var

        post_precision = prior_precision + data_precision
        post_var = 1 / post_precision
        post_mean = (prior_mean * prior_precision + data_mean * data_precision) / post_precision

        return {
            'posterior': f'Normal({post_mean:.4f}, {post_var:.6f})',
            'mean': post_mean,
            'variance': post_var,
            'std': post_var ** 0.5,
            'ci_95': (post_mean - 1.96 * post_var**0.5, post_mean + 1.96 * post_var**0.5),
            'effective_sample_weight': data_precision / post_precision
        }

    @staticmethod
    def gamma_poisson(
        prior_alpha: float, prior_beta: float,
        observed_count: int, observation_period: float
    ) -> dict:
        """
        Gamma prior + Poisson likelihood = Gamma posterior.
        Perfect for: event rates, arrival rates, incident counts.
        """
        post_alpha = prior_alpha + observed_count
        post_beta = prior_beta + observation_period

        return {
            'posterior': f'Gamma({post_alpha}, {post_beta})',
            'mean_rate': post_alpha / post_beta,
            'variance': post_alpha / post_beta**2,
            'ci_95': _gamma_ci(post_alpha, post_beta, 0.95)
        }


def _beta_ci(alpha, beta, level):
    from scipy.stats import beta as beta_dist
    low = beta_dist.ppf((1 - level) / 2, alpha, beta)
    high = beta_dist.ppf(1 - (1 - level) / 2, alpha, beta)
    return (low, high)


def _gamma_ci(alpha, beta, level):
    from scipy.stats import gamma
    low = gamma.ppf((1 - level) / 2, alpha, scale=1/beta)
    high = gamma.ppf(1 - (1 - level) / 2, alpha, scale=1/beta)
    return (low, high)

Bayesian Networks

Modeling Causal Relationships for Prediction

class BayesianNetwork:
    """Simple Bayesian network for probabilistic reasoning."""

    def __init__(self):
        self.nodes = {}  # name -> {'parents': [], 'cpt': {}}
        self.evidence = {}

    def add_node(self, name: str, parents: list = None,
                 cpt: dict = None):
        """
        Add a node with its conditional probability table.

        CPT format for binary variables:
        {
            (): 0.3,  # P(node=True) if no parents
            (True,): 0.8,  # P(node=True | parent1=True)
            (False,): 0.2,  # P(node=True | parent1=False)
            (True, True): 0.9,  # P(node=True | parent1=True, parent2=True)
        }
        """
        self.nodes[name] = {
            'parents': parents or [],
            'cpt': cpt or {}
        }

    def set_evidence(self, node: str, value: bool):
        """Set observed evidence."""
        self.evidence[node] = value

    def query(self, target: str, n_samples: int = 10000) -> float:
        """
        Estimate P(target=True | evidence) using likelihood weighting.
        """
        weighted_true = 0
        total_weight = 0

        for _ in range(n_samples):
            sample, weight = self._generate_weighted_sample()
            if weight > 0:
                total_weight += weight
                if sample.get(target, False):
                    weighted_true += weight

        return weighted_true / total_weight if total_weight > 0 else 0.5

    def _generate_weighted_sample(self) -> tuple:
        """Generate a sample using likelihood weighting."""
        sample = {}
        weight = 1.0

        # Topological order (simplified: assume nodes added in order)
        for name, node in self.nodes.items():
            parent_values = tuple(sample[p] for p in node['parents'])
            prob_true = node['cpt'].get(parent_values, 0.5)

            if name in self.evidence:
                sample[name] = self.evidence[name]
                weight *= prob_true if self.evidence[name] else (1 - prob_true)
            else:
                sample[name] = np.random.random() < prob_true

        return sample, weight


# Example: Predicting product launch success
def product_launch_model():
    """Bayesian network for product launch prediction."""
    bn = BayesianNetwork()

    bn.add_node('market_demand', cpt={(): 0.6})
    bn.add_node('competitor_launches', cpt={(): 0.4})
    bn.add_node('engineering_ready', cpt={(): 0.7})

    bn.add_node('launch_timing_good', parents=['market_demand', 'competitor_launches'], cpt={
        (True, True): 0.3,   # Good demand but competitor launching too
        (True, False): 0.9,  # Good demand, no competitor
        (False, True): 0.1,  # Bad demand + competitor
        (False, False): 0.4, # Bad demand, no competitor
    })

    bn.add_node('product_success', parents=['launch_timing_good', 'engineering_ready'], cpt={
        (True, True): 0.85,
        (True, False): 0.3,
        (False, True): 0.4,
        (False, False): 0.1,
    })

    # Set evidence: we know engineering is ready
    bn.set_evidence('engineering_ready', True)

    # Query: what is P(product success | engineering ready)?
    prob_success = bn.query('product_success')
    print(f"P(success | engineering ready) = {prob_success:.1%}")

    # Additional evidence: market demand looks strong
    bn.set_evidence('market_demand', True)
    prob_success_updated = bn.query('product_success')
    print(f"P(success | eng ready, strong demand) = {prob_success_updated:.1%}")

    return bn

Sequential Bayesian Updating

Real-Time Forecast Updates

class SequentialBayesianForecaster:
    """Update forecasts in real-time as new data arrives."""

    def __init__(self, prior_mean: float, prior_std: float):
        self.mean = prior_mean
        self.variance = prior_std ** 2
        self.update_history = [{
            'step': 0,
            'mean': prior_mean,
            'std': prior_std,
            'evidence': 'prior'
        }]

    def update(self, observation: float, observation_noise: float,
               label: str = ""):
        """
        Kalman filter style update for sequential Bayesian forecasting.
        """
        # Kalman gain
        kalman_gain = self.variance / (self.variance + observation_noise**2)

        # Update mean
        innovation = observation - self.mean
        self.mean = self.mean + kalman_gain * innovation

        # Update variance (always decreases with more data)
        self.variance = (1 - kalman_gain) * self.variance

        self.update_history.append({
            'step': len(self.update_history),
            'mean': self.mean,
            'std': self.variance ** 0.5,
            'evidence': label,
            'innovation': innovation,
            'kalman_gain': kalman_gain
        })

        return self.mean, self.variance ** 0.5

    def predict(self, confidence: float = 0.90) -> dict:
        """Get current prediction with confidence interval."""
        from scipy.stats import norm
        z = norm.ppf(1 - (1 - confidence) / 2)
        std = self.variance ** 0.5

        return {
            'point_estimate': self.mean,
            'std': std,
            'ci_lower': self.mean - z * std,
            'ci_upper': self.mean + z * std,
            'confidence': confidence,
            'n_updates': len(self.update_history) - 1
        }

    def information_value(self, potential_observation_noise: float) -> float:
        """
        Estimate how much a new observation would reduce uncertainty.
        Returns expected variance reduction as fraction.
        """
        potential_gain = self.variance / (self.variance + potential_observation_noise**2)
        potential_new_variance = (1 - potential_gain) * self.variance
        reduction = 1 - potential_new_variance / self.variance
        return reduction


# Example: Tracking an election forecast
def election_tracking_example():
    forecaster = SequentialBayesianForecaster(
        prior_mean=0.50,  # Start at 50/50
        prior_std=0.10    # Moderate initial uncertainty
    )

    # Polls come in over time
    polls = [
        (0.52, 0.03, "Ipsos Poll Jan"),
        (0.48, 0.04, "YouGov Poll Jan"),
        (0.51, 0.03, "Pew Poll Feb"),
        (0.54, 0.025, "Large national poll Mar"),
        (0.53, 0.02, "Aggregated poll Mar"),
        (0.55, 0.03, "Post-debate poll Apr"),
    ]

    for poll_result, poll_error, label in polls:
        mean, std = forecaster.update(poll_result, poll_error, label)
        print(f"After {label}: {mean:.1%} +/- {std:.1%}")

    prediction = forecaster.predict(confidence=0.90)
    print(f"\nFinal prediction: {prediction['point_estimate']:.1%}")
    print(f"90% CI: [{prediction['ci_lower']:.1%}, {prediction['ci_upper']:.1%}]")

Calibration for Bayesian Forecasters

Checking Your Calibration

class CalibrationTracker:
    """Track and improve Bayesian calibration over time."""

    def __init__(self):
        self.forecasts = []  # (predicted_prob, actual_outcome)

    def record(self, predicted_probability: float, outcome: bool):
        self.forecasts.append((predicted_probability, outcome))

    def calibration_score(self, n_bins: int = 10) -> dict:
        """Compute calibration metrics."""
        if not self.forecasts:
            return {'error': 'No forecasts recorded'}

        bins = np.linspace(0, 1, n_bins + 1)
        calibration_data = []

        for i in range(n_bins):
            bin_forecasts = [
                (p, o) for p, o in self.forecasts
                if bins[i] <= p < bins[i+1]
            ]
            if bin_forecasts:
                predicted_avg = np.mean([p for p, o in bin_forecasts])
                actual_avg = np.mean([o for p, o in bin_forecasts])
                calibration_data.append({
                    'bin_center': (bins[i] + bins[i+1]) / 2,
                    'predicted': predicted_avg,
                    'actual': actual_avg,
                    'count': len(bin_forecasts),
                    'error': abs(predicted_avg - actual_avg)
                })

        # Brier score
        brier = np.mean([
            (p - float(o))**2 for p, o in self.forecasts
        ])

        # Expected Calibration Error
        total = len(self.forecasts)
        ece = sum(
            (d['count'] / total) * d['error']
            for d in calibration_data
        )

        return {
            'brier_score': brier,
            'ece': ece,
            'calibration_curve': calibration_data,
            'is_well_calibrated': ece < 0.05,
            'overconfident': self._check_overconfidence(),
            'n_forecasts': len(self.forecasts)
        }

    def _check_overconfidence(self) -> bool:
        """Check if forecaster is systematically overconfident."""
        extreme = [(p, o) for p, o in self.forecasts if p > 0.8 or p < 0.2]
        if len(extreme) < 10:
            return False
        high_conf = [(p, o) for p, o in extreme if p > 0.8]
        if high_conf:
            actual_rate = np.mean([o for p, o in high_conf])
            return actual_rate < 0.7  # If 80%+ predictions only come true 70%
        return False

Practical Bayesian Forecasting with PyMC

Using PyMC for Complex Models

# pip install pymc arviz
import pymc as pm
import arviz as az

def bayesian_demand_forecast(historical_sales: list, n_forecast: int = 12):
    """
    Bayesian time series forecast for monthly demand.
    Captures trend, seasonality, and uncertainty.
    """
    n_obs = len(historical_sales)
    sales = np.array(historical_sales)

    with pm.Model() as model:
        # Trend
        trend_slope = pm.Normal('trend_slope', mu=0, sigma=10)
        trend_intercept = pm.Normal('trend_intercept', mu=sales.mean(), sigma=sales.std())

        # Seasonality (12-month cycle)
        season_amplitude = pm.HalfNormal('season_amp', sigma=sales.std() / 2)
        season_phase = pm.Uniform('season_phase', lower=0, upper=2 * np.pi)

        # Time index
        t = np.arange(n_obs)

        # Mean function
        mu = (trend_intercept +
              trend_slope * t +
              season_amplitude * np.sin(2 * np.pi * t / 12 + season_phase))

        # Noise
        sigma = pm.HalfNormal('sigma', sigma=sales.std())

        # Likelihood
        obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=sales)

        # Sample posterior
        trace = pm.sample(2000, tune=1000, cores=2, return_inferencedata=True)

    # Generate forecasts
    t_future = np.arange(n_obs, n_obs + n_forecast)

    with model:
        pm.set_data({})  # Clear observed data
        # Use posterior samples to generate predictions
        post = trace.posterior
        forecast_samples = (
            post['trend_intercept'].values.flatten()[:, None] +
            post['trend_slope'].values.flatten()[:, None] * t_future[None, :] +
            post['season_amp'].values.flatten()[:, None] *
            np.sin(2 * np.pi * t_future[None, :] / 12 +
                   post['season_phase'].values.flatten()[:, None])
        )

    return {
        'forecast_mean': np.mean(forecast_samples, axis=0),
        'forecast_std': np.std(forecast_samples, axis=0),
        'forecast_p10': np.percentile(forecast_samples, 10, axis=0),
        'forecast_p90': np.percentile(forecast_samples, 90, axis=0),
        'trace': trace
    }

Key Takeaways

1. Bayesian forecasting maintains full probability distributions, not just point estimates, giving calibrated uncertainty
2. Conjugate priors (Beta-Binomial, Normal-Normal, Gamma-Poisson) provide analytical updates without computation
3. Prior selection encodes domain knowledge; uninformative priors let data dominate, informative priors regularize with sparse data
4. Sequential updating via Kalman-filter-style methods enables real-time forecast refinement as evidence arrives
5. Bayesian networks model causal relationships, enabling "what if" reasoning for prediction
6. Calibration tracking is essential: regularly check whether your 80% confidence intervals contain the truth 80% of the time
7. PyMC and similar tools make complex Bayesian models accessible without deriving posteriors analytically
8. The information value calculation helps decide which data to collect next — prioritize observations that maximally reduce uncertainty