ƒ

Fugue

A type-safe, monadic probabilistic programming library for Rust — compose models in direct style, listen backwards with state-of-the-art inference.

Fugue Evo ↗Evolution as Bayesian inference — CMA-ES, NSGA-II, island models, and estimation-of-distribution algorithms on the same foundations, with its own live playground.

About Fugue

• 🧩 Monadic PPL: Compose probabilistic programs as pure Model values, then interpret them with pluggable handlers — effects are interpreted, never performed by the model itself
• 🔒 Type-Safe Distributions: 17 built-in distributions with natural return types (Bernoulli → bool, Poisson → u64, Categorical → usize) — the type system tracks what each draw is
• 📊 Multiple Inference Methods: adaptive Metropolis–Hastings, Hamiltonian Monte Carlo, Sequential Monte Carlo (with an unbiased log-evidence estimate), Variational Inference (support-matched guide families), and importance-weighted ABC-SMC (see Advanced Inference)
• 🔍 Comprehensive Diagnostics: split-R̂ convergence, autocorrelation-based effective sample size, Geweke, and closed-form validation
• 🌀 Stack-Safe Interpreter: a trampolined runtime evaluates models in O(1) stack depth — 100 000-site programs run without overflowing
• ⚡ Numerically Stable: log-space computations throughout for robust probability arithmetic
• ✨ Ergonomic Macros: Do-notation (prob!), vectorization (plate!), addressing (addr!)

Installation

[dependencies]
fugue-ppl = "0.2.0"

More Resources

• GitHub Repository
• API Reference on docs.rs
• Zotero Library
• Ask DeepWiki

Getting Started with Fugue

Welcome to Fugue, a type-safe probabilistic programming library for Rust! This guide will get you building Bayesian models in just 15-20 minutes.

Learning Path

We recommend following this path for the best learning experience:

flowchart LR
    A[Installation<br/>2 min] --> B[Your First Model<br/>5 min]
    B --> C[Understanding Models<br/>8 min]
    C --> D[Running Inference<br/>5 min]
    D --> E[Complete Tutorials<br/>45-60 min each]

Quick Start

If you're impatient and want to see Fugue in action immediately:

# Create a new project
cargo new my_bayesian_project
cd my_bayesian_project

# Add Fugue
cargo add fugue-ppl rand

# Copy our "Hello, Probabilistic World!" example into src/main.rs
# (See Installation section)

# Run it!
cargo run

What Makes Fugue Different?

🔒 Type Safety First

Unlike other probabilistic programming libraries, Fugue preserves natural types:

// In Fugue ✅
let coin: bool = sample(addr!("coin"), Bernoulli::new(0.5).unwrap());
let count: u64 = sample(addr!("events"), Poisson::new(3.0).unwrap());
let category: usize = sample(addr!("choice"), Categorical::uniform(5).unwrap());

// Other PPLs ❌
let coin: f64 = sample("coin", Bernoulli(0.5));  // Returns 0.0 or 1.0
let count: f64 = sample("events", Poisson(3.0)); // Need to cast to int
// let category: f64 = sample("choice", Categorical([...])); // Risky indexing

🚀 Zero-Cost Abstractions

Models compile to efficient code with no runtime overhead.

🧰 Composable Architecture

Separate model specification from execution strategy through handlers.

📊 Diagnostics Built In

Split-R̂, autocorrelation-based and multi-chain effective sample size, and a structured error taxonomy.

Architecture Overview

Fugue's modular design separates concerns cleanly:

graph TB
    subgraph "Your Code"
        M[Model Definition]
        D[Data & Observations]
    end

    subgraph "Core System"
        C[Distributions & Types]
        H[Handlers & Interpreters]
        T[Traces]
    end

    subgraph "Inference Engines"
        MCMC[MCMC Sampling]
        HMC[Hamiltonian Monte Carlo]
        SMC[Particle Filtering]
        VI[Variational Inference]
        ABC[ABC Methods]
    end

    M --> C
    D --> C
    C --> H
    H --> T
    T --> MCMC
    T --> HMC
    T --> SMC
    T --> VI
    T --> ABC

The Big Picture

Probabilistic Programming lets you:

1. Model uncertainty and relationships in data
2. Condition on observations to learn parameters
3. Infer posterior distributions and make predictions
4. Quantify uncertainty in your conclusions

Fugue makes this safe, fast, and composable in Rust.

Next Steps

Ready to dive in?

After completing Getting Started, explore:

• Complete Tutorials - End-to-end projects with real applications
• How-To Guides - Specific techniques and best practices
• API Documentation - Comprehensive technical reference

Prerequisites: Basic Rust knowledge (variables, functions, cargo commands)

Installation

Getting Fugue set up in your Rust project takes just 2 minutes. Let's get you running!

# Install Rust via rustup
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh

# Update to latest stable
rustup update stable

Adding Fugue to Your Project

New Project

cargo new my_probabilistic_project
cd my_probabilistic_project

Add Fugue to your Cargo.toml:

[dependencies]
fugue-ppl = "0.2.0"
rand = "0.8"  # For random number generation

Existing Project

Add Fugue to your existing Cargo.toml:

[dependencies]
fugue-ppl = "0.2.0"
rand = "0.8"

Or use cargo add:

cargo add fugue-ppl rand

Verification: "Hello, Probabilistic World!"

Let's verify your installation with a simple example that showcases Fugue's type safety.

Create or replace src/main.rs:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    println!("🎲 Hello, Probabilistic World!");

    // Create a simple model: flip a biased coin
    let coin_model = sample(addr!("coin"), Bernoulli::new(0.7).unwrap());

    // Run the model with a seeded RNG for reproducible results
    let mut rng = StdRng::seed_from_u64(42);
    let (is_heads, trace) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        coin_model,
    );

    // Print the result - notice it's a bool, not a float!
    let result = if is_heads { "Heads" } else { "Tails" };
    println!("🪙 Coin flip result: {}", result);
    println!("📊 Log probability: {:.4}", trace.total_log_weight());

    // Demonstrate type safety with different distributions
    let mut rng = StdRng::seed_from_u64(123);

    // Count events - returns u64 directly
    let (event_count, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("events"), Poisson::new(3.5).unwrap()),
    );
    println!("🎯 Event count: {} (type: u64)", event_count);

    // Choose category - returns usize for safe indexing
    let options = vec!["red", "green", "blue"];
    let (category_idx, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("color"), Categorical::uniform(3).unwrap()),
    );
    println!("🎨 Chosen color: {} (safe indexing!)", options[category_idx]);

    println!("✅ Fugue is working correctly!");
}

Run it to verify everything works:

cargo run

You should see output like:

🎲 Hello, Probabilistic World!
🪙 Coin flip result: Heads
📊 Log probability: -0.3567
🎯 Event count: 4 (type: u64)
🎨 Chosen color: blue (safe indexing!)
✅ Fugue is working correctly!

IDE Setup

VS Code

Install the rust-analyzer extension for the best development experience:

1. Open VS Code
2. Go to Extensions (Ctrl+Shift+X)
3. Search for "rust-analyzer"
4. Install the official rust-analyzer extension

Other IDEs

• IntelliJ/CLion: Install the Rust plugin
• Vim/Neovim: Use coc-rust-analyzer or native LSP
• Emacs: Use lsp-mode with rust-analyzer

Optional: Running Examples

Fugue comes with comprehensive examples to explore:

# Clone the repository to access examples
git clone https://github.com/alexnodeland/fugue
cd fugue

# List available examples
ls examples/

# Run a simple example
cargo run --example bayesian_coin_flip

# Try a more complex one
cargo run --example working_with_distributions

Troubleshooting

Common Issues

Build fails with dependency errors:

# Make sure you're using Rust 1.87+
rustc --version

# Update your dependencies
cargo update

Examples don't run:

# Make sure you're in the project root directory
pwd

# Check example names
ls examples/

IDE doesn't provide completions:

• Make sure rust-analyzer is installed and running
• Try restarting your IDE after installing dependencies
• Check that your Cargo.toml has the correct dependencies

Getting Help

If you encounter issues:

1. Check the GitHub Issues
2. Review the examples for working code
3. Read the API documentation

Next Steps

Installation complete! 🎉

Ready to build your first probabilistic model? → Your First Model

Want to explore examples first? → Complete Tutorials

Time: ~2 minutes • Next: Your First Model

Your First Model

Now that Fugue is installed, let's build your first probabilistic model step by step. We'll start simple and gradually introduce the key concepts.

Step 1: The Simplest Model

Let's start with the simplest possible model - one that always returns the same value:

use fugue::*;

fn constant_model() -> Model<f64> {
    pure(42.0)
}

This model always returns 42.0. The pure function creates a deterministic Model<f64>.

Key insight: Models are descriptions of computations, not the computations themselves.

Step 2: Adding Randomness

Now let's add some randomness by sampling from a probability distribution:

use fugue::*;

fn random_model() -> Model<f64> {
    sample(addr!("x"), Normal::new(0.0, 1.0).unwrap())
}

Step 3: Running Your Model

To actually get values from your model, you need to "run" it with a handler:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let model = sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
    
    // Create a seeded random number generator
    let mut rng = StdRng::seed_from_u64(42);
    
    // Run the model with PriorHandler (forward sampling)
    let (value, trace) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        model,
    );
    
    println!("Sampled value: {:.4}", value);
    println!("Log probability: {:.4}", trace.total_log_weight());
}

Running this outputs:

Sampled value: 1.0175
Log probability: -0.9189

Step 4: Type Safety in Action

Fugue's type safety really shines with discrete distributions:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn type_safe_examples() {
    let mut rng = StdRng::seed_from_u64(42);
    
    // Flip a coin - returns bool directly!
    let (is_heads, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("coin"), Bernoulli::new(0.6).unwrap()),
    );
    
    // Natural boolean usage - no comparisons needed!
    let outcome = if is_heads { "Heads" } else { "Tails" };
    println!("Coin flip: {}", outcome);
    
    // Count events - returns u64 directly!
    let (count, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("events"), Poisson::new(3.0).unwrap()),
    );
    
    println!("Event count: {} (no casting needed!)", count);
    
    // Choose category - returns usize for safe indexing!
    let colors = vec!["red", "green", "blue", "yellow"];
    let (idx, _) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        sample(addr!("color"), Categorical::uniform(4).unwrap()),
    );
    
    println!("Chosen color: {}", colors[idx]); // Safe indexing!
}

# Other PPLs - everything returns float
coin = sample("coin", Bernoulli(0.6))  # Returns 0.0 or 1.0 
if coin == 1.0:  # Need comparison ❌
    ...

count = sample("events", Poisson(3.0))  # Returns float
count_int = int(count)  # Need casting ❌

idx = sample("color", Categorical([...]))  # Returns float  
colors[int(idx)]  # Risky casting and indexing ❌

Step 5: Your First Bayesian Model

Now let's create a simple Bayesian model that learns from data:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn estimate_mean(observation: f64) -> Model<f64> {
    sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap())  // Prior belief
        .bind(move |mu| {
            observe(addr!("y"), Normal::new(mu, 1.0).unwrap(), observation)  // Likelihood
                .map(move |_| mu)  // Return the parameter
        })
}

fn main() {
    let observation = 3.0;  // We observed a value of 3.0
    let model = estimate_mean(observation);
    
    let mut rng = StdRng::seed_from_u64(42);
    let (estimated_mu, trace) = runtime::handler::run(
        runtime::interpreters::PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        model,
    );
    
    println!("Observation: {}", observation);
    println!("Estimated mean: {:.4}", estimated_mu);
    println!("Log probability: {:.4}", trace.total_log_weight());
}

Understanding Model Composition

Fugue models compose using two key operations:

map - Transform Values

let doubled = sample(addr!("x"), Normal::new(0.0, 1.0).unwrap())
    .map(|x| x * 2.0);  // Apply function to the result

bind - Dependent Computations

let dependent = sample(addr!("x"), Normal::new(0.0, 1.0).unwrap())
    .bind(|x| sample(addr!("y"), Normal::new(x, 0.5).unwrap()));  // y depends on x

Key Takeaways

After working through these examples, you should understand:

✅ Models are values: Model<T> represents a probabilistic computation
✅ Safe constructors: Distributions use .new().unwrap() for parameter validation
✅ Type safety: Distributions return natural types (bool, u64, f64)
✅ Addressing: addr!("name") gives names to random variables
✅ Execution: Models need handlers to run and produce values
✅ Composition: Use map and bind to build complex models from simple parts

What's Next?

You can now build and run basic probabilistic models! 🎉

Continue your journey:

Time: ~5 minutes • Next: Understanding Models

Understanding Models

Now that you can build basic models, let's understand the key concepts that make Fugue powerful. This will give you the mental framework to build sophisticated probabilistic programs.

The Big Picture: Models vs Execution

One of Fugue's key insights is separating model specification from execution:

graph LR
    subgraph "Your Code"
        M[Model Definition<br/>What to compute]
    end
    
    subgraph "Runtime System"  
        H[Handler<br/>How to execute]
        T[Trace<br/>What happened]
    end
    
    M --> H
    H --> T

Why this matters: The same model can be executed in different ways:

• Forward sampling (generate data from priors)
• Conditioning (inference given observations)
• Replay (MCMC proposals)
• Scoring (compute probabilities)

Addresses: The Key to Advanced Inference

Every sample and observe site needs a unique address:

use fugue::*;

// Good addressing ✅
let model = sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
    .bind(|mu| sample(addr!("sigma"), LogNormal::new(0.0, 0.5).unwrap())
        .bind(move |sigma| {
            observe(addr!("y1"), Normal::new(mu, sigma).unwrap(), 2.1);
            observe(addr!("y2"), Normal::new(mu, sigma).unwrap(), 1.9);
            pure((mu, sigma))
        }));

Addressing Patterns

Simple names for scalar parameters:

let mu = sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());

Indexed addresses for collections:

for i in 0..10 {
    let x_i = sample(addr!("x", i), Normal::new(mu, 1.0).unwrap());
}

Scoped addresses for hierarchical models:

let encoder_z = sample(scoped_addr!("encoder", "z"), dist);
let decoder_z = sample(scoped_addr!("decoder", "z"), dist);

// ❌ DON'T: Random or non-deterministic addresses  
let addr = format!("param_{}", rng.gen::<u32>());  // NEVER!

// ❌ DON'T: Reuse addresses for different purposes
sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
sample(addr!("x"), Bernoulli::new(0.5).unwrap());  // Collision!

// ❌ DON'T: Missing addresses in loops
for i in 0..10 {
    sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());  // All same address!
}

Model Composition: The Monadic Structure

Fugue models follow a monadic pattern that makes complex models composable:

Three Fundamental Operations

graph TB
    subgraph "Model Building Blocks"
        P[pure: A → Model⟨A⟩<br/>Lift values into models]
        M[map: Model⟨A⟩ → ⟨A → B⟩ → Model⟨B⟩<br/>Transform outputs]  
        B[bind: Model⟨A⟩ → ⟨A → Model⟨B⟩⟩ → Model⟨B⟩<br/>Chain dependent computations]
    end

pure - Lift Values

Use pure to inject deterministic values into the probabilistic context:

use fugue::*;

// Lift a constant
let constant = pure(42.0);

// Lift computed values
let computed = pure(data.iter().sum::<f64>() / data.len() as f64);

map - Transform Outputs

Use map when you want to transform the result without adding randomness:

use fugue::*;

// Transform a single sample
let squared = sample(addr!("x"), Normal::new(0.0, 1.0).unwrap())
    .map(|x| x * x);

// Combine multiple values  
let sum = zip(
    sample(addr!("a"), Normal::new(0.0, 1.0).unwrap()),
    sample(addr!("b"), Normal::new(0.0, 1.0).unwrap())
).map(|(a, b)| a + b);

bind - Chain Dependent Computations

Use bind when the next random choice depends on a previous one:

use fugue::*;

// Dependent sampling: variance depends on mean
let hierarchical = sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
    .bind(|mu| sample(addr!("sigma"), LogNormal::new(mu.abs().ln(), 0.1).unwrap())
        .bind(move |sigma| sample(addr!("y"), Normal::new(mu, sigma).unwrap())));

// Conditional branching: choice affects distribution  
let mixture = sample(addr!("component"), Bernoulli::new(0.5).unwrap())
    .bind(|component| {
        if component {
            sample(addr!("value"), Normal::new(-2.0, 1.0).unwrap())
        } else {
            sample(addr!("value"), Normal::new(2.0, 1.0).unwrap())  
        }
    });

Advanced Composition Patterns

Building Collections

use fugue::*;

// Fixed-size collection
let samples = traverse_vec((0..10).collect(), |i| {
    sample(addr!("x", i), Normal::new(0.0, 1.0).unwrap())
});

// Data-driven collection
let observations = traverse_vec(data, |datum| {
    observe(addr!("y", datum.id), Normal::new(datum.mu, 1.0).unwrap(), datum.value)
});

Conditional Models

use fugue::*;

fn model_selection(data: &[f64]) -> Model<String> {
    sample(addr!("use_robust"), Bernoulli::new(0.2).unwrap())
        .bind(|use_robust| {
            if use_robust {
                // Robust model with t-distribution
                sample(addr!("df"), LogNormal::new(1.0, 0.5).unwrap())
                    .bind(|df| {
                        // Hypothetical t-distribution sampling
                        pure("robust".to_string())
                    })
            } else {
                // Standard normal model
                sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
                    .map(|_| "normal".to_string())
            }
        })
}

Hierarchical Structure

use fugue::*;

fn hierarchical_model(groups: Vec<Vec<f64>>) -> Model<Vec<f64>> {
    // Global hyperparameters
    sample(addr!("global_mu"), Normal::new(0.0, 2.0).unwrap())
        .bind(|global_mu| {
            sample(addr!("global_sigma"), LogNormal::new(0.0, 0.5).unwrap())
                .bind(move |global_sigma| {
                    // Group-level parameters  
                    traverse_vec(groups.into_iter().enumerate().collect(), move |(g, data)| {
                        sample(addr!("group_mu", g), Normal::new(global_mu, global_sigma).unwrap())
                            .bind(move |group_mu| {
                                // Individual observations
                                traverse_vec(data.into_iter().enumerate().collect(), move |(i, y)| {
                                    observe(addr!("y", g, i), Normal::new(group_mu, 1.0).unwrap(), y)
                                }).map(move |_| group_mu)
                            })
                    })
                })
        })
}

Mental Models for Success

Think in Terms of Generative Stories

Ask yourself: "How could this data have been generated?"

// Story: "Each person has a skill level, and their performance 
//         on each task reflects that skill plus task-specific noise"

let model = sample(addr!("skill"), Normal::new(100.0, 15.0).unwrap())  // Person's skill
    .bind(|skill| {
        traverse_vec(tasks, move |task| {
            let difficulty = task.difficulty;
            let expected_score = skill - difficulty;
            observe(addr!("score", task.id), Normal::new(expected_score, 5.0).unwrap(), task.actual_score)
        }).map(move |_| skill)
    });

Separate What from How

• What: Model describes relationships and distributions
• How: Handler determines execution strategy (sampling, inference, etc.)

// The SAME model...
let model = build_regression_model(&data);

// Can be executed different ways:
let (sample, _) = run(PriorHandler { /*...*/ }, model.clone());      // Forward sampling
let (_, scored_trace) = run(ScoreGivenTrace { /*...*/ }, model.clone()); // Compute likelihood  
let mcmc_chain = adaptive_mcmc_chain(rng, || model.clone(), 1000, 500); // MCMC inference

Key Takeaways

You now understand the foundational concepts:

✅ Separation of Concerns: Models describe computations, handlers execute them
✅ Addressing Strategy: Unique, stable addresses enable advanced inference
✅ Monadic Composition: pure, map, bind build complex models from simple parts
✅ Compositional Patterns: Collections, conditionals, and hierarchical structures
✅ Generative Thinking: Model the data generation process

What's Next?

You have the conceptual foundation to build sophisticated models! 🎉

Time: ~8 minutes • Next: Running Inference

Running Inference

You now know how to build probabilistic models. But models alone don't give you answers - you need inference to extract insights from them. Let's explore Fugue's inference algorithms!

What is Inference?

Inference is the process of learning about model parameters after seeing data. In Bayesian terms:

$$\text{Posterior}=\frac{\text{Prior}\times \text{Likelihood}}{\text{Evidence}}$$

graph LR
    subgraph "Before Data"
        P["Prior Beliefs<br/>p(theta)"]
    end

    subgraph "Observing Data"
        L["Likelihood<br/>p(y|theta)"]
        D["Data<br/>y₁, y₂, ..."]
    end

    subgraph "After Data"
        Post["Posterior Beliefs<br/>p(theta|y)"]
    end

    P --> Post
    L --> Post
    D --> Post

The Challenge

Most real models don't have analytical solutions. We need algorithms to approximate the posterior distribution.

Fugue's Inference Arsenal

1. MCMC (Markov Chain Monte Carlo) 🥇

Best for: Most general-purpose Bayesian inference

How it works: Generates samples that approximate the posterior distribution

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn coin_bias_model(heads: u64, total: u64) -> Model<f64> {
    sample(addr!("bias"), Beta::new(1.0, 1.0).unwrap())  // Prior
        .bind(move |bias| {
            observe(addr!("heads"), Binomial::new(total, bias).unwrap(), heads)  // Likelihood
                .map(move |_| bias)
        })
}

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // Run adaptive MCMC
    let samples = inference::mh::adaptive_mcmc_chain(
        &mut rng,
        || coin_bias_model(7, 10),  // 7 heads out of 10 flips
        1000,  // number of samples
        500,   // warmup samples
    );

    // Extract bias estimates
    let bias_samples: Vec<f64> = samples.iter()
        .filter_map(|(_, trace)| trace.get_f64(&addr!("bias")))
        .collect();

    let mean_bias = bias_samples.iter().sum::<f64>() / bias_samples.len() as f64;
    println!("Estimated bias: {:.3}", mean_bias);

    // Compute 90% credible interval
    let mut sorted = bias_samples.clone();
    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let lower = sorted[(0.05 * sorted.len() as f64) as usize];
    let upper = sorted[(0.95 * sorted.len() as f64) as usize];
    println!("90% credible interval: [{:.3}, {:.3}]", lower, upper);
}

When to use MCMC:

• ✅ Want exact posterior samples
• ✅ Moderate number of parameters (< 100)
• ✅ Can afford computation time
• ✅ Model evaluation is reasonably fast

2. HMC (Hamiltonian Monte Carlo) 🎢

Best for: Continuous, correlated, or higher-dimensional posteriors where single-site MCMC mixes slowly

How it works: Treats the negative log-posterior as a landscape and rolls a simulated ball across it using its gradient, moving every continuous parameter together instead of one at a time

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // Same conjugate model, but HMC moves `bias` using gradient information
    // instead of a random-walk proposal.
    let samples = inference::hmc::hmc_chain(
        &mut rng,
        || coin_bias_model(7, 10),
        1000,               // number of samples
        500,                // warmup iterations (step-size adaptation)
        HMCConfig::default(), // 16 leapfrog steps, target 80% acceptance
    );

    let bias_samples: Vec<f64> = samples.iter()
        .filter_map(|(_, trace)| trace.get_f64(&addr!("bias")))
        .collect();

    let mean_bias = bias_samples.iter().sum::<f64>() / bias_samples.len() as f64;
    println!("HMC estimated bias: {:.3}", mean_bias);
}

Fugue's HMC computes the gradient with central finite differences (models are plain Rust closures, not auto-diff traces), then uses a leapfrog integrator and an exact Metropolis correction against the true log-density — the finite-difference force only affects efficiency, never correctness. Step size is tuned automatically during warmup via dual averaging toward an 80% target acceptance rate (Hoffman & Gelman 2014); discrete sites are held fixed for the duration of the HMC update (Metropolis-within-Gibbs), so compose with adaptive_mcmc_chain when a model mixes continuous and discrete latents.

When to use HMC:

• ✅ All (or most) latent variables are continuous
• ✅ Parameters are correlated (HMC exploits gradient direction; single-site MH can't)
• ✅ You want fewer iterations to reach the same effective sample size
• ⚠️ Purely discrete models get no benefit — HMC degenerates to prior draws when there are no continuous sites

3. SMC (Sequential Monte Carlo) 🎯

Best for: Sequential data and online learning

How it works: Uses particles to approximate the posterior, good for streaming data

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // Generate particles from prior
    let particles = inference::smc::smc_prior_particles(
        &mut rng,
        1000,  // number of particles
        || coin_bias_model(7, 10),
    );

    println!("Generated {} particles", particles.len());

    // Compute weighted posterior mean
    let total_weight: f64 = particles.iter().map(|p| p.weight).sum();
    let weighted_mean: f64 = particles.iter()
        .filter_map(|p| {
            p.trace.get_f64(&addr!("bias"))
                .map(|bias| bias * p.weight)
        })
        .sum::<f64>() / total_weight;

    println!("Weighted posterior mean: {:.3}", weighted_mean);

    // Check effective sample size
    let weights: Vec<f64> = particles.iter().map(|p| p.weight).collect();
    let ess = 1.0 / weights.iter().map(|w| w * w).sum::<f64>();
    println!("Effective sample size: {:.1}", ess);
}

When to use SMC:

• ✅ Sequential/streaming data
• ✅ Online inference needed
• ✅ Many discrete latent variables
• ✅ Want to visualize inference process

4. Variational Inference (VI) ⚡

Best for: Fast approximate inference with many parameters

How it works: Finds the best approximation within a family of simple distributions

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // Estimate ELBO (Evidence Lower BOund)
    let elbo = inference::vi::estimate_elbo(
        &mut rng,
        || coin_bias_model(7, 10),
        100,  // number of samples for estimation
    );

    println!("ELBO estimate: {:.3}", elbo);

    // For more sophisticated VI, you'd set up a variational guide
    // and optimize it (see the VI tutorial for details)
}

When to use VI:

• ✅ Need fast approximate inference
• ✅ Many parameters (> 100)
• ✅ Can accept approximation error
• ✅ Want predictable runtime

5. ABC (Approximate Bayesian Computation) 🎲

Best for: Models where likelihood is intractable or expensive

How it works: Simulation-based inference using distance between simulated and observed data

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let mut rng = StdRng::seed_from_u64(42);

    // ABC with summary statistics
    let observed_summary = 0.7; // 7/10 = 0.7 success rate
    let samples = inference::abc::abc_scalar_summary(
        &mut rng,
        || sample(addr!("bias"), Beta::new(1.0, 1.0).unwrap()), // Prior only
        |trace| trace.get_f64(&addr!("bias")).unwrap_or(0.0), // Extract bias
        observed_summary,  // Target summary statistic
        0.1,              // Tolerance
        1000,             // Max samples to try
    );

    println!("ABC accepted {} samples", samples.len());

    if !samples.is_empty() {
        let abc_estimates: Vec<f64> = samples.iter()
            .filter_map(|trace| trace.get_f64(&addr!("bias")))
            .collect();
        let abc_mean = abc_estimates.iter().sum::<f64>() / abc_estimates.len() as f64;
        println!("ABC estimated bias: {:.3}", abc_mean);
    }
}

When to use ABC:

• ✅ Likelihood is intractable or very expensive
• ✅ Can simulate from the model easily
• ✅ Have good summary statistics
• ✅ Can tolerate approximation error

Algorithm Comparison

Practical Inference Workflow

Here's a typical workflow for real inference:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn inference_workflow() {
    let mut rng = StdRng::seed_from_u64(42);

    // 1. Define your model
    let model = || coin_bias_model(17, 25);  // 17 heads out of 25 flips

    // 2. Run inference (adaptive MCMC is often a good default)
    let samples = inference::mh::adaptive_mcmc_chain(
        &mut rng,
        model,
        2000,  // samples
        1000,  // warmup
    );

    // 3. Extract parameter values
    let bias_samples: Vec<f64> = samples.iter()
        .filter_map(|(_, trace)| trace.get_f64(&addr!("bias")))
        .collect();

    // 4. Compute summary statistics
    let mean = bias_samples.iter().sum::<f64>() / bias_samples.len() as f64;
    let variance = bias_samples.iter()
        .map(|&x| (x - mean).powi(2))
        .sum::<f64>() / (bias_samples.len() - 1) as f64;
    let std_dev = variance.sqrt();

    println!("Posterior Summary:");
    println!("  Mean: {:.3}", mean);
    println!("  Std Dev: {:.3}", std_dev);

    // 5. Check convergence (autocorrelation-based effective sample size)
    let ess = inference::diagnostics::effective_sample_size(&bias_samples);
    println!("  Effective Sample Size: {:.1}", ess);
    // For multiple chains, also check split-R̂ (inference::diagnostics::r_hat_f64)
    // and inference::mcmc_utils::effective_sample_size_multichain.

    if ess > 100.0 {
        println!("  ✅ Good mixing!");
    } else {
        println!("  ⚠️ Poor mixing - consider more samples");
    }

    // 6. Make predictions
    println!("\nPredictions:");
    println!("  P(bias > 0.5) = {:.2}",
        bias_samples.iter().filter(|&&b| b > 0.5).count() as f64 / bias_samples.len() as f64);
}

Choosing the Right Algorithm

Decision Tree

graph TD
    A[Need inference?] -->|Yes| B[Real-time/online?]
    A -->|No| Z[Use PriorHandler<br/>for forward sampling]

    B -->|Yes| SMC[SMC]
    B -->|No| C[Likelihood tractable?]

    C -->|No| ABC[ABC]
    C -->|Yes| D[Many parameters?]

    D -->|Yes > 100| VI[Variational Inference]
    D -->|No < 100| E[Need exact samples?]

    E -->|Yes| F[Mostly continuous & correlated?]
    E -->|No| VI2[VI for speed]

    F -->|Yes| HMC[HMC]
    F -->|No, mostly discrete| MCMC[MCMC]

Rules of Thumb

1. Start with MCMC for most problems - it's the most general
2. Reach for HMC when parameters are continuous and correlated - it mixes far faster than single-site MCMC by using the gradient
3. Use SMC if you have sequential/streaming data
4. Use VI if you need speed and can accept approximation
5. Use ABC only when likelihood is truly intractable

Key Takeaways

You now know how to extract insights from your models:

✅ Inference Purpose: Learn parameters from data using Bayesian updating
✅ Algorithm Options: MCMC, HMC, SMC, VI, ABC each have their strengths
✅ Practical Workflow: Define model → Run inference → Extract parameters → Check diagnostics
✅ Algorithm Selection: Choose based on problem characteristics and requirements

What's Next?

You've completed Getting Started! 🎉

Time: ~5 minutes • Next: Complete Tutorials

Explorables

Most explanations of probabilistic programming ask you to read first and understand later. These pages invert that. Every figure here is a small machine: drag the numbers, click the data, step the algorithms, and watch the mathematics respond. You will have felt why a sampler stalls or a particle filter collapses before you meet the equation that says so. The genre is borrowed with gratitude from Bret Victor's explorable explanations and 3Blue1Brown's visual mathematics.

Two habits worth forming:

• Touch everything. If a number looks interesting, try dragging it — dashed underlines mark the ones that respond. Canvases are clickable more often than not.
• Mind the seed. Wherever randomness appears, a scrubbable seed appears with it. The same seed always replays the same run — which is not a gimmick, it is fugue's worldview: a recorded run is a trace, and traces can be replayed. You will meet this idea everywhere.

One color language runs through every page: prior × likelihood = posterior, with coral for the current sample and violet for momentum and structure. Once you know it, every canvas and every equation on this site reads at a glance.

The six machines

1. Anatomy of a Probabilistic Program — a coin, a prior you can bend, and data you can click into existence. Bayes' rule as something your hands learn first.
2. The Model Is a Score — five observations you can drag along an axis, a prior, and an exact posterior that follows your hand. Underneath, fugue's Model monad performed one effect at a time — then handed to a different performer. The page that explains fugue itself.
3. Random Walks in Posterior Space — a real regression: drag the data points and watch the posterior heatmap deform under the sampler's feet, while accepted samples paint fit-lines through your data. Live split-R̂ and ESS referee the whole thing.
4. Rolling, Not Guessing: Hamiltonian Monte Carlo — the same regression, but the sample gets momentum and rolls. Leapfrog trajectories, divergences, and a side-by-side race against the random walk at a matched budget — visible in the data, not just the parameters.
5. Particles That Tell Stories — Sequential Monte Carlo as a population: propagate, weight, resample. Watch lineages go extinct and understand degeneracy by witnessing it.
6. A Field Guide to Distributions — every distribution fugue ships, with real parameterizations, a sampler racing its own density, and the natural return types that make fugue's models type-safe.

Read them in order for a course, or jump to what you came for — each page stands alone and links onward. When a page convinces you, the matching tutorial turns the intuition into working Rust; the code on every explorable compiles against fugue-ppl 0.2.0 as written.

One more honesty note, and a point of pride: the samplers behind the MH, HMC, and SMC machines are not JavaScript imitations — they are the fugue crate itself, compiled to WebAssembly, driven one transition per animation frame. Same kernels, same diagnostics, same seeded reproducibility as production; only the drawing is JavaScript. (Browsers without the WASM package fall back to a faithful JS mirror.) When you want to type your own model at these kernels, the Playground is exactly that.

Anatomy of a Probabilistic Program

A probabilistic program has three moving parts: a prior (what you believe before data), a likelihood (how data is explained), and a posterior (what you believe after). Here they are, side by side, for the oldest question in statistics: is this coin fair? Flip a chip, drag the prior, and watch the green curve reshape itself in real time — each change leaving a fading ghost of the belief it just replaced.

You start with a prior belief about the coin's bias p: Beta(α = 2, β = 2). The blue curve is that belief. Each coin chip is one observation; click any chip to flip it between heads and tails. The posterior — the green curve — is the prior times the likelihood, renormalized. Press Replay to watch the flips arrive one at a time and the green curve walk from the prior to the posterior — Bayesian updating, animated — or press Step to take that walk one flip at a time yourself. Press Deal to draw a fresh dozen flips from seed 11; the same seed always deals the same coins, because in fugue a seeded run is a replayable recording.

Things to try

1. Add tails until the green curve's peak crosses left of p = 0.5 — that is the moment the data outvotes a "fair coin" prior.
2. Drag β up to 20 with only a few flips: a stubborn prior barely moves. Now add twenty flips — data eventually wins, no matter the prior.
3. Set α = 0.5, β = 0.5 (the Jeffreys prior). The blue curve bends up at both ends: it says "this coin is probably rigged one way or the other."
4. Turn show likelihood on and watch the yellow curve. The green posterior always sits between blue and yellow — pulled toward whichever is sharper.
5. Press Replay: the posterior starts as the blue prior, then each flip nudges it — heads pull right, tails pull left. The fading green ghosts are the beliefs it held along the way. This is Bayesian updating, one datum at a time.
6. Scrub the seed while Dealing: every value gives a different but fully reproducible dataset. That reproducibility is the whole point of a trace.

What you just saw

The green curve is not drawn by a formula you have to trust — it is the literal product of the other two, normalized so its area is one:

$$p(p\mid \mathcal{D})\propto p(p)\times p(\mathcal{D}\mid p)$$

For a coin this product has a closed form. A Beta prior multiplied by h heads and t tails of Bernoulli likelihood is again a Beta — the two are conjugate — so the posterior is exact:

$$\underset{\text{prior}}{\underbrace{\mathrm{Beta}(\alpha ,\beta )}}\times \underset{\text{likelihood}}{\underbrace{p^h(1-p{)}^t}}\propto \underset{\text{posterior}}{\underbrace{\mathrm{Beta}(\alpha +h,\beta +t)}}$$

Every heads slides one unit of belief into α; every tails, into β. That is why Replay works: feeding the flips one at a time, using each posterior as the next prior, lands on exactly the same green curve as folding them in all at once. Bayesian updating is associative. The posterior mean is the updated ratio, and the readouts track it live:

$$\mathbb{E}[p\mid \mathcal{D}]=\frac{\alpha +h}{\alpha +\beta +h+t}$$

The coral marker is the MAP — the posterior's most probable bias, its mode. The shaded band is the 90% credible interval: the model's honest "I'm 90% sure the bias is in here." Conjugacy is a lucky gift of the coin; most models have no such shortcut. That is why fugue exists.

The fugue code

The widget uses conjugacy because it can. fugue does not need to — it runs Metropolis–Hastings on the same model and lands on the same number.

use fugue::*;
use fugue::inference::mh::adaptive_mcmc_chain;
use rand::{SeedableRng, rngs::StdRng};

// The score: a prior over the bias, then one observe per flip.
fn coin(data: Vec<bool>) -> Model<f64> {
    prob!(
        // PRIOR — the blue curve: belief about the bias p before any flip.
        let p <- sample(addr!("p"), Beta::new(2.0, 2.0).unwrap());

        // LIKELIHOOD — each yellow chip is one `observe`: a flip explained by p.
        let _obs <- plate!(i in 0..data.len() => {
            observe(addr!("flip", i), Bernoulli::new(p).unwrap(), data[i])
        });

        pure(p) // return the inferred bias
    )
}

fn main() {
    // 7 heads, 3 tails — the widget's starting data.
    let data = vec![true, false, true, true, false, true, true, false, true, true];
    let mut rng = StdRng::seed_from_u64(11);

    // Metropolis–Hastings — no conjugacy assumed, adaptive step size.
    let samples = adaptive_mcmc_chain(&mut rng, || coin(data.clone()), 4000, 1000);
    let ps: Vec<f64> = samples
        .iter()
        .filter_map(|(_, t)| t.get_f64(&addr!("p")))
        .collect();
    let mean = ps.iter().sum::<f64>() / ps.len() as f64;

    // Conjugacy says the posterior is Beta(2+7, 2+3) = Beta(9, 5),
    // whose mean is 9/14 ≈ 0.643. MH agrees without ever knowing that.
    println!("MH posterior mean ≈ {mean:.3}  (analytic 9/14 = {:.3})", 9.0 / 14.0);
}

sample records a choice at an address; observe scores the data against the current bias and folds a log-weight into the trace; pure returns a value. The prob! and plate! macros are sugar over the Model monad — the score is written once, then performed by a handler. That separation is the next explorable.

Go deeper

• Tutorial: Bayesian Coin Flip walks the same model without the pictures.
• Next explorable: The Model Is a Score — step through the interpreter that actually runs this program.
• API: Beta, adaptive_mcmc_chain.

Next: The Model Is a Score

The Model Is a Score

A Model in fugue is not a running program. It is a score — pure notation that describes what could happen, note by note, but makes no sound on its own. A handler is the performer that reads the score and decides what each note means: improvise a fresh value, or replay one from a recording.

Below is a real inference problem. Five data points say where the world landed; a prior says what you believed about the mean μ beforehand; the posterior is the answer — updated live. Grab a yellow dot and drag it. The green curve follows in real time. The strip below the picture is the machinery fugue actually steps through to get there.

This run uses seed 11. A seeded run is a replayable trace — the same seed always draws the same μ.

Things to try

1. Drag any yellow dot to the right. The green posterior slides after it and the posterior mean readout climbs — the data pulled your belief.
2. Press Step five times. Watch the SampleF64 chip fire first (it sets μ), then each ObserveF64 chip light up and its data dot pulse in the picture — chip and dot are the same event, seen twice.
3. Press Perform ×200 with the PriorHandler active: 200 fresh μ draws rain down as a blue cloud that traces the prior curve. The handler is improvising.
4. Flip on the Replay handler and press Perform ×200 again: all 200 draws stack on one value — a single coral spike. The handler is performing a fixed recording, not drawing.
5. Scrub the seed in prior mode — the coral μ tick jumps to a new draw. Scrub it in replay mode — it does not move. Same score, different performer.

What you just saw

The score is a tiny probabilistic program: sample a mean, then observe five data points drawn from it.

$$\mu \sim \mathrm{Normal}(0,2),y_i\sim \mathrm{Normal}(\mu ,1)i=1,\dots ,5.$$

Because both the prior and the likelihood are Gaussian, the posterior over μ is Gaussian too, in closed form — this is the exact green curve you were dragging:

$$\mu \mid y\sim \mathrm{Normal}({\mu }_n,{\sigma }_n^2),{\sigma }_n^2={\left(\frac{1}{2^2}+\frac{n}{1^2}\right)}^{-1},{\mu }_n={\sigma }_n^2\left(\frac{0}{2^2}+\frac{{\sum }_i{y}_i}{1^2}\right).$$

With the five default points (∑yᵢ = 6.0), that is Normal(1.143, 0.436²) — the numbers in the readout. The score never computes this. The handler does, one node at a time, accumulating a log-weight:

$$\underset{\text{total log-weight}}{\underbrace{\log p(\mu ,y)}}=\underset{\text{log\_prior}}{\underbrace{\log p(\mu )}}+\underset{\text{log\_likelihood}}{\underbrace{{\sum }_i\log p(y_i\mid \mu )}}.$$

The key idea: effects are interpreted, not performed. The score names a sample site; it does not draw. Whether a value is improvised (PriorHandler) or replayed (ReplayHandler) is the handler's decision, made later. That is what the Perform ×200 button shows: the same score becomes a cloud of guesses or a single fixed spike, depending only on who performs it. This separation is the machinery every inference algorithm in fugue is built from — and fugue's MCMC reaches that same green curve without knowing the conjugate formula, by walking this chain over and over.

Monads, demystified in one paragraph

A Model is built from two operations. bind is "and then": run this node, then feed its value to a function that produces the next node. pure is "done": wrap a plain value as a finished model. That is the entire monad — and then and done. Each node stores its continuation k, the rest of the score. The score is a linked list of "and then"s ending in a "done", and it does not play itself.

The fugue code

Here is the exact score from the widget. The sample line is the SampleF64 chip; each observe in the fold is an ObserveF64 chip lighting its yellow dot; .map(move |_| mu) is the trailing Pure(μ).

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

// The score: sample the mean, then observe five data points given it.
fn model(data: Vec<f64>) -> Model<f64> {
    sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap()).bind(move |mu| {
        // Fold the observations into the chain: each one is an ObserveF64 node.
        let mut m = pure(mu);
        for (i, y) in data.into_iter().enumerate() {
            m = m.bind(move |mu| {
                observe(addr!("y", i), Normal::new(mu, 1.0).unwrap(), y)
                    .map(move |_| mu)
            });
        }
        m // ends in Pure(mu): the score terminates, returning mu
    })
}

fn main() {
    let data = vec![1.3, 0.7, 2.1, 0.4, 1.5];

    // PriorHandler improvises: it draws a fresh mu from the RNG.
    let mut rng = StdRng::seed_from_u64(11);
    let (mu, recording) = runtime::handler::run(
        PriorHandler { rng: &mut rng, trace: Trace::default() },
        model(data.clone()),
    );
    println!("prior draw: mu = {mu:.3}, total log-weight = {:.3}",
             recording.total_log_weight());

    // ReplayHandler performs that recording: the same mu, no fresh draw.
    let (mu_again, _) = runtime::handler::run(
        ReplayHandler { rng: &mut rng, base: recording, trace: Trace::default() },
        model(data),
    );
    assert_eq!(mu, mu_again); // same score, replayed exactly
}

Under the hood, the score is the Model enum. Every effect variant carries its address, its distribution, and a boxed continuation k — the rest of the score (abridged from src/core/model.rs):

pub enum Model<A> {
    Pure(A),
    SampleF64 {
        addr: Address,
        dist: Box<dyn Distribution<f64>>,
        k: Box<dyn FnOnce(f64) -> Model<A> + Send + 'static>, // "and then"
    },
    ObserveF64 {
        addr: Address,
        dist: Box<dyn Distribution<f64>>,
        value: f64,
        k: Box<dyn FnOnce(()) -> Model<A> + Send + 'static>,
    },
    Factor { logw: f64, k: Box<dyn FnOnce(()) -> Model<A> + Send + 'static> },
    // ... plus SampleBool/U64/Usize/I64 and their Observe twins, one per
    // return type — this is fugue's type-safety story: a Bernoulli site yields
    // a bool, a Poisson site a u64, never an untyped float.
}

The handler you toggled is any type that answers those effects. run is the interpreter — a flat trampoline, not recursion (abridged from src/runtime/handler.rs):

pub fn run<A>(mut h: impl Handler, m: Model<A>) -> (A, Trace) {
    let mut m = m;
    let a = loop {
        m = match m {
            Model::Pure(a) => break a,                    // "done": return the value
            Model::SampleF64 { addr, dist, k } => {
                let x = h.on_sample_f64(&addr, &*dist);    // ask the handler
                k(x)                                       // advance to the rest
            }
            Model::ObserveF64 { addr, dist, value, k } => {
                h.on_observe_f64(&addr, &*dist, value);    // score the data
                k(())
            }
            Model::Factor { logw, k } => { h.on_factor(logw); k(()) }
            // ... one arm per variant ...
        };
    };
    (a, h.finish())
}

Each node hands back its continuation k(value) and the loop goes around again. Nothing recurses, so the interpreter runs in constant stack depth — a model with 100 000 sample sites is interpreted without overflowing the stack.

Go deeper

• Understanding Models — the conceptual tour of Model, bind, and pure.
• Custom Handlers — write your own performer.
• Trace Manipulation — read and edit the recordings a handler produces.
• Model and Handler on docs.rs.
• Next explorable: Random Walks in Posterior Space — now that the score has a value and a weight, how do we listen our way back to the posterior you just watched form?

Next: Random Walks in Posterior Space

Random Walks in Posterior Space

You have data and a line you want to fit through it — but you want the whole posterior over slope and intercept, not a single answer. That posterior has no formula you can read off. So you walk. Stand at some (slope, intercept), propose a small random step, and keep it more often when it explains the data better. Do this long enough and the places you linger are the posterior.

That walk is Metropolis-Hastings. The score names a shape — here, the posterior over two numbers; the walker is a performer who never sees the whole shape at once, only whether the next fit sounds better or worse than the last.

Two linked spaces. On the left, the data — twelve points you can grab and drag. Every green line is a recently accepted fit; the coral line is where each chain stands now; rejected proposals flash coral-dashed and vanish. On the right, the same walk seen in parameter space: a live posterior heatmap over (slope, intercept), with the coral dot threading it and blue proposal arrows firing each tick. A point in the right panel is a line in the left panel — watch them move together.

The model is honest Bayesian linear regression: $y\sim \mathcal{N}(ax+b,{\sigma }_{\text{obs}})$ with ${\sigma }_{\text{obs}}$ fixed at $0.8$, and priors $a,b\sim \mathcal{N}(0,2.5)$. The split-R̂ and ESS readouts are fugue's real convergence diagnostics, computed live on the samples accruing on screen.

Things to try

1. The chains are already walking — drag a point far off the line. The right-hand heatmap morphs and the whole chain migrates to the new best fit — same frame. This is the moment: the data is the posterior, and you are reshaping it with your cursor.
2. Drag PROPOSAL σ down to 0.02. Acceptance climbs toward 100% — every tiny step is safe — yet the coral dot barely crawls and R̂ stays stubbornly above 1. High acceptance is not the goal.
3. Now drag σ up to 4. Steps overshoot the tight posterior ridge and almost every proposal is rejected. The chain freezes, twitching in place. Too bold is as stuck as too timid.
4. Find the Goldilocks band (σ near 0.3–0.6). Acceptance settles around 25–45%, the green spaghetti fans tightly around the true line, and R̂ falls toward 1.0.
5. Set CHAINS to 4. They start from dispersed corners of parameter space. When σ is small, watch them stay marooned apart — R̂ stays high because they disagree. One chain alone could never have told you it was stuck.

What you just saw

Metropolis-Hastings builds a Markov chain whose stationary distribution is the posterior. Write the parameters as $\theta =(a,b)$. From the current state you draw a proposal ${\theta }^{\prime }$ from a symmetric Gaussian kernel, ${\theta }^{\prime }\sim \mathcal{N}(\theta ,{\sigma }^{2}I)$, and accept it with probability

$$\alpha =\min \left(1,\frac{p({\theta }^{\prime })p(\mathcal{D}\mid {\theta }^{\prime })}{p(\theta )p(\mathcal{D}\mid \theta )}\right).$$

The ratio is prior times likelihood, new over old — the evidence cancels, which is the whole trick: you never need the intractable normalizer. Because the proposal is symmetric, $q({\theta }^{\prime }\mid \theta )=q(\theta \mid {\theta }^{\prime })$ drops out too. Here the prior is $\mathcal{N}(0,2.5)$ on each of $a,b$ and the likelihood is the product of the twelve $\mathcal{N}(a{x}_i+b,0.8)$ terms — one per yellow point. Computed in log space, where the widget lives, acceptance is a subtraction:

$$\log \alpha =\min (0,\log p({\theta }^{\prime })-\log p(\theta )).$$

Uphill moves ($\log \alpha =0$) are always taken; downhill moves are taken with probability $e^{\log \alpha }$. That occasional downhill step is what lets the chain explore the full spread of plausible lines instead of collapsing onto the single best fit — which is exactly the green spaghetti you see fanning around the data.

Why σ is a dial, not a detail. The proposal scale trades off two failures. Too small and consecutive samples are nearly identical — high autocorrelation, so your 2000 draws carry the information of a handful. Too large and you reject constantly, so the chain sits still — again few effective draws. The effective sample size (ESS) measures exactly this: how many independent draws your correlated chain is worth.

$$\mathrm{ESS}=\frac{mn}{\hat{\tau }},\hat{\tau }=1+2\sum _{k\ge 1}{\hat{\rho }}_k,$$

where ${\hat{\rho }}_k$ is the autocorrelation at lag $k$ and $\hat{\tau }$ the integrated autocorrelation time. Fugue estimates $\hat{\tau }$ with Geyer's initial positive sequence, pooled across chains (Vehtari et al. 2021).

Why more than one chain. A single walker stuck in one corner of parameter space looks perfectly converged from the inside. Split-R̂ compares the variance between chains to the variance within them, after splitting each chain in half so a slow within-chain drift can't hide:

$$\hat{R}=\sqrt{\frac{{\widehat{\mathrm{var}}}^{+}}{W}},{\widehat{\mathrm{var}}}^{+}=\frac{n-1}{n}W+\frac{1}{n}B.$$

$W$ is the within-chain variance, $B$ the between-chain variance. When the chains agree, $B\to 0$ and $\hat{R}\to 1$. Anything above $1.01$ means they haven't mixed. This is the same split statistic fugue reports from r_hat_f64 — the widget ports its arithmetic verbatim.

The fugue code

The widget tunes σ by hand. Fugue's adaptive_mcmc_chain does it for you, nudging each site's proposal scale toward the 0.44 acceptance rate you just discovered is healthy — then you referee convergence with the very diagnostics on screen. This is the same regression model as the widget: slope and intercept sampled from priors, one observe per data point.

use fugue::*;
use fugue::inference::mh::adaptive_mcmc_chain;
use fugue::inference::diagnostics::r_hat_f64;
use fugue::inference::mcmc_utils::effective_sample_size_mcmc;
use rand::{rngs::StdRng, SeedableRng};

// Bayesian linear regression — the same model the widget samples.
// slope (a) and intercept (b) each get a Normal(0, 2.5) prior; every data
// point contributes one yellow observation y_i ~ Normal(a*x_i + b, 0.8).
fn regression(x: Vec<f64>, y: Vec<f64>) -> impl Fn() -> Model<(f64, f64)> {
    move || {
        let x = x.clone();
        let y = y.clone();
        prob! {
            let a <- sample(addr!("slope"),     Normal::new(0.0, 2.5).unwrap());
            let b <- sample(addr!("intercept"), Normal::new(0.0, 2.5).unwrap());
            // One observe per point — this loop IS the yellow data panel.
            let _obs <- plate!(i in 0..x.len() => {
                observe(addr!("y", i), Normal::new(a * x[i] + b, 0.8).unwrap(), y[i])
            });
            pure((a, b))
        }
    }
}

fn main() {
    // Synthetic data from the true line y = 0.8x - 0.4 + noise, same as the widget.
    let x: Vec<f64> = (0..12).map(|i| -3.0 + 6.0 * i as f64 / 11.0).collect();
    let mut dgen = StdRng::seed_from_u64(11);
    let y: Vec<f64> = x
        .iter()
        .map(|&xi| Normal::new(0.8 * xi - 0.4, 0.8).unwrap().sample(&mut dgen))
        .collect();

    // Reproducibility is a fugue value: a seed *is* the chain. Run four chains
    // and let split-R-hat referee whether they agree.
    let mut chains: Vec<Vec<Trace>> = Vec::new();
    for seed in [11u64, 12, 13, 14] {
        let mut rng = StdRng::seed_from_u64(seed);
        // 2000 kept draws, 500 warmup steps of proposal adaptation.
        let draws = adaptive_mcmc_chain(&mut rng, regression(x.clone(), y.clone()), 2_000, 500);
        chains.push(draws.into_iter().map(|(_, trace)| trace).collect());
    }

    // Split-R-hat (Vehtari et al. 2021) and ESS — the same convergence
    // diagnostics the widget shows.
    let rhat = r_hat_f64(&chains, &addr!("slope"));
    let slope_chain0: Vec<f64> = chains[0]
        .iter()
        .filter_map(|t| t.get_f64(&addr!("slope")))
        .collect();
    let ess = effective_sample_size_mcmc(&slope_chain0);

    println!("split-R-hat(slope) = {rhat:.3}   ESS(slope, chain 0) = {ess:.0}");
}

adaptive_mcmc_chain(&mut rng, model_fn, 2_000, 500) runs single-site Metropolis with a diminishing-adaptation schedule targeting 0.44 acceptance — the middle of the Goldilocks band from Things to try #4. It returns Vec<(A, Trace)>: the returned value and the full recording behind it. r_hat_f64 is the coral/green R̂ readout; effective_sample_size_mcmc is the ESS readout.

Go deeper

• Next explorable: Rolling, Not Guessing: HMC — the same regression, the same twelve points, but the walker rolls downhill with momentum instead of guessing, and covers the posterior in a fraction of the steps.
• Tutorial: Basic Inference — running MCMC end to end on a real model.
• How-to: Debugging Models — reading R̂ and ESS when a chain misbehaves.
• API: adaptive_mcmc_chain · r_hat_f64

Next: Rolling, Not Guessing: Hamiltonian Monte Carlo

Rolling, Not Guessing: Hamiltonian Monte Carlo

Random-walk Metropolis proposes blindly and hopes. Hamiltonian Monte Carlo does something smarter: it gives the sampler momentum and lets it roll across the posterior like a ball on a landscape, following the slope instead of guessing against it. One good roll crosses ground that a random walk needs hundreds of timid steps to cover.

Same problem as the Metropolis explorable: fit a straight line y = a·x + b to twelve noisy points. Left is data space — the yellow points and a fan of candidate fit lines. Right is parameter space — the posterior over (slope, intercept), with the sampler rolling through it. A point on the right is a line on the left: when the coral ball moves right, its coral line swings on the left.

These are the same twelve seeded points as the Metropolis page — run them next to each other. Every run is a replayable trace: fix the 11 seed and you get the exact same momenta and trajectories every time.

Things to try

1. Watch the violet leapfrog trajectory roll across parameter space — a single proposal travelling much farther than a random-walk hop — while its coral fit line swings across the data on the left. (It is already rolling; the controls let you pause, step, and steer.)
2. Drag the rightmost yellow point far up. The posterior heatmap tilts toward steeper slopes and the coral ball rolls after it within a few transitions. This linked deformation is the whole point of the page.
3. Turn on MH side-by-side. Both samplers get the same number of gradient-budget evaluations; the random walk's fits are the dim green spaghetti, HMC's the bright green — HMC fans across the plausible lines noticeably faster.
4. Push STEP ε up past 0.25. Trajectories start glowing coral and the divergence counter ticks — the integrator has gone unstable and every such proposal is rejected.
5. Set LEAPFROG L to 1. HMC collapses toward a random walk — the momentum never gets to carry the ball anywhere. Watch the energy strip lurch.

What you just saw

The target is an honest Bayesian linear regression — no banana, no toy density. Each observation is Gaussian around the line, with fixed noise σ_obs = 0.8, and both parameters get a Normal(0, 2.5) prior:

$$\log \pi (a,b)=\sum _i\log \mathcal{N}(y_i\mid a{x}_i+b,0.8)+\log \mathcal{N}(a\mid 0,2.5)+\log \mathcal{N}(b\mid 0,2.5)$$

The likelihood (yellow) pulls the line through the points; the prior (blue) keeps the coefficients sane; together they make the posterior (green) — the heatmap on the right.

HMC augments the position q = (a, b) with a fresh momentum p drawn from a Gaussian, then simulates a physical system whose total energy is the Hamiltonian:

$$H(q,p)=U(q)+K(p),U(q)=-\log \pi (q),K(p)=\frac{1}{2}{p}^{\top }{M}^{-1}p$$

The potential energy is the negative log-posterior — the green surface the ball rolls on. The kinetic energy is the momentum you flick it with. Low posterior density means high potential, so the ball is pulled toward the high-probability valley, exactly where you want samples.

The trajectory is integrated by the leapfrog scheme — a half-kick to momentum, a full drift in position, another half-kick:

$$p_{t+\frac{1}{2}}=p_t+\frac{\epsilon }{2}{\nabla }_q\log \pi (q_t)$$ $$q_{t+1}=q_t+\epsilon M^{-1}{p}_{t+\frac{1}{2}}$$ $$p_{t+1}=p_{t+\frac{1}{2}}+\frac{\epsilon }{2}{\nabla }_q\log \pi (q_{t+1})$$

Here the gradient ∇ log π is available in closed form (a Gaussian model), and the widget uses it exactly — that is why the leapfrog force points straight at the posterior mode. Leapfrog is reversible and volume-preserving, so the Metropolis correction at the end of the trajectory has no Jacobian term. It reduces to a comparison of total energy at the endpoints:

$$\alpha =\min (1,\exp (H(q,p)-H(q^{\prime },p^{\prime })\left)\right)$$

If the integrator were exact, H would be conserved and every proposal accepted. It is not exact, so a small energy error ΔH remains — that is the number on the strip chart. Keep it small (tune ε) and acceptance stays high. Let it explode and the proposal diverges: the trajectory shoots off, ΔH blows past any sane bound, and the sample is thrown away. Divergences are not a bug to hide — real samplers count and report them as a warning that the geometry is too sharp for the current step size.

The fugue code

Fugue ships HMC as hmc_chain. Models are ordinary Rust closures with no autodiff, so the force ∇ log π is computed by central finite differences — an approximate force, but an exact accept/reject against the true log-density, so the stationary distribution is exactly the posterior (the finite-difference step only costs efficiency, never correctness).

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

// The widget's model: fit y = a·x + b with fixed observation noise. The green
// heatmap is exactly this posterior over (slope, intercept); the coral ball is
// the current (a, b) rolling toward the values the yellow points support.
fn regression(xs: Vec<f64>, ys: Vec<f64>) -> Model<(f64, f64)> {
    prob! {
        let a <- sample(addr!("slope"),     Normal::new(0.0, 2.5).unwrap()); // blue prior
        let b <- sample(addr!("intercept"), Normal::new(0.0, 2.5).unwrap()); // blue prior
        // one yellow datum per point; sigma_obs fixed at 0.8
        let _obs <- plate!(i in xs.iter().zip(ys.iter()).enumerate() => {
            let (k, (x, y)) = i;
            observe(addr!("y", k), Normal::new(a * x + b, 0.8).unwrap(), *y)
        });
        pure((a, b))
    }
}

fn main() {
    // The same twelve seeded points the widget draws.
    let xs: Vec<f64> = (0..12).map(|i| -3.0 + 6.0 * i as f64 / 11.0).collect();
    let ys: Vec<f64> = xs.iter().map(|&x| 0.8 * x - 0.4).collect(); // + noise in practice

    let mut rng = StdRng::seed_from_u64(11); // a seed is a replayable trace

    // Mirror the widget's sliders: L leapfrog steps, an initial step size eps.
    let config = HMCConfig {
        n_leapfrog: 25,               // the LEAPFROG L slider
        init_step_size: Some(0.08),   // the STEP ε slider (warmup still tunes it
                                      // by dual averaging toward target_accept)
        target_accept: 0.8,           // Hoffman & Gelman's recommended target
        ..HMCConfig::default()        // finite_diff_eps = 1e-5, adapt_mass = false
    };

    // hmc_chain(rng, model_fn, n_samples, n_warmup, config) -> Vec<(A, Trace)>
    let model_fn = move || regression(xs.clone(), ys.clone());
    let samples = hmc_chain(&mut rng, model_fn, 1000, 500, config);

    let slopes: Vec<f64> = samples.iter()
        .filter_map(|(_, trace)| trace.get_f64(&addr!("slope")))
        .collect();
    let mean = slopes.iter().sum::<f64>() / slopes.len() as f64;
    println!("posterior mean slope ≈ {mean:.3}"); // near the true 0.8
}

The same ESS you watch converge in the widget is a real diagnostic. Run a few chains and combine them with the multi-chain estimators:

use fugue::*;

// After collecting `chains: Vec<Vec<Trace>>` from several hmc_chain runs:
let slope_series: Vec<Vec<f64>> = chains.iter()
    .map(|c| c.iter().filter_map(|t| t.get_f64(&addr!("slope"))).collect())
    .collect();

let ess = effective_sample_size(&slope_series[0]); // per-chain ESS
let rhat = r_hat_f64(&chains, &addr!("slope"));     // split-R̂: want < 1.01
println!("ESS = {ess:.0},  R̂ = {rhat:.3}");

hmc_chain holds any discrete sites fixed during the roll (a Metropolis-within-Gibbs treatment), so compose it with adaptive_mcmc_chain when a model mixes continuous and discrete latents.

Go deeper

• Compare by hand: Random Walks in Posterior Space — the same twelve points, sampled the slow way.
• Tutorial: Basic Inference covers when to reach for HMC over MH or SMC.
• API docs: hmc_chain and HMCConfig.
• Next: Particles That Tell Stories — inference that moves through time instead of space.

Next: Particles That Tell Stories

Particles That Tell Stories

A hidden thing moves through time. You never see it directly — only noisy glimpses, one per tick. How do you track it? You send a swarm of guesses forward, reward the ones that fit each new glimpse, and let the rest die. That swarm is a particle filter, and every guess is a little story about where the truth went.

The score is a state-space model; the particles are performers improvising the hidden path; each observation is a critic that reweights them. Watch time flow left to right.

The yellow dots are what you observe. The dashed ink line is the truth (you never get to see it during real inference). The blue circles are particles — each radius is that particle's weight. The particles are the machinery; the green is the answer. At the current time step, the green violin is the filtering distribution — the weighted swarm's full belief about \( x_t \) right now, drawn as a weighted-particle density. The green line traces the filtered mean over time, and the green band around it is that belief's ±1σ spread — how sure the filter is, tick by tick.

Things to try

1. Press Step a few times. Each tick does three things: propagate the particles forward, weight them against the new yellow dot, then resample if they have grown too uneven.
2. Drag particles down to 10 and press Play. The cloud collapses onto one or two lineages within a few steps — degeneracy, the failure mode SMC exists to fight.
3. Turn adaptive resample off and play. With no resampling, a single particle swallows almost all the weight; ESS / N crashes toward 1/N and the green estimate goes deaf.
4. Push obs noise up. Flatter likelihoods mean gentler reweighting, so ESS stays high and lineages survive longer — and watch the green band fatten: looser data means a less certain filter, and the filtering violin spreads to match.
5. Change the seed scrub, then reset. A seeded run is a replayable trace: same seed, same particles, same story, every time.

What you just saw

The hidden path is a random walk; each observation is a noisy read of it:

$$x_t\sim \mathcal{N}(x_{t-1},\sigma )y_t\sim \mathcal{N}(x_t,\tau )$$

The blue line is the transition — how a particle proposes its next position. The yellow line is the likelihood — how well that position explains the observation.

Propagate. Every particle draws its next state from the transition. On the canvas this is the drift step: the whole cloud slides one column right.

Weight. Each particle's weight is multiplied by how likely the new observation was under it:

$${\tilde{w}}_t^{(i)}\propto W_{t-1}^{(i)}p\left(y_t\mid x_t^{(i)}\right)$$

Particles near the yellow dot fatten; particles far from it shrink. Normalize and you have the new weights.

The filtering distribution. The weighted swarm is an estimate of one distribution: \( p(x_t \mid y_{1:t}) \), your belief about the hidden state given every observation so far. That is what the green violin draws — a weighted-particle kernel density of \( {x_t^{(i)}, W_t^{(i)}} \) at the current column. Its mean and ±1σ band are the one-number summaries you actually report:

$${\hat{\mu }}_t=\sum _i{W}_t^{(i)}{x}_t^{(i)}{\hat{\sigma }}_t^2=\sum _i{W}_t^{(i)}{\left(x_t^{(i)}-{\hat{\mu }}_t\right)}^2$$

Measure health. The effective sample size counts how many particles are really doing work:

$$\mathrm{ESS}=\frac{1}{{\sum }_i{\left(W_t^{(i)}\right)}^2}$$

All weight on one particle gives ESS = 1; perfectly even weights give ESS = N. The readout shows ESS / N, green above 0.5, coral below.

Resample. When ESS / N falls under the threshold, draw a fresh population by sampling parents in proportion to weight, then reset weights to uniform. Fat particles spawn duplicates (the violet fan); starved particles go extinct. This is the moment of the page: the swarm forgets its dead ends and concentrates where the evidence is. fugue's default threshold is exactly this — resample when ESS / N < 0.5.

Evidence, for free. Each weighting step also hands you a piece of the marginal likelihood. Summing the log of each step's mean weight gives an unbiased estimate of the log-evidence:

$$\log p(y_{1:T})=\sum _t\log \sum _i{W}_{t-1}^{(i)}p\left(y_t\mid x_t^{(i)}\right)$$

That number — the log-evidence readout — is what makes SMC a model-comparison tool, not merely a sampler.

The fugue code

fugue expresses the same state-space model as one Model. Each latent state depends on the previous one; each observe is a yellow dot.

use fugue::*;

const STEP: f64 = 0.7; // latent random-walk step σ
const OBS: f64 = 0.6; // observation noise τ

// x_0 ~ N(0,1); x_t ~ N(x_{t-1}, STEP); observe each y_t ~ N(x_t, OBS).
// Returns the whole latent path.
fn state_space(ys: Vec<f64>) -> Model<Vec<f64>> {
    let y0 = ys[0];
    let mut m: Model<Vec<f64>> = sample(addr!("x", 0), Normal::new(0.0, 1.0).unwrap())
        .bind(move |x0| {
            observe(addr!("y", 0), Normal::new(x0, OBS).unwrap(), y0).map(move |_| vec![x0])
        });
    for t in 1..ys.len() {
        let yt = ys[t];
        m = m.bind(move |xs| {
            let prev = *xs.last().unwrap(); // depend on the previous state
            sample(addr!("x", t), Normal::new(prev, STEP).unwrap()).bind(move |xt| {
                observe(addr!("y", t), Normal::new(xt, OBS).unwrap(), yt) // the yellow dot
                    .map(move |_| {
                        let mut xs = xs;
                        xs.push(xt);
                        xs
                    })
            })
        });
    }
    m
}

Run Sequential Monte Carlo over it. The config is the same knobs you just played with:

use fugue::*;
use rand::rngs::StdRng;
use rand::SeedableRng;

fn main() {
    let ys: Vec<f64> = vec![/* your observation series */];

    let mut rng = StdRng::seed_from_u64(42);
    let config = SMCConfig {
        resampling_method: ResamplingMethod::Systematic, // fugue's default
        ess_threshold: 0.5,      // resample when ESS / N < 0.5 — the green/coral band
        rejuvenation_steps: 1,   // MCMC moves that restore diversity after a resample
    };

    let result = adaptive_smc(&mut rng, 200, || state_space(ys.clone()), config);

    println!("particles:    {}", result.particles.len());
    println!("ESS:          {:.1}", effective_sample_size(&result)); // the ESS readout
    println!("log-evidence: {:.3}", result.log_evidence); // unbiased log p(y_1:T)
}

Every name here is real: adaptive_smc, SMCConfig, ResamplingMethod::Systematic, effective_sample_size, and SMCResult::log_evidence all live in src/inference/smc.rs. SMCResult dereferences to Vec<Particle>, so effective_sample_size(&result) works directly on the returned population.

Go deeper

• Tutorial: Sequential Monte Carlo — the full API, resampling methods, and rejuvenation.
• API: adaptive_smc and effective_sample_size.
• Next explorable: A Field Guide to Distributions — the building blocks every model above is made of.

Next: A Field Guide to Distributions

A Field Guide to Distributions

A distribution is a shape and a rule for drawing from it. The shape is its density — where values are likely to land. The rule is its sampler — how a stream of random numbers becomes draws. Every distribution below is one note in fugue's vocabulary; pick one and hear it.

The blue curve is the exact law. The green bars are samples piling up. Watch them converge — that is the law of large numbers, live.

Things to try

1. Start on Normal and drag σ up to 3 — the bell flattens and spreads, but the green samples still trace it.
2. Switch to Beta and set α and β both below 1 — the density turns into a U, piling mass at the two edges.
3. Open Cauchy. Its mean and variance readouts both read —: the tails are so heavy that neither integral converges. The samples wander wildly.
4. Pick Bernoulli and read the sample → badge: it says bool. Fugue gives you a real boolean, not a float you have to compare against 1.0.
5. Drag Uniform's high below its low. The canvas turns red — this is exactly the Err that Uniform::new returns for an invalid interval.
6. On any distribution, drop the seed back to a value you already used. The green histogram redraws the identical run: a seeded stream is a replayable trace.

What you just saw

For a continuous distribution the blue curve is the probability density f(x). It is not a probability — it is a density, so it can exceed 1. What integrates to 1 is the area:

$${\int }_{\mathcal{S}}f(x)dx=1$$

For a discrete distribution the blue stems are the probability mass P(x), one bar per outcome, and the bars sum to 1:

$$\sum _{x\in \mathcal{S}}P(x)=1$$

The green bars are an empirical estimate of that same law built from samples. As the sample count grows they converge on the blue — the law of large numbers. The coral line is a single query: it reports $\log f(x)$ at the x you drag to. Inference works in log-space because a product of thousands of these densities underflows to zero in ordinary floating point; a sum of their logs does not.

The markers name three summaries of the shape: the mean (the balance point), the median (half the mass on each side), and the mode (the peak). For a skewed law like LogNormal they sit in different places; for a symmetric one they coincide.

Fugue's type story

Most PPLs make every draw a f64, so a coin flip comes back as 1.0 and you compare floats; a count comes back as 4.0 and you round it. Fugue draws return their natural type, checked at compile time.

use fugue::*;
use rand::thread_rng;

fn main() {
    let mut rng = thread_rng();

    // Bernoulli -> bool. The coral query line is `log_prob` on one outcome.
    let coin = Bernoulli::new(0.5).unwrap();
    let heads: bool = coin.sample(&mut rng);
    let lp_true: f64 = coin.log_prob(&true); // ln(0.5)
    if heads { /* no `== 1.0`; it's already a bool */ }

    // Poisson -> u64. Counts are counts, not rounded floats.
    let arrivals: u64 = Poisson::new(4.0).unwrap().sample(&mut rng);
    let total_wait = arrivals * 10; // integer arithmetic, no cast

    // Categorical -> usize. Safe to index an array with, by construction.
    let pick: usize = Categorical::new(vec![0.5, 0.3, 0.2]).unwrap().sample(&mut rng);
    let labels = ["a", "b", "c"];
    let chosen = labels[pick];

    // Continuous laws return f64, as expected.
    let x: f64 = Normal::new(0.0, 1.0).unwrap().sample(&mut rng);
    let density: f64 = Normal::new(0.0, 1.0).unwrap().log_prob(&x);

    println!("{heads} {arrivals} {total_wait} {chosen} {x:.3} {density:.3} {lp_true:.3}");
}

Every constructor is fallible: Normal::new, Beta::new, Uniform::new and the rest return Result, so an invalid parameter (a negative σ, a low ≥ high) is an Err you handle, never a silent NaN. That is the red screen in try #5.

Inside a model, a draw is a sample at a named address; its distribution rides along and fugue scores it for you:

#![allow(unused)]
fn main() {
use fugue::*;

// Estimate a coin's bias, then predict the next flip's count over 10 tosses.
let model = prob!(
    let p <- sample(addr!("bias"), Beta::new(2.0, 2.0).unwrap());
    let hits <- sample(addr!("hits"), Binomial::new(10, p).unwrap());
    pure(hits)
);
}

The full field guide

Every distribution in fugue. The sample → column is the natural return type.

Continuous

Discrete

Go deeper

• Tutorial: Working with Distributions — the same catalog in prose, with modeling patterns.
• API: fugue::core::distribution — every constructor, its constraints, and its log_prob.
• Next explorable: Anatomy of a Probabilistic Program — put a Beta prior and Bernoulli data together and watch Bayes multiply them.

Next: Explorables

Playground

Everything on this page runs the real fugue crate, compiled to WebAssembly. The editor below speaks a subset of the prob! macro language; when you press Run, your model text is parsed in Rust, folded into actual Model combinators, and handed to the same inference kernels this book documents — adaptive Metropolis–Hastings one transition at a time, Hamiltonian Monte Carlo with dual-averaging warmup, adaptive tempered SMC with a log-evidence estimate. The draws you watch stream in are fugue's draws, not a JavaScript imitation of them.

What the editor understands

The playground accepts the statement forms of prob!, with data arrays provided as JSON (an object of named arrays, or a bare array bound to data):

let p <- sample(addr!("p"), Beta(2.0, 2.0));   // sample a latent
let m = 2.0 * p - 1.0;                          // deterministic let
observe(addr!("y"), Normal(m, 0.8), 1.4);       // condition on data
factor(-0.5 * m * m);                           // add a log-weight
for i in 0..data.len() { ... }                  // plates
pure(p)                                         // the model's return value

Rust spellings from the docs paste in unchanged — Normal::new(0.0, 1.0).unwrap() parses the same as Normal(0.0, 1.0). Available distributions: Normal, Uniform, LogNormal, Exponential, Beta, Gamma, InverseGamma, StudentT, Cauchy, Laplace, Weibull, ChiSquared, Bernoulli, Binomial, Poisson, Categorical, DiscreteUniform — the same constructors, parameterizations, and support checks as the crate, because they are the crate. Expressions know + - * /, exp, ln, sqrt, abs, pow, min, max, floor, sin, cos, tanh, array indexing y[i], and .len().

Two honest limitations: parameter regions that would make a distribution invalid (say, a negative scale reached through your arithmetic) score as impossible (-inf) rather than erroring, exactly how the samplers treat leaving a support; and the interpreter covers the statement subset above, not arbitrary Rust — for that, there is cargo add fugue-ppl.

Things to try

1. Watch R̂ earn its keep. Run the coin-flip preset with 3 chains: the chains start from independent prior draws, disagree, then R̂ falls toward 1.00 as their histograms merge — the same split-R̂ computation fugue::inference::diagnostics runs in production.
2. Break a model, kindly. Change a flip observation to 2.0 (a coin that landed on its edge?). Bernoulli observes a boolean, so anything nonzero coerces to true — then make the prior Beta(0.5, 0.5) and watch the posterior bend toward the edges.
3. Compare kernels on the same posterior. Run the eight-schools preset under MH, then under HMC. Single-site MH updates one coordinate at a time and mixes slowly through the funnel; HMC moves all ten coordinates jointly. The ESS readout is the receipt.
4. Ask for evidence. Switch to SMC on any preset: you get weighted posterior particles and log Z, the marginal likelihood MH and HMC never give you.
5. Make it yours. Replace the data JSON with your own numbers. The model recompiles as you type — the ✓/✗ line is fugue's parser talking.

How this works

crates/fugue-wasm (in the fugue repository) exposes wasm-bindgen entry points over the crate: compile a model source + data payload, then drive step(n) per animation frame and read draws and diagnostics back as typed arrays. Sampler state lives in Rust — traces, adaptation state, tuned step sizes — and every run is seeded, so a seed is a replayable recording here too. The interactive figures in the explorables chapters use the same package for their samplers; only their canvas rendering is JavaScript.

Go deeper

• Explorables: Random Walks in Posterior Space, Rolling, Not Guessing: HMC, Particles That Tell Stories — the guided tours of the kernels this playground exposes.
• The real thing: Installation — cargo add fugue-ppl and the full prob! language, typed traces, custom handlers.
• API: docs.rs/fugue-ppl.

How-To Guides

The How-To guides provide practical, task-oriented instructions for accomplishing specific goals with Fugue. Unlike tutorials that teach concepts step-by-step, these guides assume you understand the basics and want to solve particular problems efficiently.

Guide Overview

These guides are designed to be example-first and immediately actionable. Each guide includes comprehensive, executable code examples that serve as the canonical source of truth for the patterns they demonstrate.

📊 Working with Distributions

When to use: You need to understand Fugue's type-safe distribution system, parameter validation, or probability calculations.

What you'll learn:

• Type-safe distribution usage (bool, u64, usize, f64 return types)
• Parameter validation and error handling
• Continuous vs discrete distribution patterns
• Categorical distributions and safe indexing
• Distribution composition and practical modeling
• Probability calculations and testing strategies

Key patterns: Natural return types, parameter validation, distribution testing

🏗️ Building Complex Models

When to use: You want to compose sophisticated probabilistic models using Fugue's macro system and advanced patterns.

What you'll learn:

• prob! macro for do-notation style probabilistic programming
• plate! macro for vectorized operations and array processing
• scoped_addr! macro for hierarchical address management
• Model composition and sequential dependencies
• Hierarchical modeling patterns
• Bayesian linear regression and mixture models

Key patterns: Monadic composition, vectorization, hierarchical structure

⚡ Optimizing Performance

When to use: Your models need to run efficiently in production or handle large-scale inference workloads.

What you'll learn:

• Numerical stability with log-space computations
• Batch processing and vectorized model patterns
• Performance monitoring and measurement

Key patterns: Numerical stability, batch processing

🔍 Debugging Models

When to use: Your probabilistic models aren't behaving as expected, or you need to diagnose issues in inference.

What you'll learn:

• Comprehensive trace inspection and analysis
• Type-safe value access with proper error handling
• Model validation against analytical solutions
• Safe vs strict handler usage for error resilience
• MCMC diagnostics and convergence assessment
• Performance and memory debugging techniques

Key patterns: Trace analysis, validation testing, diagnostic metrics

🎛️ Custom Handlers

When to use: You need to extend Fugue's execution model with custom behavior, logging, or specialized inference algorithms.

What you'll learn:

• Complete Handler trait implementation
• Decorator pattern for cross-cutting concerns
• Stateful handlers for analytics and monitoring
• Conditional filtering and value modification
• Performance monitoring integration
• Custom inference algorithms (MCMC-like patterns)
• Handler composition and chaining

Key patterns: Algebraic effects, decorator composition, custom inference

🚀 Production Deployment

When to use: You're deploying probabilistic models to production environments and need reliability, monitoring, and operational excellence.

What you'll learn:

• Error handling and graceful degradation patterns
• Circuit breaker implementation for fault tolerance
• Configuration management for multiple environments
• Comprehensive metrics and Prometheus integration
• Automated health checks and system validation
• Input validation and security best practices
• Deployment strategies (blue-green, canary, rolling)

Key patterns: Fault tolerance, observability, operational readiness

How to Use These Guides

🎯 Task-Oriented Approach

Each guide focuses on solving specific problems:

• Need to understand distributions? → Start with "Working with Distributions"
• Building complex models? → "Building Complex Models" has the macros and patterns
• Performance issues? → "Optimizing Performance" covers memory and numerical techniques
• Models not working? → "Debugging Models" provides diagnostic approaches
• Need custom behavior? → "Custom Handlers" shows how to extend the system
• Going to production? → "Production Deployment" covers operational concerns

📚 Progressive Complexity

The guides are ordered by increasing complexity:

graph TD
    A["📊 Working with<br/>Distributions"] --> B["🏗️ Building Complex<br/>Models"]
    B --> C["⚡ Optimizing<br/>Performance"]
    A --> D["🚀 Production<br/>Deployment"]
    E["🔍 Debugging<br/>Models"] --> F["🎛️ Custom<br/>Handlers"]
    F --> D
    C --> D

• Start with distributions and model building
• Add performance optimization when needed
• Use debugging when things go wrong
• Extend with custom handlers for specialized needs
• Deploy with production patterns for real applications

🔄 Cross-References

Guides frequently reference each other:

• Performance optimization builds on complex models
• Debugging techniques apply to all model types
• Custom handlers can incorporate performance patterns
• Production deployment uses patterns from all previous guides

📝 Example-First Philosophy

Every guide follows the same structure:

1. Executable examples as the source of truth
2. Comprehensive code snippets with anchor tags
3. Practical explanations of when and how to use patterns
4. Testing strategies to verify correctness
5. Best practices learned from real-world usage

Code Examples

All code examples in these guides are:

• ✅ Executable: Run with cargo run --example <guide_name>
• ✅ Tested: Verified with cargo test --examples
• ✅ Documented: Included via {{#include}} from example files
• ✅ Comprehensive: Cover real-world usage patterns
• ✅ Type-Safe: Leverage Rust's type system throughout

Quick Reference

Integration with Other Documentation

🚀 Getting Started → How-To Guides

After completing the Getting Started tutorials, use these guides to solve specific problems in your own projects.

📖 How-To Guides → Tutorials

For deeper conceptual understanding, see the Tutorials section which provides comprehensive examples of complete applications.

🔧 How-To Guides → API Reference

For detailed API documentation, consult the API Reference section for specific functions and types mentioned in these guides.

Contributing

When adding new How-To guides:

1. Create executable examples first in examples/
2. Use anchor tags to mark code sections for inclusion
3. Write the guide using {{#include}} for all code snippets
4. Test thoroughly with both cargo test and mdbook test
5. Update this README with the new guide information

Each guide should solve specific, practical problems that users commonly encounter when working with Fugue in real applications.

Working with Distributions

Fugue's type-safe distribution system represents a principled approach to probabilistic programming, eliminating entire classes of runtime errors through rigorous type theory while preserving the full expressiveness of statistical modeling. This guide demonstrates the mathematical foundations and practical applications of Fugue's distribution architecture.

Fugue ships 17 distributions in fugue::core::distribution: 12 continuous (Normal, Uniform, LogNormal, Exponential, Beta, Gamma, StudentT, Cauchy, Laplace, Weibull, ChiSquared, InverseGamma) and 5 discrete/other (Bernoulli, Categorical, Binomial, Poisson, DiscreteUniform).

Type Safety in Practice

Traditional probabilistic programming libraries return f64 for everything, leading to casting overhead and runtime errors. Fugue distributions return their natural types:

    // Demonstrate natural return types
    let coin = Bernoulli::new(0.5).unwrap();
    let flip: bool = coin.sample(&mut rng); // Natural boolean

    if flip {
        println!("🪙 Coin flip result: Heads!");
    } else {
        println!("🪙 Coin flip result: Tails!");
    }

No casting, no comparisons with floating-point values—just natural boolean logic.

Continuous Distributions

Continuous distributions in Fugue model phenomena over uncountable domains $\mathcal{S}\subseteq {\mathbb{R}}^n$. The probability density function $f_X(x)$ satisfies the normalization condition:

$${\int }_{\mathcal{S}}{f}_X(x)dx=1$$

For computational stability, Fugue operates in log-space by default, computing $\log f_X(x)$ to avoid numerical underflow:

    // Working with continuous distributions
    let standard_normal = Normal::new(0.0, 1.0).unwrap();
    let sample: f64 = standard_normal.sample(&mut rng);
    println!("📊 Standard normal sample: {:.3}", sample);

    // Compute log-probability density
    let log_density = standard_normal.log_prob(&0.0); // Peak of standard normal
    println!(
        "📈 Log-density at x=0: {:.3} (peak of standard normal)",
        log_density
    );

    // Custom parameters
    let measurement_model = Normal::new(10.0, 0.5).unwrap();
    let measurement = measurement_model.sample(&mut rng);
    println!("🔬 Sensor measurement (μ=10.0, σ=0.5): {:.3}", measurement);

Key Points:

• sample() returns f64 for direct arithmetic
• log_prob() computes log-density (avoids numerical underflow)
• Parameter validation happens at construction time

Discrete Distributions

Discrete distributions operate over countable support sets $\mathcal{S}\subseteq {\mathbb{Z}}^n$ or finite sets. The probability mass function satisfies:

$$\sum _{x\in \mathcal{S}}P(X=x)=1$$

Fugue enforces this constraint at construction time and leverages natural integer types to eliminate precision loss from floating-point representation:

    // Working with discrete distributions

    // Count data
    let event_rate = Poisson::new(3.0).unwrap();
    let count: u64 = event_rate.sample(&mut rng);
    println!("📅 Event count (λ=3.0): {} events", count);

    // Log-probability mass
    let prob_3_events = event_rate.log_prob(&3);
    println!(
        "🎯 Log-probability of exactly 3 events: {:.3}",
        prob_3_events
    );

    // Use counts directly in calculations
    let total_cost = count * 50; // Direct arithmetic with u64
    println!("💰 Total cost ({} events × $50): ${}", count, total_cost);

Benefits:

• u64 counts support direct arithmetic without casting
• No precision loss from floating-point representations
• Natural integration with Rust's type system

Safe Categorical Sampling

Categorical distributions return usize for safe array indexing:

    // Safe categorical sampling
    let choices = vec![0.3, 0.5, 0.2]; // Three categories
    let categorical = Categorical::new(choices).unwrap();
    let selected: usize = categorical.sample(&mut rng);

    // Safe array indexing (no bounds checking needed)
    let options = ["Option A", "Option B", "Option C"];
    println!(
        "🎲 Categorical choice (weights: 0.3, 0.5, 0.2): {}",
        options[selected]
    );

    // Uniform categorical
    let uniform_choice = Categorical::uniform(5).unwrap();
    let idx: usize = uniform_choice.sample(&mut rng);
    println!("🎯 Uniform random index (0-4): {}", idx);

Parameter Validation

Fugue enforces mathematical constraints through compile-time and runtime validation. Each distribution family ${\mathcal{D}}_{\theta }$ has a parameter space $\Theta$ defining valid configurations:

graph TD
    A[Parameter Input θ] --> B{θ ∈ Θ?}
    B -->|Yes| C[Distribution Construction]
    B -->|No| D[ValidationError]
    C --> E[Type-Safe Sampling]
    D --> F[Early Failure Detection]

Constraint Examples:

• Normal Distribution: $\mathcal{N}(\mu ,{\sigma }^2)$ requires $\sigma >0$
• Beta Distribution: $\text{Beta}(\alpha ,\beta )$ requires $\alpha ,\beta >0$
• Categorical Distribution: ${\sum }_i{p}_i=1$ and $p_i\ge 0\forall i$

    // Distribution parameter validation

    // This will return an error
    match Normal::new(0.0, -1.0) {
        Ok(_) => println!("✅ Normal(μ=0.0, σ=-1.0) created successfully"),
        Err(e) => println!("❌ Normal(μ=0.0, σ=-1.0) failed: {:?}", e),
    }

    // Beta distribution parameters must be positive
    match Beta::new(0.0, 1.0) {
        Ok(_) => println!("✅ Beta(α=0.0, β=1.0) created successfully"),
        Err(e) => println!("❌ Beta(α=0.0, β=1.0) failed: {:?}", e),
    }

    // Poisson rate must be non-negative
    match Poisson::new(-1.0) {
        Ok(_) => println!("✅ Poisson(λ=-1.0) created successfully"),
        Err(e) => println!("❌ Poisson(λ=-1.0) failed: {:?}", e),
    }

Validation Rules:

• Normal: σ > 0
• Beta: α > 0, β > 0
• Poisson: λ ≥ 0
• Categorical: probabilities sum to 1, all non-negative

Storing Mixed Distributions

Use trait objects for collections of distributions with the same return type:

    // Storing different distributions together
    let continuous_dists: Vec<Box<dyn Distribution<f64>>> = vec![
        Normal::new(0.0, 1.0).unwrap().clone_box(),
        Beta::new(2.0, 5.0).unwrap().clone_box(),
        Uniform::new(-1.0, 1.0).unwrap().clone_box(),
    ];

    // Sample from each
    for (i, dist) in continuous_dists.iter().enumerate() {
        let sample = dist.sample(&mut rng);
        let dist_name = match i {
            0 => "Normal(0,1)",
            1 => "Beta(2,5)",
            2 => "Uniform(-1,1)",
            _ => "Unknown",
        };
        println!("📦 {} sample: {:.3}", dist_name, sample);
    }

This enables dynamic distribution selection and model composition patterns.

Practical Modeling Patterns

Common modeling scenarios demonstrate natural type usage:

    // Practical modeling examples

    // Model a sensor with noise
    let true_temperature = 20.5; // True value
    let sensor_noise = Normal::new(0.0, 0.2).unwrap(); // Measurement error
    let measured_temp = true_temperature + sensor_noise.sample(&mut rng);
    println!(
        "🌡️  True temperature: {:.2}°C → Measured: {:.2}°C",
        true_temperature, measured_temp
    );

    // Count model for arrivals
    let arrival_rate = Poisson::new(2.5).unwrap(); // 2.5 arrivals per hour
    let hourly_arrivals = arrival_rate.sample(&mut rng);
    println!(
        "🚪 Expected arrivals: 2.5/hour → Actual: {} arrivals",
        hourly_arrivals
    );

    // Decision model
    let decision_prob = 0.7;
    let decision = Bernoulli::new(decision_prob).unwrap();
    let will_buy = decision.sample(&mut rng);
    if will_buy {
        println!("🛒 Customer decision (p=0.7): Will make a purchase");
    } else {
        println!("🚶 Customer decision (p=0.7): Will not purchase");
    }

Each distribution serves its natural domain without artificial conversions.

Working with Log-Probabilities

Logarithmic probability computation is essential for numerical stability in probabilistic programming. Consider the log-sum-exp operation for computing:

$$\log \left(\sum _{i=1}^n{e}^{x_i}\right)=x_{\max }+\log \left(\sum _{i=1}^n{e}^{x_i-x_{\max }}\right)$$

where $x_{\max }={\max }_i{x}_i$. This formulation prevents overflow when $|x_i|$ is large:

    // Working with log-probabilities

    let normal = Normal::new(100.0, 15.0).unwrap();

    // Multiple observations
    let observations = vec![98.5, 102.1, 99.8, 101.5, 97.2];
    let mut total_log_prob = 0.0;

    println!(
        "📋 Evaluating {} observations under Normal(μ=100.0, σ=15.0):",
        observations.len()
    );
    for obs in &observations {
        let log_p = normal.log_prob(obs);
        total_log_prob += log_p;
        println!("   x={}: log P(x) = {:.3}", obs, log_p);
    }

    println!("🔢 Joint log-probability: {:.3}", total_log_prob);

    // Convert back to probability (be careful with underflow!)
    if total_log_prob > -700.0 {
        // Avoid underflow
        let probability = total_log_prob.exp();
        println!("📊 Joint probability: {:.2e}", probability);
    } else {
        println!("⚠️  Joint probability too small to represent as f64");
    }

Advanced Patterns

For complex modeling scenarios, see these patterns:

Hierarchical Models

    // Hierarchical prior structure
    let global_mean = Normal::new(0.0, 10.0).unwrap();
    let mu = global_mean.sample(&mut rng);

    let group_precision = Gamma::new(2.0, 0.5).unwrap();
    let tau = group_precision.sample(&mut rng);
    let sigma = (1.0 / tau).sqrt(); // Convert precision to std dev

    // Individual observations from hierarchical model
    let individual = Normal::new(mu, sigma).unwrap();
    let observation = individual.sample(&mut rng);

    println!("🌐 Global mean: {:.3}", mu);
    println!("📊 Group std dev: {:.3}", sigma);
    println!("👤 Individual observation: {:.3}", observation);

Mixture Components

    // Mixture model components
    let mixture_weights = vec![0.6, 0.3, 0.1];
    let component_selector = Categorical::new(mixture_weights).unwrap();
    let selected_component: usize = component_selector.sample(&mut rng);

    // Different components
    let components = [
        Normal::new(-2.0, 0.5).unwrap(),
        Normal::new(0.0, 1.0).unwrap(),
        Normal::new(3.0, 0.8).unwrap(),
    ];

    let sample = components[selected_component].sample(&mut rng);
    println!(
        "🎯 Selected component {}: sample = {:.3}",
        selected_component, sample
    );

Conjugate Priors

    // Beta-Bernoulli conjugacy
    let prior_alpha = 2.0;
    let prior_beta = 8.0;
    let prior = Beta::new(prior_alpha, prior_beta).unwrap();
    let p: f64 = prior.sample(&mut rng);

    // Simulate some trials
    let trials = 20;
    let mut successes = 0;
    let bernoulli = Bernoulli::new(p).unwrap();

    for _ in 0..trials {
        if bernoulli.sample(&mut rng) {
            successes += 1;
        }
    }

    // Posterior parameters (conjugate update)
    let posterior_alpha = prior_alpha + successes as f64;
    let posterior_beta = prior_beta + (trials - successes) as f64;
    let posterior = Beta::new(posterior_alpha, posterior_beta).unwrap();
    let updated_p = posterior.sample(&mut rng);

    println!("🎲 Prior p: {:.3}", p);
    println!("📈 Observed: {}/{} successes", successes, trials);
    println!("🔄 Posterior p: {:.3}", updated_p);

Testing Your Distributions

Always test distribution properties and parameter validation:

    #[test]
    fn test_distribution_properties() {
        let mut rng = thread_rng();

        // Test type safety
        let coin = Bernoulli::new(0.5).unwrap();
        let flip: bool = coin.sample(&mut rng);
        assert!(flip == true || flip == false); // Must be boolean

        // Test parameter validation
        assert!(Normal::new(0.0, -1.0).is_err());
        assert!(Beta::new(0.0, 1.0).is_err());
        assert!(Poisson::new(-1.0).is_err());

        // Test valid distributions
        let normal = Normal::new(0.0, 1.0).unwrap();
        let sample = normal.sample(&mut rng);
        let log_prob = normal.log_prob(&sample);
        assert!(log_prob.is_finite());
    }

Common Pitfalls

1. Underflow in probability space: Use log-probabilities for accumulation
2. Parameter validation: Check constructor errors, don't assume success
3. Precision with counts: Use u64 return types directly, avoid f64 conversion
4. Categorical indexing: Trust the usize return—it's guaranteed valid

Next Steps

• Complex Models: See Building Complex Models for compositional patterns
• Debugging: Check out Debugging Models for troubleshooting
• Custom Logic: Learn Custom Handlers for specialized inference

The type-safe distribution system eliminates entire classes of runtime errors while making statistical code more readable and maintainable.

Building Complex Models

Fugue's compositional architecture is grounded in category theory and monadic structures, enabling the systematic construction of sophisticated probabilistic models through principled composition operators. This guide explores the mathematical foundations and practical applications of Fugue's macro system for building complex probabilistic programs.

Do-Notation with prob!

The prob! macro implements monadic do-notation for probabilistic computations, providing a natural syntax for sequential dependence. Formally, it translates:

$$\begin{array}{rl}\text{prob!}\{& \\ & x\leftarrow {\mathcal{M}}_1\\ & y\leftarrow {\mathcal{M}}_2(x)\\ & \text{pure}(f(x,y))\end{array}\}$$

into the monadic composition ${\mathcal{M}}_1\gg =\lambda x.{\mathcal{M}}_2(x)\gg =\lambda y.\eta (f(x,y))$:

graph LR
    A[M₁] -->|bind| B[λx.M₂⁽ˣ⁾]
    B -->|bind| C[λy.η⁽f⁽ˣ'ʸ⁾⁾]
    C --> D[Result⁽ᶻ⁾]

    // Simple do-notation style probabilistic program
    let _simple_model = prob!(
        let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
        let y <- sample(addr!("y"), Normal::new(x, 0.5).unwrap());
        let sum = x + y;  // Regular variable assignment
        pure(sum)
    );
    println!("✅ Created simple model with prob! macro");

Key Features:

• <- for probabilistic binding (monadic bind)
• = for regular variable assignment
• pure() to lift deterministic values
• Natural control flow without callback nesting

Vectorized Operations with plate!

The plate! macro implements plate notation from graphical models, representing conditionally independent replications. Given $N$ independent observations, plate notation expresses:

$$P(\mathbf{x}\mid \theta )=\prod _{i=1}^{N}P(x_i\mid \theta )$$

The computational graph shows the independence structure:

graph TB
    subgraph "Plate: i ∈ {1..N}"
        A[θ] --> B1[x₁]
        A --> B2[x₂]
        A --> B3[...]
        A --> BN[xₙ]
    end

    // Independent samples using plate notation
    let _vector_model = plate!(i in 0..5 => {
        sample(addr!("sample", i), Normal::new(0.0, 1.0).unwrap())
    });
    println!("✅ Created vectorized model with {} samples", 5);

    // Plate with observations
    let observations = [1.2, -0.5, 2.1, 0.8, -1.0];
    let n_obs = observations.len();
    let _observed_model = plate!(i in 0..n_obs => {
        observe(addr!("obs", i), Normal::new(0.0, 1.0).unwrap(), observations[i])
    });
    println!("✅ Created observation model for {} data points", n_obs);

Benefits:

• Automatic address indexing prevents conflicts
• Natural iteration over data structures
• Vectorized likelihood computations
• Clear intent for independent operations

Hierarchical Address Management

Complex models require systematic parameter organization following a tree-structured address space. The address hierarchy $\mathcal{A}$ forms a prefix tree where each node represents a scope:

$$\mathcal{A}=\{a_1,a_1.a_2,a_1.a_2.a_3,\dots \}$$

This hierarchical structure prevents address collisions and enables efficient parameter lookup:

    // Hierarchical model using scoped addresses
    let _hierarchical_model = prob!(
        let global_mu <- sample(addr!("global_mu"), Normal::new(0.0, 10.0).unwrap());
        let group_mu <- sample(scoped_addr!("group", "mu", "{}", 0),
                              Normal::new(global_mu, 1.0).unwrap());
        pure((global_mu, group_mu))
    );
    println!("✅ Created hierarchical model with scoped addresses");

Address Strategy:

• scoped_addr! prevents parameter name collisions
• Hierarchical structure mirrors model dependencies
• Systematic naming aids debugging and introspection
• Indices enable parameter arrays

Sequential Dependencies

Sequential models exhibit temporal dependence where the state at time $t$ depends on previous states. This creates a Markov chain structure:

$$P(x_{1:T})=P(x_1)\prod _{t=2}^{T}P(x_t\mid x_{t-1})$$

The computational challenge lies in maintaining state consistency while enabling efficient inference:

    // Sequential model with dependencies
    let _sequential_model = prob! {
        let states <- plate!(t in 0..3 => {
            sample(addr!("x", t), Normal::new(0.0, 1.0).unwrap())
                .bind(move |x_t| {
                    observe(addr!("y", t), Normal::new(x_t, 0.5).unwrap(), 1.0 + t as f64)
                        .map(move |_| x_t)
                })
        });

        pure(states)
    };
    println!("✅ Created sequential model with observations");

Patterns:

• Explicit state threading through computations
• Observation conditioning at each time step
• Autoregressive dependencies
• Mixed probabilistic and deterministic updates

Composable Model Functions

Build reusable model components:

    // Helper function to create a component model
    fn create_normal_component(name: &str, mean: f64, std: f64) -> Model<f64> {
        sample(addr!(name), Normal::new(mean, std).unwrap())
    }

    // Compose multiple components
    let _composition_model = prob! {
        let param1 <- create_normal_component("param1", 0.0, 1.0);
        let param2 <- create_normal_component("param2", 2.0, 0.5);
        let combined = param1 * param2;
        pure(combined)
    };
    println!("✅ Created composed model with reusable components");

Design Principles:

• Functions return Model<T> for composability
• Pattern matching enables model selection
• Pure functions for deterministic transformations
• Higher-order functions for model templates

Advanced Address Patterns

For large-scale models like neural networks:

    // Complex addressing for large models
    let _neural_layer_model = plate!(layer in 0..3 => {
        let layer_size = match layer {
            0 => 4,
            1 => 8,
            2 => 1,
            _ => 1,
        };

        plate!(i in 0..layer_size => {
            sample(
                scoped_addr!("layer", "weight", "{}_{}", layer, i),
                Normal::new(0.0, 0.1).unwrap()
            )
        })
    });
    println!("✅ Created neural network parameter structure");

Scaling Strategies:

• Systematic parameter naming conventions
• Multi-level scoping for complex architectures
• Consistent indexing schemes
• Hierarchical parameter organization

Mixing Styles for Flexibility

Combine macros with traditional function composition:

    // Mixture model with component selection
    let _mixture_model = prob! {
        let component <- sample(addr!("component"), Bernoulli::new(0.3).unwrap());
        let mu = if component { -2.0 } else { 2.0 };
        let x <- sample(addr!("x"), Normal::new(mu, 1.0).unwrap());
        pure((component, x))
    };
    println!("✅ Created mixture model with 2 components");

Best Practices:

• Use functions for reusable components
• Use macros for readable composition
• Separate concerns (priors, likelihood, observations)
• Document parameter dependencies

Real-World Applications

Bayesian Linear Regression

Bayesian linear regression models the relationship $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$ with uncertainty quantification:

$$\begin{array}{rl}\mathbf{\beta }& \sim \mathcal{N}({\mathbf{\mu }}_0,{\mathbf{\Sigma }}_0)\\ {\sigma }^2& \sim \text{InverseGamma}(\alpha ,\beta )\\ y_i\mid {\mathbf{x}}_i,\mathbf{\beta },{\sigma }^2& \sim \mathcal{N}({\mathbf{x}}_i^T\mathbf{\beta },{\sigma }^2)\end{array}$$

    // Complete Bayesian linear regression
    let x_data = [1.0, 2.0, 3.0, 4.0, 5.0];
    let y_data = [2.1, 3.9, 6.2, 8.1, 9.8];
    let n = x_data.len();

    let _regression_model = prob! {
        let intercept <- sample(addr!("intercept"), Normal::new(0.0, 10.0).unwrap());
        let slope <- sample(addr!("slope"), Normal::new(0.0, 10.0).unwrap());
        let precision <- sample(addr!("precision"), Gamma::new(1.0, 1.0).unwrap());
        let sigma = (1.0 / precision).sqrt();

        let _likelihood <- plate!(i in 0..n => {
            let predicted = intercept + slope * x_data[i];
            observe(addr!("y", i), Normal::new(predicted, sigma).unwrap(), y_data[i])
        });

        pure((intercept, slope, sigma))
    };
    println!("✅ Created Bayesian linear regression model");

Hierarchical Clustering

Hierarchical models implement partial pooling through multi-level parameter structures. The hierarchy enables information sharing across groups while maintaining group-specific effects:

$$\begin{array}{rl}{\mu }_{\text{pop}}& \sim \mathcal{N}(0,{\tau }_{\text{pop}}^2)\\ {\sigma }_{\text{group}}& \sim \text{HalfCauchy}({\sigma }_0)\\ {\mu }_j\mid {\mu }_{\text{pop}},{\sigma }_{\text{group}}& \sim \mathcal{N}({\mu }_{\text{pop}},{\sigma }_{\text{group}}^2)\\ y_{ij}\mid {\mu }_j,{\sigma }_j& \sim \mathcal{N}({\mu }_j,{\sigma }_j^2)\end{array}$$

    // Simplified hierarchy to avoid nested macro issues
    let _multilevel_model = prob!(
        let pop_mean <- sample(addr!("pop_mean"), Normal::new(0.0, 10.0).unwrap());
        let _pop_precision <- sample(addr!("pop_precision"), Gamma::new(2.0, 0.5).unwrap());
        let group_mean <- sample(scoped_addr!("group", "mean", "{}", 0),
                                Normal::new(pop_mean, 1.0).unwrap());
        pure((pop_mean, group_mean))
    );
    println!("✅ Created hierarchical model structure");

State Space Models

Sequential latent variable models:

    // Sequential model with dependencies
    let _sequential_model = prob! {
        let states <- plate!(t in 0..3 => {
            sample(addr!("x", t), Normal::new(0.0, 1.0).unwrap())
                .bind(move |x_t| {
                    observe(addr!("y", t), Normal::new(x_t, 0.5).unwrap(), 1.0 + t as f64)
                        .map(move |_| x_t)
                })
        });

        pure(states)
    };
    println!("✅ Created sequential model with observations");

Multi-Level Hierarchies

Population → Groups → Individuals structure:

    // Simplified hierarchy to avoid nested macro issues
    let _multilevel_model = prob!(
        let pop_mean <- sample(addr!("pop_mean"), Normal::new(0.0, 10.0).unwrap());
        let _pop_precision <- sample(addr!("pop_precision"), Gamma::new(2.0, 0.5).unwrap());
        let group_mean <- sample(scoped_addr!("group", "mean", "{}", 0),
                                Normal::new(pop_mean, 1.0).unwrap());
        pure((pop_mean, group_mean))
    );
    println!("✅ Created hierarchical model structure");

Key Features:

• Partial pooling across hierarchy levels
• Systematic parameter organization
• Natural shrinkage properties
• Scalable to large group structures

Configurable Model Factories

Dynamic model construction:

    // Helper function to create a component model
    fn create_normal_component(name: &str, mean: f64, std: f64) -> Model<f64> {
        sample(addr!(name), Normal::new(mean, std).unwrap())
    }

    // Compose multiple components
    let _composition_model = prob! {
        let param1 <- create_normal_component("param1", 0.0, 1.0);
        let param2 <- create_normal_component("param2", 2.0, 0.5);
        let combined = param1 * param2;
        pure(combined)
    };
    println!("✅ Created composed model with reusable components");

Flexibility Benefits:

• Runtime model configuration
• Conditional model components
• A/B testing different model structures
• Experiment management

Testing Complex Models

Model validation requires systematic testing across multiple dimensions: syntactic correctness, semantic validity, and statistical consistency:

graph TD
    A[Model M] --> B[Syntactic Tests]
    A --> C[Semantic Tests]
    A --> D[Statistical Tests]

    B --> E[Type Checking]
    B --> F[Address Uniqueness]

    C --> G[Trace Validity]
    C --> H[Parameter Bounds]

    D --> I[Prior Predictive]
    D --> J[Posterior Consistency]

    E --> K{All Pass?}
    F --> K
    G --> K
    H --> K
    I --> K
    J --> K

    K -->|Yes| L[Model Validated]
    K -->|No| M[Refinement Required]

Testing Hierarchy:

1. Unit Tests: Individual model components
2. Integration Tests: Model composition correctness
3. Statistical Tests: Distributional properties
4. Performance Tests: Scalability and efficiency

    #[test]
    fn test_model_composition() {
        // Test that models construct without errors
        let _simple = prob! {
            let x <- sample(addr!("test_x"), Normal::new(0.0, 1.0).unwrap());
            pure(x)
        };

        // Test plate notation
        let _plate_model = plate!(i in 0..3 => {
            sample(addr!("plate_test", i), Normal::new(0.0, 1.0).unwrap())
        });

        // Test scoped addresses
        let addr1 = scoped_addr!("test", "param");
        let addr2 = scoped_addr!("test", "param", "{}", 42);

        // Addresses should be different
        assert_ne!(addr1.as_str(), addr2.as_str());
        assert!(addr2.as_str().contains("42"));

        // Test hierarchical model construction
        let _hierarchical = prob! {
            let global <- sample(addr!("global"), Normal::new(0.0, 1.0).unwrap());
            let locals <- plate!(i in 0..2 => {
                sample(scoped_addr!("local", "param", "{}", i),
                       Normal::new(global, 0.1).unwrap())
            });
            pure((global, locals))
        };

        // All models should construct successfully
        // (Actual execution would require handlers)
    }

Common Pitfalls

1. Address Conflicts: Use scoped_addr! for complex models
2. Memory Usage: Large plate operations can create big traces
3. Sequential Dependencies: Explicit state management required
4. Type Inference: Sometimes need explicit type annotations

Performance Considerations

• Plate Size: Very large plates may exceed memory limits
• Nesting Depth: Deep hierarchies increase trace size
• Address Complexity: Simple addresses are more efficient
• Function Composition: Pure functions are optimized away

Next Steps

• Optimization: See Optimizing Performance for efficiency techniques
• Debugging: Check Debugging Models for troubleshooting complex models
• Production: Learn Production Deployment for scaling

Complex models become tractable and maintainable through systematic composition, principled addressing, and mathematical abstraction. Fugue's macro system provides elegant syntactic sugar while preserving the underlying categorical structure that enables powerful inference algorithms and compositional reasoning about probabilistic programs.

Optimizing Performance

Performance optimization in probabilistic programming requires understanding both computational complexity and numerical analysis. This guide explores Fugue's systematic approach to numerical stability and algorithmic efficiency for production-scale probabilistic workloads.

Numerical Stability

Numerical stability in probabilistic computing requires careful analysis of condition numbers and floating-point precision. The log-sum-exp operation is fundamental:

$$\text{LSE}(\mathbf{x})=\log \left(\sum _{i=1}^n{e}^{x_i}\right)=x_{\max }+\log \left(\sum _{i=1}^n{e}^{x_i-x_{\max }}\right)$$

Stability Analysis: Direct computation of $\sum e^{x_i}$ has condition number $\kappa \approx e^{x_{\max }-x_{\min }}$, which becomes ill-conditioned when $x_{\max }-x_{\min }\gg \log ({ϵ}_{\text{machine}})$.

    // Demonstrate stable log-probability computations
    let extreme_log_probs = vec![700.0, 701.0, 699.0, 698.0]; // Would overflow in linear space

    // Safe log-sum-exp prevents overflow
    let log_normalizer = log_sum_exp(&extreme_log_probs);
    let normalized_probs = normalize_log_probs(&extreme_log_probs);

    println!("✅ Stable computation with extreme log-probabilities");
    println!("   - Log normalizer: {:.2}", log_normalizer);
    println!(
        "   - Probabilities sum to: {:.10}",
        normalized_probs.iter().sum::<f64>()
    );

    // Weighted log-sum-exp for importance sampling
    let log_values = vec![-1.0, -2.0, -3.0, -4.0];
    let weights = vec![0.4, 0.3, 0.2, 0.1];
    let weighted_result = weighted_log_sum_exp(&log_values, &weights);

    println!("   - Weighted log-sum-exp: {:.4}", weighted_result);

    // Safe logarithm handling
    let safe_results: Vec<f64> = [1.0, 0.0, -1.0].iter().map(|&x| safe_ln(x)).collect();
    println!("   - Safe ln results: {:?}", safe_results);

Stability Features:

• log_sum_exp prevents overflow in mixture computations
• weighted_log_sum_exp for importance sampling
• safe_ln handles edge cases gracefully
• All operations maintain numerical precision across scales

Vectorized Model Patterns

Structure models for efficient batch processing:

    // Pre-allocate data structures for repeated use
    let observations: Vec<f64> = (0..100).map(|i| i as f64 * 0.1).collect();
    let n = observations.len();

    // Efficient vectorized model
    let vectorized_model = || {
        prob!(
            let mu <- sample(addr!("global_mu"), Normal::new(0.0, 10.0).unwrap());
            let precision <- sample(addr!("precision"), Gamma::new(2.0, 1.0).unwrap());
            let sigma = (1.0 / precision).sqrt();

            // Use plate for efficient vectorized operations
            let _likelihoods <- plate!(i in 0..n => {
                observe(addr!("obs", i), Normal::new(mu, sigma).unwrap(), observations[i])
            });

            pure((mu, sigma))
        )
    };

    let start = Instant::now();
    let (_result, _trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        vectorized_model(),
    );
    let vectorized_time = start.elapsed();

    println!("✅ Optimized vectorized model");
    println!("   - Processed {} observations in {:?}", n, vectorized_time);

Vectorization Strategy:

• Pre-allocate data collections
• Use plate! for independent parallel operations
• Minimize dynamic allocations in hot paths
• Leverage compiler optimizations with static sizing

Performance Monitoring

Systematic performance monitoring requires tracking multiple performance metrics with their theoretical bounds:

$$\begin{array}{rl}\text{Throughput}& =\frac{\text{samples}}{\text{time}}\le \frac{1}{{\tau }_{\min }}\\ \text{Latency}& =\text{time per sample}\ge {\tau }_{\min }\\ \text{Memory Efficiency}& =\frac{\text{useful allocations}}{\text{total allocations}}\to 1\end{array}$$

where ${\tau }_{\min }$ is the theoretical minimum execution time per sample.

    // Monitor trace characteristics for optimization insights
    #[derive(Debug)]
    struct TraceMetrics {
        num_choices: usize,
        log_weight: f64,
        is_valid: bool,
        memory_size_estimate: usize,
    }

    impl TraceMetrics {
        fn from_trace(trace: &Trace) -> Self {
            let num_choices = trace.choices.len();
            let log_weight = trace.total_log_weight();
            let is_valid = log_weight.is_finite();

            // Rough memory estimate (actual implementation would be more precise)
            let memory_size_estimate = num_choices * 64; // Rough bytes per choice

            Self {
                num_choices,
                log_weight,
                is_valid,
                memory_size_estimate,
            }
        }
    }

    // Example: Monitor a complex model's performance
    let complex_model = || {
        prob!(
            let components <- plate!(c in 0..5 => {
                sample(addr!("weight", c), Gamma::new(1.0, 1.0).unwrap())
                    .bind(move |weight| {
                        sample(addr!("mu", c), Normal::new(0.0, 2.0).unwrap())
                            .map(move |mu| (weight, mu))
                    })
            });

            let selector <- sample(addr!("selector"),
                                  Categorical::new(vec![0.2, 0.2, 0.2, 0.2, 0.2]).unwrap());

            pure((components, selector))
        )
    };

    let (_result, trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        complex_model(),
    );

    let metrics = TraceMetrics::from_trace(&trace);
    println!("✅ Performance monitoring active");
    println!("   - Trace choices: {}", metrics.num_choices);
    println!("   - Log weight: {:.2}", metrics.log_weight);
    println!("   - Valid: {}", metrics.is_valid);
    println!(
        "   - Memory estimate: {} bytes",
        metrics.memory_size_estimate
    );

Monitoring Approach:

• Collect trace characteristics for optimization insights
• Track memory usage patterns
• Validate numerical stability
• Profile execution bottlenecks

Numerical Precision Testing

Validate stability across different computational scales:

    // Test numerical stability across different scales
    let test_scales = vec![1e-10, 1e-5, 1.0, 1e5, 1e10];

    for &scale in &test_scales {
        let scale: f64 = scale;
        let log_vals = vec![scale.ln() + 1.0, scale.ln() + 2.0, scale.ln() + 0.5];

        let stable_sum = log_sum_exp(&log_vals);
        let log1p_result = log1p_exp(scale.ln());

        println!(
            "   Scale {:.0e}: log_sum_exp={:.4}, log1p_exp={:.4}",
            scale, stable_sum, log1p_result
        );
    }
    println!("✅ Numerical stability verified across scales");

Testing Strategy:

• Verify stability across extreme value ranges
• Test edge cases and boundary conditions
• Validate consistency of numerical operations
• Profile precision vs. performance trade-offs

Performance Testing

Implement systematic performance validation:

    #[test]
    fn test_numerical_stability() {
        // Test log_sum_exp with extreme values
        let extreme_vals = vec![700.0, 701.0, 699.0];
        let result = log_sum_exp(&extreme_vals);
        assert!(
            result.is_finite(),
            "log_sum_exp should handle extreme values"
        );

        // Test normalization
        let normalized = normalize_log_probs(&extreme_vals);
        let sum: f64 = normalized.iter().sum();
        assert!(
            (sum - 1.0).abs() < 1e-10,
            "Normalized probabilities should sum to 1"
        );

        // Test weighted computation
        let weights = vec![0.5, 0.3, 0.2];
        let weighted_result = weighted_log_sum_exp(&extreme_vals, &weights);
        assert!(
            weighted_result.is_finite(),
            "Weighted log_sum_exp should be finite"
        );
    }

Testing Framework:

• Numerical stability regression tests
• End-to-end inference benchmarking (cargo bench --bench f_perf)

Production Deployment

Numerical Strategies

• Always use log-space for probability computations
• Validate extreme value handling in testing
• Monitor for numerical instabilities in production
• Use stable algorithms for critical computations

Monitoring and Alerting

• Track inference latency and memory usage
• Alert on numerical instabilities or performance degradation
• Profile hot paths for optimization opportunities

Common Performance Patterns

1. Log Always: Work in log-space for numerical stability
2. Batch Everything: Amortize costs across multiple samples
3. Monitor Continuously: Track performance metrics in production
4. Test Extremes: Validate stability with extreme values

These optimization strategies enable Fugue to handle production-scale probabilistic programming workloads with consistent performance and numerical reliability.

Debugging Models

Debugging probabilistic models presents unique challenges due to their stochastic nature and high-dimensional parameter spaces. Unlike deterministic programs, probabilistic models require statistical validation, convergence analysis, and distributional testing. This guide establishes a systematic methodology for probabilistic model debugging using Fugue's comprehensive diagnostic framework.

Trace Inspection and Analysis

Execution traces form the foundation of probabilistic model debugging. Each trace $\mathcal{T}$ contains a complete record of the program's stochastic execution:

$$\mathcal{T}=\{(a_i,v_i,w_i){\}}_{i=1}^n$$

where $a_i$ is the address, $v_i$ is the sampled value, and $w_i$ is the log-weight contribution.

graph TD
    subgraph "Trace Analysis Workflow"
        A[Execution Trace T] --> B{Finite Log-Weight?}
        B -->|No| C[Numerical Instability]
        B -->|Yes| D[Choice Analysis]
        D --> E[Parameter Extraction]
        E --> F[Statistical Validation]
        F --> G{Passes Tests?}
        G -->|No| H[Model Refinement]
        G -->|Yes| I[Model Validated]
        C --> J[Debug Constraints]
        H --> A
        J --> A
    end

Mathematical Properties: A valid trace must satisfy the weight consistency condition: $$\log P(\mathcal{T})=\sum _{i=1}^n{w}_i<\infty$$

    // Execute a model and examine its trace structure
    let mut rng = thread_rng();

    let diagnostic_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap());
            let sigma <- sample(addr!("sigma"), Gamma::new(2.0, 1.0).unwrap());
            observe(addr!("obs1"), Normal::new(mu, sigma).unwrap(), 1.5);
            observe(addr!("obs2"), Normal::new(mu, sigma).unwrap(), 1.2);
            factor(if mu.abs() < 3.0 { 0.0 } else { f64::NEG_INFINITY });
            pure((mu, sigma))
        )
    };

    let ((mu_val, sigma_val), trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        diagnostic_model(),
    );

    println!("✅ Model execution complete");
    println!("   - Result: mu = {:.3}, sigma = {:.3}", mu_val, sigma_val);
    println!("   - Choices recorded: {}", trace.choices.len());
    println!("   - Prior log-weight: {:.6}", trace.log_prior);
    println!("   - Likelihood log-weight: {:.6}", trace.log_likelihood);
    println!("   - Factor log-weight: {:.6}", trace.log_factors);
    println!("   - Total log-weight: {:.6}", trace.total_log_weight());

    // Per-choice breakdown
    println!("   - Choice breakdown:");
    for (addr, choice) in &trace.choices {
        println!(
            "     {}: {:?} (logp: {:.6})",
            addr, choice.value, choice.logp
        );
    }

Key Debugging Insights:

• Choice count reveals model complexity and structure
• Log-weight decomposition identifies prior vs. likelihood vs. factor issues
• Per-choice analysis shows individual parameter contributions
• Finite log-weights indicate valid model execution

Type-Safe Value Access

Fugue provides robust access patterns that handle type mismatches gracefully:

    // Safe access patterns that handle type mismatches gracefully

    // Option-based access (returns None on mismatch)
    match trace.get_f64(&addr!("mu")) {
        Some(mu) => println!("✅ Retrieved mu = {:.3}", mu),
        None => println!("❌ Failed to get mu as f64"),
    }

    // Result-based access (returns detailed error info)
    match trace.get_f64_result(&addr!("sigma")) {
        Ok(sigma) => println!("✅ Retrieved sigma = {:.3}", sigma),
        Err(e) => println!("❌ Error getting sigma: {}", e),
    }

    // Handle missing addresses
    match trace.get_f64_result(&addr!("missing_param")) {
        Ok(_) => unreachable!(),
        Err(e) => println!("✅ Correctly caught missing address: {}", e),
    }

    // Handle type mismatches
    match trace.get_bool_result(&addr!("mu")) {
        Ok(_) => unreachable!(),
        Err(e) => println!("✅ Correctly caught type mismatch: {}", e),
    }

    // Iterate through all choices for debugging
    println!("   - All choices and their types:");
    for (addr, choice) in &trace.choices {
        let type_info = match &choice.value {
            ChoiceValue::F64(_) => "f64",
            ChoiceValue::Bool(_) => "bool",
            ChoiceValue::U64(_) => "u64",
            ChoiceValue::I64(_) => "i64",
            ChoiceValue::Usize(_) => "usize",
        };
        println!("     {} ({}): {:?}", addr, type_info, choice.value);
    }

Error Handling Strategies:

• Use get_*_result() for detailed error information
• Use get_*() for simple None-handling
• Always check for missing addresses before assuming success
• Iterate through all choices to understand model structure

Model Validation and Testing

Systematic validation ensures your model behaves as expected:

    // Test a simple conjugate model against analytical solution
    let conjugate_model = || {
        prob!(
            let theta <- sample(addr!("theta"), Beta::new(1.0, 1.0).unwrap());
            observe(addr!("successes"), Binomial::new(10, theta).unwrap(), 7u64);
            pure(theta)
        )
    };

    // Run a few samples to test basic functionality
    let mut theta_samples = Vec::new();
    for _ in 0..20 {
        let (theta, test_trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            conjugate_model(),
        );

        // Validate trace structure
        assert!(test_trace.choices.contains_key(&addr!("theta")));
        assert!(
            test_trace.total_log_weight().is_finite(),
            "Trace should have finite log-weight"
        );
        assert!(
            test_trace.log_likelihood.is_finite(),
            "Likelihood should be finite"
        );

        theta_samples.push(theta);
    }

    // Basic statistical checks
    let sample_mean = theta_samples.iter().sum::<f64>() / theta_samples.len() as f64;
    println!("✅ Validation tests passed");
    println!("   - Generated {} samples", theta_samples.len());
    println!(
        "   - Sample mean: {:.3} (expected ~0.7 for Beta-Binomial)",
        sample_mean
    );
    println!("   - All traces had finite log-weights");

Validation Best Practices:

• Test against known analytical solutions
• Verify all traces have finite log-weights
• Check basic statistical properties (means, variances)
• Test edge cases and boundary conditions

Safe vs Strict Error Handling

Fugue provides both strict (fail-fast) and safe (error-resilient) execution modes:

    // Create a trace with known structure for replay testing
    let mut base_trace = Trace::default();
    base_trace.insert_choice(addr!("param"), ChoiceValue::F64(1.5), -0.5);

    let test_model = || sample(addr!("param"), Normal::new(0.0, 1.0).unwrap());

    // Strict replay - will panic on mismatch (commented out for safety)
    // let strict_replay = ReplayHandler { base_trace: &base_trace };
    // let (strict_result, strict_trace) = runtime::handler::run(strict_replay, test_model());

    // Safe replay - handles errors gracefully
    let safe_replay = SafeReplayHandler {
        rng: &mut rng,
        base: base_trace.clone(),
        trace: Trace::default(),
        warn_on_mismatch: true,
    };
    let (safe_result, safe_trace) = runtime::handler::run(safe_replay, test_model());

    println!("✅ Safe replay succeeded");
    println!("   - Result: {:.3}", safe_result);
    println!(
        "   - Retrieved value: {:?}",
        safe_trace.get_f64(&addr!("param"))
    );

    // Test scoring with safe handler
    let safe_score = SafeScoreGivenTrace {
        base: base_trace,
        trace: Trace::default(),
        warn_on_error: false,
    };
    let (_, score_trace) = runtime::handler::run(safe_score, test_model());

    println!(
        "   - Score trace log-weight: {:.3}",
        score_trace.total_log_weight()
    );

When to Use Each:

• Strict handlers (ReplayHandler, ScoreGivenTrace): Development and testing
• Safe handlers (SafeReplayHandler, SafeScoreGivenTrace): Production systems
• Safe handlers log warnings instead of panicking on mismatches

MCMC Diagnostics

Markov Chain Monte Carlo convergence assessment requires statistical hypothesis testing and diagnostic metrics. The fundamental question is whether the chain has reached its stationary distribution $\pi (\theta )$.

Split-R-hat Diagnostic

The potential scale reduction factor $\hat{R}$ compares within-chain and between-chain variance:

$$\hat{R}=\sqrt{\frac{\hat{V}}{W}}$$

where:

• $W=\frac{1}{m}{\sum }_{j=1}^m{s}_j^2$ (within-chain variance)
• $B=\frac{n}{m-1}{\sum }_{j=1}^m({\bar{\theta }}_{j\cdot }-{\bar{\theta }}_{\cdot \cdot }{)}^2$ (between-chain variance)
• $\hat{V}=\frac{n-1}{n}W+\frac{1}{n}B$ (marginal posterior variance estimate)

Effective Sample Size

The effective sample size accounts for autocorrelation in MCMC samples:

$$\text{ESS}=\frac{mn}{1+2{\sum }_{t=1}^{\infty }{\rho }_t}$$

where ${\rho }_t$ is the lag-$t$ autocorrelation and $mn$ is the total number of samples.

    // Generate simple MCMC chains for diagnostic testing
    let mcmc_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());
            observe(addr!("y"), Normal::new(mu, 0.5).unwrap(), 1.0);
            pure(mu)
        )
    };

    // Generate two short chains for R-hat calculation
    let n_samples = 50;
    let n_warmup = 10;

    let mut chain1_samples = Vec::new();
    let mut chain2_samples = Vec::new();

    // Chain 1
    let mut rng1 = rand::rngs::StdRng::seed_from_u64(42);
    let chain1 = adaptive_mcmc_chain(&mut rng1, mcmc_model, n_samples, n_warmup);
    for (_, trace) in &chain1 {
        if let Some(mu) = trace.get_f64(&addr!("mu")) {
            chain1_samples.push(mu);
        }
    }

    // Chain 2
    let mut rng2 = rand::rngs::StdRng::seed_from_u64(123);
    let chain2 = adaptive_mcmc_chain(&mut rng2, mcmc_model, n_samples, n_warmup);
    for (_, trace) in &chain2 {
        if let Some(mu) = trace.get_f64(&addr!("mu")) {
            chain2_samples.push(mu);
        }
    }

    // Compute diagnostics
    if !chain1_samples.is_empty() && !chain2_samples.is_empty() {
        // Extract traces for R-hat calculation
        let chain1_traces: Vec<Trace> = chain1.into_iter().map(|(_, trace)| trace).collect();
        let chain2_traces: Vec<Trace> = chain2.into_iter().map(|(_, trace)| trace).collect();
        let r_hat = r_hat_f64(&[chain1_traces, chain2_traces], &addr!("mu"));
        let ess1 = effective_sample_size_mcmc(&chain1_samples);
        let ess2 = effective_sample_size_mcmc(&chain2_samples);

        println!("✅ MCMC diagnostics computed");
        println!(
            "   - Chain 1: {} samples, ESS = {:.1}",
            chain1_samples.len(),
            ess1
        );
        println!(
            "   - Chain 2: {} samples, ESS = {:.1}",
            chain2_samples.len(),
            ess2
        );
        println!("   - R-hat: {:.4} (< 1.1 indicates convergence)", r_hat);

        if r_hat < 1.1 {
            println!("   - ✅ Chains appear to have converged");
        } else {
            println!("   - ⚠️  Chains may not have converged - run longer");
        }
    }

effective_sample_size_mcmc (in fugue::inference::mcmc_utils, used above) computes ESS for a single chain. Fugue also provides effective_sample_size_multichain, which pools autocorrelation and between-chain variance across all chains at once — prefer it over averaging single-chain ESS values when you have more than one chain.

Convergence Indicators:

• Split-R-hat < 1.1: Chains have converged
• High ESS: Efficient sampling without excessive correlation
• Multiple chains: Essential for reliable convergence assessment
• Visual inspection: Always examine trace plots when possible

Model Structure Analysis

Model structure analysis reveals the computational graph and parameter dependencies. This analysis is crucial for understanding model complexity and identifying potential issues:

graph TD
    subgraph "Model Structure Hierarchy"
        A[Model M] --> B[Parameter Groups]
        B --> C1[Hyperpriors θ₁]
        B --> C2[Primary Parameters θ₂] 
        B --> C3[Observations y]
        C1 --> D1[Constraint Analysis]
        C2 --> D2[Dependency Graph]
        C3 --> D3[Likelihood Terms]
        D1 --> E[Structure Validation]
        D2 --> E
        D3 --> E
    end

Structural Invariants to validate:

1. Address Uniqueness: $|\{a_i\}|=n$ (no collisions)
2. Parameter Hierarchy: $\forall i,j:a_i⪯a_j⟹\text{dependency}(i,j)$
3. Choice Count Consistency: Expected vs. actual parameter count
4. Type Safety: Each address maps to consistent value types

    // Create a complex model to demonstrate structure analysis
    let complex_model = || {
        prob!(
            // Hierarchical structure
            let global_scale <- sample(addr!("global_scale"), Gamma::new(2.0, 1.0).unwrap());

            let group_params <- plate!(g in 0..3 => {
                sample(addr!("group_mean", g), Normal::new(0.0, global_scale).unwrap())
                    .bind(move |mean| {
                        sample(addr!("group_precision", g), Gamma::new(2.0, 1.0).unwrap())
                            .map(move |prec| (mean, prec))
                    })
            });

            // Individual observations (simplified to avoid move issues)
            let observations = [1.2, 1.5, 0.8];
            let likelihoods <- plate!(i in 0..observations.len() => {
                // Use fixed parameters for demonstration
                observe(addr!("obs", i), Normal::new(0.0, 1.0).unwrap(), observations[i])
            });

            pure((global_scale, group_params, likelihoods))
        )
    };

    let (_result, complex_trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        complex_model(),
    );

    // Analyze model structure
    let mut address_analysis = BTreeMap::new();
    for (addr, choice) in &complex_trace.choices {
        let addr_str = addr.as_str().to_string();
        let category = if addr_str.contains("global") {
            "Global Parameters"
        } else if addr_str.contains("group") {
            "Group Parameters"
        } else if addr_str.contains("obs") {
            "Observations"
        } else {
            "Other"
        };

        address_analysis
            .entry(category)
            .or_insert(Vec::new())
            .push((addr_str, choice.logp));
    }

    println!("✅ Complex model structure analysis");
    println!("   - Total choices: {}", complex_trace.choices.len());
    println!("   - Address structure:");
    for (category, addresses) in address_analysis {
        println!("     {}: {} choices", category, addresses.len());
        for (addr, logp) in addresses.iter().take(3) {
            // Show first 3
            println!("       {} (logp: {:.3})", addr, logp);
        }
        if addresses.len() > 3 {
            println!("       ... and {} more", addresses.len() - 3);
        }
    }

Structure Analysis Benefits:

• Understand parameter organization and hierarchies
• Detect unexpected address patterns
• Verify choice counts match model expectations
• Identify bottlenecks in complex models

Performance Diagnostics

Monitor computational efficiency and identify bottlenecks:

    use std::time::Instant;

    // Benchmark model execution and trace construction
    let benchmark_model = || {
        prob!(
            let params <- plate!(i in 0..100 => {
                sample(addr!("param", i), Normal::new(0.0, 1.0).unwrap())
            });
            pure(params)
        )
    };

    let start = Instant::now();
    let (_, bench_trace) = runtime::handler::run(
        PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        },
        benchmark_model(),
    );
    let execution_time = start.elapsed();

    // Analyze trace characteristics
    let choice_count = bench_trace.choices.len();
    let memory_estimate = choice_count * 64; // Rough estimate
    let log_weight_is_finite = bench_trace.total_log_weight().is_finite();

    println!("✅ Performance diagnostics");
    println!("   - Execution time: {:?}", execution_time);
    println!("   - Choices created: {}", choice_count);
    println!("   - Memory estimate: ~{} bytes", memory_estimate);
    println!("   - Log-weight valid: {}", log_weight_is_finite);

    // Check for potential issues
    if choice_count == 0 {
        println!("   - ⚠️  No choices recorded - possible model issue");
    }
    if !log_weight_is_finite {
        println!("   - ⚠️  Invalid log-weight - check factors and observations");
    }
    if execution_time.as_millis() > 100 {
        println!("   - ⚠️  Slow execution - consider optimization");
    }

Performance Warning Signs:

• Zero choices recorded (model execution failure)
• Infinite log-weights (constraint violations)
• Excessive execution time (optimization needed)
• Large memory footprint (consider streaming approaches)

Common Debugging Patterns

Systematic debugging follows a hierarchical validation strategy from basic correctness to statistical validity:

graph TD
    subgraph "Debugging Methodology"
        A[Model Implementation] --> B{Syntax Valid?}
        B -->|No| C[Fix Code Structure]
        B -->|Yes| D{Types Consistent?}
        D -->|No| E[Fix Type Errors]
        D -->|Yes| F{Finite Log-Weights?}
        F -->|No| G[Fix Constraints]
        F -->|Yes| H{Statistical Properties?}
        H -->|No| I[Validate Distributions]
        H -->|Yes| J{Convergence?}
        J -->|No| K[Tune Inference]
        J -->|Yes| L[Model Validated]
        
        C --> A
        E --> A
        G --> A
        I --> A
        K --> A
    end

Debug Level Hierarchy:

1. Syntactic: Code compiles and types check
2. Semantic: Model executes without runtime errors
3. Numerical: Computations remain stable and finite
4. Statistical: Results match theoretical expectations
5. Convergence: Inference algorithms reach stationarity

    // Pattern 1: Systematic model testing
    fn test_model_basic_properties<F, T>(
        model_fn: F,
        expected_choice_count: usize,
        description: &str,
    ) where
        F: Fn() -> Model<T>,
    {
        let mut rng = thread_rng();
        let (_, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            model_fn(),
        );

        println!("Testing {}", description);

        // Basic trace validity
        assert!(
            trace.total_log_weight().is_finite(),
            "Log-weight should be finite"
        );
        assert_eq!(
            trace.choices.len(),
            expected_choice_count,
            "Choice count mismatch"
        );

        // Check for common issues
        if trace.log_prior.is_infinite() {
            println!("  - ⚠️  Infinite prior - check parameter ranges");
        }
        if trace.log_likelihood.is_infinite() {
            println!("  - ⚠️  Infinite likelihood - check observations");
        }
        if trace.log_factors.is_infinite() {
            println!("  - ⚠️  Infinite factors - check constraint satisfaction");
        }

        println!("  - ✅ {} passed basic tests", description);
    }

    // Test simple models
    test_model_basic_properties(
        || sample(addr!("x"), Normal::new(0.0, 1.0).unwrap()),
        1,
        "Simple normal sampling",
    );

    test_model_basic_properties(
        || {
            prob!(
                let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
                observe(addr!("y"), Normal::new(x, 0.5).unwrap(), 1.0);
                pure(x)
            )
        },
        1,
        "Normal model with observation",
    );

    // Pattern 2: Address collision detection
    fn check_address_collisions(trace: &Trace) -> Vec<String> {
        let mut collisions = Vec::new();
        let addresses: Vec<&str> = trace.choices.keys().map(|addr| addr.as_str()).collect();

        for (i, addr1) in addresses.iter().enumerate() {
            for addr2 in addresses.iter().skip(i + 1) {
                if addr1 == addr2 {
                    collisions.push(format!("Duplicate address: {}", addr1));
                }
            }
        }
        collisions
    }

    let test_trace = complex_trace; // Use trace from earlier
    let collisions = check_address_collisions(&test_trace);
    if collisions.is_empty() {
        println!("  - ✅ No address collisions detected");
    } else {
        for collision in collisions {
            println!("  - ⚠️  {}", collision);
        }
    }

    println!("✅ Debugging patterns demonstration complete");

Debugging Workflow:

1. Start Simple: Test individual components before complex composition
2. Validate Incrementally: Add complexity one piece at a time
3. Check Address Uniqueness: Prevent parameter collision bugs
4. Monitor Log-Weights: Track prior, likelihood, and factor contributions
5. Use Systematic Testing: Automated validation for all model components

Testing Framework Integration

Embed debugging checks in your test suite:

    #[test]
    fn test_trace_inspection_patterns() {
        let mut rng = thread_rng();

        let model = prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            let y <- sample(addr!("y"), Beta::new(1.0, 1.0).unwrap());
            observe(addr!("obs"), Normal::new(x, 0.1).unwrap(), 1.5);
            pure((x, y))
        );

        let (_result, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            model,
        );

        // Basic trace properties
        assert_eq!(trace.choices.len(), 2); // x and y samples
        assert!(trace.total_log_weight().is_finite());
        assert!(trace.log_likelihood.is_finite());

        // Type-safe access
        assert!(trace.get_f64(&addr!("x")).is_some());
        assert!(trace.get_f64(&addr!("y")).is_some());
        assert!(trace.get_bool(&addr!("x")).is_none()); // Type mismatch

        // Result access patterns
        assert!(trace.get_f64_result(&addr!("x")).is_ok());
        assert!(trace.get_f64_result(&addr!("missing")).is_err());
    }

    #[test]
    fn test_safe_vs_strict_handlers() {
        let mut rng = thread_rng();

        // Create base trace
        let mut base_trace = Trace::default();
        base_trace.insert_choice(addr!("param"), ChoiceValue::F64(2.5), -1.0);

        let model = sample(addr!("param"), Normal::new(0.0, 1.0).unwrap());

        // Safe replay should work
        let safe_handler = SafeReplayHandler {
            rng: &mut rng,
            base: base_trace,
            trace: Trace::default(),
            warn_on_mismatch: false,
        };
        let (result, trace) = runtime::handler::run(safe_handler, model);

        assert_eq!(result, 2.5);
        assert_eq!(trace.get_f64(&addr!("param")), Some(2.5));
    }

    #[test]
    fn test_model_structure_analysis() {
        let mut rng = thread_rng();

        let hierarchical_model = || {
            prob!(
                let global <- sample(addr!("global"), Normal::new(0.0, 1.0).unwrap());
                let locals <- plate!(i in 0..3 => {
                    sample(addr!("local", i), Normal::new(global, 0.1).unwrap())
                });
                pure((global, locals))
            )
        };

        let (_, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            hierarchical_model(),
        );

        // Should have global + 3 local parameters
        assert_eq!(trace.choices.len(), 4);

        // Check address structure
        assert!(trace.choices.contains_key(&addr!("global")));
        assert!(trace.choices.contains_key(&addr!("local", 0)));
        assert!(trace.choices.contains_key(&addr!("local", 1)));
        assert!(trace.choices.contains_key(&addr!("local", 2)));
    }

    #[test]
    fn test_performance_diagnostics() {
        use std::time::Instant;
        let mut rng = thread_rng();

        let large_model = || {
            plate!(i in 0..50 => {
                sample(addr!("x", i), Normal::new(0.0, 1.0).unwrap())
            })
        };

        let start = Instant::now();
        let (_, trace) = runtime::handler::run(
            PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            },
            large_model(),
        );
        let duration = start.elapsed();

        assert_eq!(trace.choices.len(), 50);
        assert!(trace.total_log_weight().is_finite());

        // Performance should be reasonable
        assert!(duration.as_millis() < 1000, "Model execution too slow");
    }

Testing Strategy:

• Unit tests for individual model components
• Integration tests for complete workflows
• Performance regression tests
• Statistical validation against known results

Common Issues and Solutions

Issue: Infinite Log-Weights

Symptoms: trace.total_log_weight().is_infinite()

Causes:

• Factor statements with impossible constraints
• Parameters outside valid ranges
• Numerical overflow in likelihood computations

Solutions:

• Check factor conditions carefully
• Validate parameter ranges in constructors
• Use log-space computations for numerical stability

Issue: Missing or Wrong Parameter Values

Symptoms: get_*() returns None or wrong types

Causes:

• Address typos or inconsistencies
• Model structure doesn't match expectations
• Type mismatches in trace replay

Solutions:

• Use consistent address naming conventions
• Print all addresses for verification
• Use safe handlers for production resilience

Issue: Poor MCMC Convergence

Symptoms: High R-hat values, low ESS

Causes:

• Inappropriate step sizes
• Poor model parameterization
• Insufficient warm-up periods

Solutions:

• Increase warm-up iterations
• Reparameterize for better geometry
• Use adaptive algorithms with proper tuning

Issue: Slow Model Execution

Symptoms: High execution times, memory usage

Causes:

• Inefficient model structure
• Excessive address creation
• Large trace construction overhead

Solutions:

• Use plate! for vectorized operations
• Pre-allocate data structures when possible
• Profile with performance diagnostics

Best Practices Summary

1. Debug Incrementally: Start simple and add complexity systematically
2. Use All Tools: Combine trace inspection, validation, and diagnostics
3. Test Edge Cases: Verify behavior at parameter boundaries
4. Monitor Performance: Track execution time and memory usage
5. Validate Statistically: Compare against known theoretical results
6. Handle Errors Gracefully: Use safe handlers in production
7. Document Assumptions: Clear model specifications aid debugging

Effective debugging transforms probabilistic programming from guesswork into systematic model development. Fugue's comprehensive debugging toolkit enables confident deployment of complex probabilistic systems.

Custom Handlers

Fugue's handler system is grounded in algebraic effect theory, providing a principled approach to effect interpretation and computational extension. Custom handlers enable specialized execution strategies, monitoring systems, and novel inference algorithms through systematic effect handling and handler composition.

Understanding the Handler Trait

The Handler trait provides the algebraic signature for probabilistic effects. Each method represents an effect operation with its semantic interpretation:

graph TD
    subgraph "Handler Architecture"
        A[Effect E] --> B{Effect Type}
        B -->|sample| C[on_sample_T]
        B -->|observe| D[on_observe_T]  
        B -->|factor| E[on_factor]
        C --> F[Handler State H]
        D --> F
        E --> F
        F --> G[Updated State H']
        G --> H[Continue Execution]
    end

Effect Algebra: Each handler interprets the probabilistic effect signature:

$$\begin{array}{rl}\text{sample}& :\text{Dist}[T]\to T\\ \text{observe}& :\text{Dist}[T]\times T\to \text{Unit}\\ \text{factor}& :\mathbb{R}\to \text{Unit}\end{array}$$

where the carrier type varies by handler implementation.

/// Simple handler that just samples from priors (similar to PriorHandler)
struct BasicHandler<R: Rng> {
    rng: R,
    trace: Trace,
}

impl<R: Rng> Handler for BasicHandler<R> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = dist.sample(&mut self.rng);
        let log_prob = dist.log_prob(&value);

        // Store in trace
        self.trace.log_prior += log_prob;
        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::F64(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = dist.sample(&mut self.rng);
        let log_prob = dist.log_prob(&value);

        self.trace.log_prior += log_prob;
        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Bool(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        let value = dist.sample(&mut self.rng);
        let log_prob = dist.log_prob(&value);

        self.trace.log_prior += log_prob;
        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::U64(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = dist.sample(&mut self.rng);
        let log_prob = dist.log_prob(&value);

        self.trace.log_prior += log_prob;
        self.trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Usize(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_observe_f64(&mut self, _addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_bool(&mut self, _addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_u64(&mut self, _addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_usize(&mut self, _addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.trace.log_factors += log_weight;
    }

    fn finish(self) -> Trace {
        self.trace
    }
}

Handler Responsibilities:

• Type-specific sampling: Handle f64, bool, u64, and usize distributions appropriately
• Observation handling: Process observed values and update likelihood components
• Factor management: Accumulate constraint and penalty terms
• Trace construction: Build execution traces with choices and log-weights
• Resource cleanup: Properly finalize and return traces

Decorator Pattern for Handler Composition

The decorator pattern implements handler composition through effect forwarding with computational augmentation. This pattern follows the mathematical principle of function composition:

$$(f\circ g)(x)=f(g(x))$$

Applied to handlers: $h_{\text{decorated}}=h_{\text{decorator}}\circ h_{\text{base}}$

graph LR
    subgraph "Handler Composition Chain"
        A[Effect] --> B[Decorator₁]
        B --> C[Decorator₂]
        C --> D[Base Handler]
        D --> E[Result]

        B -.->|"Log, Monitor"| F[Side Effects]
        C -.->|"Transform, Filter"| G[Modifications]
    end

Compositional Properties:

• Associativity: $(h_1\circ h_2)\circ h_3=h_1\circ (h_2\circ h_3)$
• Identity: $\text{id}\circ h=h\circ \text{id}=h$
• Effect Preservation: Core semantics remain unchanged

/// Handler decorator that logs all operations
struct LoggingHandler<H: Handler> {
    inner: H,
    log: Vec<String>,
    verbose: bool,
}

impl<H: Handler> LoggingHandler<H> {
    fn new(inner: H, verbose: bool) -> Self {
        Self {
            inner,
            log: Vec::new(),
            verbose,
        }
    }

    fn log_operation(&mut self, operation: String) {
        if self.verbose {
            println!("LOG: {}", operation);
        }
        self.log.push(operation);
    }
}

impl<H: Handler> Handler for LoggingHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = self.inner.on_sample_f64(addr, dist);
        self.log_operation(format!("Sample f64 at {}: {:.3}", addr, value));
        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = self.inner.on_sample_bool(addr, dist);
        self.log_operation(format!("Sample bool at {}: {}", addr, value));
        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        let value = self.inner.on_sample_u64(addr, dist);
        self.log_operation(format!("Sample u64 at {}: {}", addr, value));
        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = self.inner.on_sample_usize(addr, dist);
        self.log_operation(format!("Sample usize at {}: {}", addr, value));
        value
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.log_operation(format!("Observe f64 at {}: {:.3}", addr, value));
        self.inner.on_observe_f64(addr, dist, value);
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.log_operation(format!("Observe bool at {}: {}", addr, value));
        self.inner.on_observe_bool(addr, dist, value);
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.log_operation(format!("Observe u64 at {}: {}", addr, value));
        self.inner.on_observe_u64(addr, dist, value);
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.log_operation(format!("Observe usize at {}: {}", addr, value));
        self.inner.on_observe_usize(addr, dist, value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.log_operation(format!("Factor: {:.3}", log_weight));
        self.inner.on_factor(log_weight);
    }

    fn finish(self) -> Trace {
        let trace = self.inner.finish();
        println!("✅ Logged {} operations total", self.log.len());
        trace
    }
}

Decorator Benefits:

• Non-invasive functionality addition
• Composable and reusable components
• Separation of concerns between core logic and cross-cutting features
• Easy to enable/disable features dynamically

Stateful Handlers for Analytics

Handlers can maintain state to accumulate statistics and monitor model behavior:

/// Handler that accumulates statistics about model execution
#[derive(Debug)]
struct ExecutionStats {
    sample_counts: HashMap<String, u32>, // Type -> count
    observe_counts: HashMap<String, u32>,
    factor_count: u32,
    total_log_weight: f64,
    parameter_ranges: HashMap<String, (f64, f64)>, // Address -> (min, max) for f64 params
}

impl Default for ExecutionStats {
    fn default() -> Self {
        Self {
            sample_counts: HashMap::new(),
            observe_counts: HashMap::new(),
            factor_count: 0,
            total_log_weight: 0.0,
            parameter_ranges: HashMap::new(),
        }
    }
}

struct StatisticsHandler<H: Handler> {
    inner: H,
    stats: ExecutionStats,
}

impl<H: Handler> StatisticsHandler<H> {
    fn new(inner: H) -> Self {
        Self {
            inner,
            stats: ExecutionStats::default(),
        }
    }

    fn update_f64_range(&mut self, addr: &Address, value: f64) {
        let key = addr.as_str().to_string();
        self.stats
            .parameter_ranges
            .entry(key)
            .and_modify(|(min, max)| {
                *min = min.min(value);
                *max = max.max(value);
            })
            .or_insert((value, value));
    }
}

impl<H: Handler> Handler for StatisticsHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = self.inner.on_sample_f64(addr, dist);
        *self
            .stats
            .sample_counts
            .entry("f64".to_string())
            .or_insert(0) += 1;
        self.update_f64_range(addr, value);
        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = self.inner.on_sample_bool(addr, dist);
        *self
            .stats
            .sample_counts
            .entry("bool".to_string())
            .or_insert(0) += 1;
        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        let value = self.inner.on_sample_u64(addr, dist);
        *self
            .stats
            .sample_counts
            .entry("u64".to_string())
            .or_insert(0) += 1;
        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = self.inner.on_sample_usize(addr, dist);
        *self
            .stats
            .sample_counts
            .entry("usize".to_string())
            .or_insert(0) += 1;
        value
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        *self
            .stats
            .observe_counts
            .entry("f64".to_string())
            .or_insert(0) += 1;
        self.inner.on_observe_f64(addr, dist, value);
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        *self
            .stats
            .observe_counts
            .entry("bool".to_string())
            .or_insert(0) += 1;
        self.inner.on_observe_bool(addr, dist, value);
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        *self
            .stats
            .observe_counts
            .entry("u64".to_string())
            .or_insert(0) += 1;
        self.inner.on_observe_u64(addr, dist, value);
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        *self
            .stats
            .observe_counts
            .entry("usize".to_string())
            .or_insert(0) += 1;
        self.inner.on_observe_usize(addr, dist, value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.stats.factor_count += 1;
        self.stats.total_log_weight += log_weight;
        self.inner.on_factor(log_weight);
    }

    fn finish(self) -> Trace {
        println!("✅ Execution Statistics:");
        println!("   - Samples by type: {:?}", self.stats.sample_counts);
        println!("   - Observations by type: {:?}", self.stats.observe_counts);
        println!("   - Factor operations: {}", self.stats.factor_count);
        println!("   - Parameter ranges:");
        for (addr, (min, max)) in &self.stats.parameter_ranges {
            println!("     {}: [{:.3}, {:.3}]", addr, min, max);
        }
        self.inner.finish()
    }
}

Analytics Applications:

• Model complexity analysis (parameter counts by type)
• Execution profiling and bottleneck identification
• Parameter range monitoring for numerical stability
• Distribution usage patterns for optimization

Conditional and Filtering Handlers

Implement business logic and constraints through conditional handling:

/// Handler that filters/modifies values based on conditions
struct FilteringHandler<H: Handler> {
    inner: H,
    f64_clamp_range: Option<(f64, f64)>,
    bool_flip_probability: f64,
    rng: rand::rngs::ThreadRng,
}

impl<H: Handler> FilteringHandler<H> {
    fn new(inner: H, f64_clamp_range: Option<(f64, f64)>, bool_flip_probability: f64) -> Self {
        Self {
            inner,
            f64_clamp_range,
            bool_flip_probability,
            rng: thread_rng(),
        }
    }
}

impl<H: Handler> Handler for FilteringHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let mut value = self.inner.on_sample_f64(addr, dist);

        // Apply clamping if specified
        if let Some((min, max)) = self.f64_clamp_range {
            value = value.clamp(min, max);
        }

        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let mut value = self.inner.on_sample_bool(addr, dist);

        // Flip boolean with specified probability
        if self.rng.gen::<f64>() < self.bool_flip_probability {
            value = !value;
        }

        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        self.inner.on_sample_u64(addr, dist)
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        self.inner.on_sample_usize(addr, dist)
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.inner.on_observe_f64(addr, dist, value);
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.inner.on_observe_bool(addr, dist, value);
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.inner.on_observe_u64(addr, dist, value);
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.inner.on_observe_usize(addr, dist, value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.inner.on_factor(log_weight);
    }

    fn finish(self) -> Trace {
        self.inner.finish()
    }
}

Filtering Use Cases:

• Parameter clamping for numerical stability
• Outlier detection and handling
• Domain-specific constraints enforcement
• Robustness testing through perturbations

Performance Monitoring

Track and optimize computational characteristics with monitoring handlers:

use std::time::{Duration, Instant};

/// Handler that monitors performance characteristics
struct PerformanceHandler<H: Handler> {
    inner: H,
    start_time: Instant,
    operation_times: Vec<Duration>,
    sample_count: u32,
    observe_count: u32,
}

impl<H: Handler> PerformanceHandler<H> {
    fn new(inner: H) -> Self {
        Self {
            inner,
            start_time: Instant::now(),
            operation_times: Vec::new(),
            sample_count: 0,
            observe_count: 0,
        }
    }

    fn time_operation<F, R>(&mut self, operation: F) -> R
    where
        F: FnOnce(&mut H) -> R,
    {
        let start = Instant::now();
        let result = operation(&mut self.inner);
        let duration = start.elapsed();
        self.operation_times.push(duration);
        result
    }
}

impl<H: Handler> Handler for PerformanceHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        self.sample_count += 1;
        self.time_operation(|inner| inner.on_sample_f64(addr, dist))
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        self.sample_count += 1;
        self.time_operation(|inner| inner.on_sample_bool(addr, dist))
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        self.sample_count += 1;
        self.time_operation(|inner| inner.on_sample_u64(addr, dist))
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        self.sample_count += 1;
        self.time_operation(|inner| inner.on_sample_usize(addr, dist))
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.observe_count += 1;
        self.time_operation(|inner| inner.on_observe_f64(addr, dist, value))
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.observe_count += 1;
        self.time_operation(|inner| inner.on_observe_bool(addr, dist, value))
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.observe_count += 1;
        self.time_operation(|inner| inner.on_observe_u64(addr, dist, value))
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.observe_count += 1;
        self.time_operation(|inner| inner.on_observe_usize(addr, dist, value))
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.time_operation(|inner| inner.on_factor(log_weight))
    }

    fn finish(self) -> Trace {
        let total_time = self.start_time.elapsed();
        let avg_op_time = if !self.operation_times.is_empty() {
            self.operation_times.iter().sum::<Duration>() / self.operation_times.len() as u32
        } else {
            Duration::ZERO
        };

        println!("✅ Performance Monitoring Results:");
        println!("   - Total execution time: {:?}", total_time);
        println!("   - Operations performed: {}", self.operation_times.len());
        println!("   - Sample operations: {}", self.sample_count);
        println!("   - Observe operations: {}", self.observe_count);
        println!("   - Average operation time: {:?}", avg_op_time);

        self.inner.finish()
    }
}

Performance Insights:

• Operation timing and bottleneck identification
• Memory allocation patterns
• Execution hotspots and optimization opportunities
• Scalability analysis for production deployment

Custom Inference Algorithms

Custom inference algorithms extend Fugue's effect interpretation to implement novel sampling strategies and approximate inference methods. Each algorithm provides a unique semantic mapping from probabilistic effects to computational actions:

graph TD
    subgraph "Inference Algorithm Architecture"
        A[Model M] --> B[Effect Sequence]
        B --> C{Handler Type}
        C -->|MCMC| D[Markov Chain<br/>Sampling]
        C -->|VI| E[Variational<br/>Approximation]
        C -->|SMC| F[Sequential<br/>Monte Carlo]
        C -->|ABC| G[Approximate<br/>Bayesian Computation]

        D --> H[Posterior Samples]
        E --> I[Approximate<br/>Distribution]
        F --> J[Weighted<br/>Particles]
        G --> K[Likelihood-Free<br/>Samples]
    end

Algorithm Design Principles:

1. Effect Consistency: $\forall e\in \mathcal{E}:h(e)$ preserves probabilistic semantics
2. Convergence Guarantees: Algorithm converges to target distribution under regularity conditions
3. Computational Tractability: Runtime complexity is polynomial in problem dimensions
4. Statistical Efficiency: Effective sample size scales appropriately with computational cost

Mathematical Framework: Each inference handler implements a stochastic operator $T:\mathcal{P}(\Theta )\to \mathcal{P}(\Theta )$ with fixed point $\pi$ such that $T\pi =\pi$.

/// Simple custom MCMC-like handler that perturbs existing values
struct SimpleMCMCHandler<R: Rng> {
    rng: R,
    base_trace: Trace,
    current_trace: Trace,
    perturbation_scale: f64,
}

impl<R: Rng> SimpleMCMCHandler<R> {
    fn new(rng: R, base_trace: Trace, perturbation_scale: f64) -> Self {
        Self {
            rng,
            base_trace,
            current_trace: Trace::default(),
            perturbation_scale,
        }
    }
}

impl<R: Rng> Handler for SimpleMCMCHandler<R> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let value = if let Some(base_value) = self.base_trace.get_f64(addr) {
            // Perturb existing value
            let perturbation = Normal::new(0.0, self.perturbation_scale).unwrap();
            base_value + perturbation.sample(&mut self.rng)
        } else {
            // Sample fresh if not in base trace
            dist.sample(&mut self.rng)
        };

        let log_prob = dist.log_prob(&value);
        self.current_trace.log_prior += log_prob;
        self.current_trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::F64(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        let value = if let Some(base_value) = self.base_trace.get_bool(addr) {
            // Maybe flip the boolean with small probability
            if self.rng.gen::<f64>() < 0.1 {
                !base_value
            } else {
                base_value
            }
        } else {
            dist.sample(&mut self.rng)
        };

        let log_prob = dist.log_prob(&value);
        self.current_trace.log_prior += log_prob;
        self.current_trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Bool(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        // For simplicity, just use base value or sample fresh
        let value = self
            .base_trace
            .get_u64(addr)
            .unwrap_or_else(|| dist.sample(&mut self.rng));

        let log_prob = dist.log_prob(&value);
        self.current_trace.log_prior += log_prob;
        self.current_trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::U64(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        let value = self
            .base_trace
            .get_usize(addr)
            .unwrap_or_else(|| dist.sample(&mut self.rng));

        let log_prob = dist.log_prob(&value);
        self.current_trace.log_prior += log_prob;
        self.current_trace.choices.insert(
            addr.clone(),
            Choice {
                addr: addr.clone(),
                value: ChoiceValue::Usize(value),
                logp: log_prob,
            },
        );

        value
    }

    fn on_observe_f64(&mut self, _addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.current_trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_bool(&mut self, _addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.current_trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_u64(&mut self, _addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.current_trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_observe_usize(&mut self, _addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.current_trace.log_likelihood += dist.log_prob(&value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.current_trace.log_factors += log_weight;
    }

    fn finish(self) -> Trace {
        self.current_trace
    }
}

Inference Handler Patterns:

• MCMC variants: Custom proposal mechanisms and acceptance criteria
• Variational methods: Gradient-based optimization with custom families
• Rejection sampling: Domain-specific acceptance/rejection logic
• Importance sampling: Custom proposal distributions and weight calculations

Handler Composition and Chaining

Handler chaining implements multi-stage effect processing through systematic composition operators. The composition forms a computational pipeline with well-defined data flow and effect propagation:

graph TD
    subgraph "Handler Composition Pipeline"
        A[Raw Effect E] --> B[Statistics Handler]
        B --> C[Logging Handler]
        C --> D[Performance Handler]
        D --> E[Base Handler]
        E --> F[Result + Trace]

        B -.->|Metrics| G[(Statistics DB)]
        C -.->|Events| H[(Log Stream)]
        D -.->|Timing| I[(Performance Monitor)]

        F --> J{Validation}
        J -->|Pass| K[Success]
        J -->|Fail| L[Error Recovery]
    end

Composition Laws:

1. Preservation: $h_n\circ \dots \circ h_1$ preserves effect semantics
2. Associativity: Composition order affects performance but not correctness
3. Commutativity: Decorators with disjoint side effects commute
4. Distributivity: $h\circ (g_1+g_2)=h\circ g_1+h\circ g_2$ for effect unions

Performance Analysis: Handler chain depth $d$ introduces overhead $\mathcal{O}(d\cdot c)$ where $c$ is the per-handler cost. Optimization strategies include handler fusion and effect batching.

fn main() {
    println!("=== Custom Handlers in Fugue ===\n");

    println!("1. Basic Custom Handler Implementation");
    println!("------------------------------------");

    // Test the basic handler
    let mut rng = thread_rng();
    let handler = BasicHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };

    let test_model = || sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
    let (result, trace) = runtime::handler::run(handler, test_model());

    println!("✅ Basic handler executed");
    println!("   - Result: {:.3}", result);
    println!("   - Trace choices: {}", trace.choices.len());
    println!("   - Total log-weight: {:.3}", trace.total_log_weight());
    println!();

    println!("2. Logging Handler - Decorator Pattern");
    println!("-------------------------------------");

    // Test the logging handler
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };
    let logging_handler = LoggingHandler::new(base_handler, false); // Non-verbose

    let logged_model = || {
        prob!(
            let x <- sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());
            observe(addr!("obs"), Normal::new(x, 0.5).unwrap(), 1.2);
            factor(-0.5);
            pure(x)
        )
    };

    let (result, _trace) = runtime::handler::run(logging_handler, logged_model());
    println!("   - Logged execution result: {:.3}", result);
    println!();

    println!("3. Statistics Accumulating Handler");
    println!("--------------------------------");

    // Test the statistics handler
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };
    let stats_handler = StatisticsHandler::new(base_handler);

    let complex_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 2.0).unwrap());
            let is_outlier <- sample(addr!("outlier"), Bernoulli::new(0.1).unwrap());
            let count <- sample(addr!("count"), Poisson::new(3.0).unwrap());
            let category <- sample(addr!("category"), Categorical::new(vec![0.3, 0.4, 0.3]).unwrap());

            observe(addr!("y1"), Normal::new(mu, 1.0).unwrap(), 1.5);
            observe(addr!("y2"), Normal::new(mu, 1.0).unwrap(), 2.1);
            factor(if is_outlier { -2.0 } else { 0.0 });

            pure((mu, is_outlier, count, category))
        )
    };

    let (result, _trace) = runtime::handler::run(stats_handler, complex_model());
    println!(
        "   - Complex model result: {:?}",
        (result.0.round(), result.1, result.2, result.3)
    );
    println!();

    println!("4. Conditional Filtering Handler");
    println!("-------------------------------");

    // Test the filtering handler
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };
    let filtering_handler = FilteringHandler::new(
        base_handler,
        Some((-2.0, 2.0)), // Clamp f64 values to [-2, 2]
        0.1,               // 10% chance to flip booleans
    );

    let filter_test_model = || {
        prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 5.0).unwrap()); // Wide distribution
            let flag <- sample(addr!("flag"), Bernoulli::new(0.8).unwrap());
            pure((x, flag))
        )
    };

    let (result, _trace) = runtime::handler::run(filtering_handler, filter_test_model());
    println!("✅ Filtering handler executed");
    println!("   - Clamped value: {:.3} (should be in [-2, 2])", result.0);
    println!(
        "   - Boolean value: {} (may be flipped from original)",
        result.1
    );
    println!();

    println!("5. Performance Monitoring Handler");
    println!("--------------------------------");

    // Test the performance handler
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };
    let perf_handler = PerformanceHandler::new(base_handler);

    let perf_test_model = || {
        plate!(i in 0..10 => {
            sample(addr!("param", i), Normal::new(0.0, 1.0).unwrap())
        })
    };

    let (_result, _trace) = runtime::handler::run(perf_handler, perf_test_model());
    println!();

    println!("6. Custom Inference Handler");
    println!("---------------------------");

    // Test the custom inference handler
    let mut rng1 = thread_rng();
    let rng2 = thread_rng();

    // First get a base trace
    let base_handler = PriorHandler {
        rng: &mut rng1,
        trace: Trace::default(),
    };

    let inference_model = || {
        prob!(
            let mu <- sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap());
            observe(addr!("y"), Normal::new(mu, 0.5).unwrap(), 1.0);
            pure(mu)
        )
    };

    let (base_result, base_trace) = runtime::handler::run(base_handler, inference_model());

    // Now use custom MCMC handler to perturb it
    let base_log_weight = base_trace.total_log_weight();
    let mcmc_handler = SimpleMCMCHandler::new(rng2, base_trace, 0.1);
    let (mcmc_result, mcmc_trace) = runtime::handler::run(mcmc_handler, inference_model());

    println!("✅ Custom MCMC-like inference:");
    println!("   - Base result: {:.3}", base_result);
    println!("   - MCMC result: {:.3}", mcmc_result);
    println!("   - Base log-weight: {:.3}", base_log_weight);
    println!("   - MCMC log-weight: {:.3}", mcmc_trace.total_log_weight());
    println!();

    println!("7. Handler Composition and Chaining");
    println!("----------------------------------");

    // Demonstrate composing multiple handler decorators
    let mut rng = thread_rng();
    let base_handler = PriorHandler {
        rng: &mut rng,
        trace: Trace::default(),
    };

    // Chain multiple decorators: Statistics -> Logging -> Performance -> Base
    let stats_handler = StatisticsHandler::new(base_handler);
    let logging_handler = LoggingHandler::new(stats_handler, false);
    let performance_handler = PerformanceHandler::new(logging_handler);

    let composition_model = || {
        prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            let y <- sample(addr!("y"), Bernoulli::new(0.7).unwrap());
            observe(addr!("obs"), Normal::new(x, 0.2).unwrap(), 0.5);
            factor(-0.3);
            pure((x, y))
        )
    };

    println!("✅ Handler composition example:");
    let (_result, _trace) = runtime::handler::run(performance_handler, composition_model());
    println!("   - Multiple handler layers executed successfully");
    println!();

    println!("=== Custom Handler Patterns Demonstrated! ===");
}

Composition Strategies:

• Layered approach: Statistics → Logging → Performance → Base
• Conditional activation: Enable decorators based on environment/configuration
• Feature flags: Runtime selection of handler combinations
• Pipeline optimization: Order decorators for minimal overhead

Advanced Handler Patterns

Caching Handler

struct CachingHandler<H: Handler> {
    inner: H,
    cache: HashMap<(Address, String), ChoiceValue>, // Address + dist info -> cached value
}

Distributed Handler

struct DistributedHandler<H: Handler> {
    inner: H,
    worker_id: usize,
    coordinator: Arc<Mutex<SharedState>>,
}

Fault-Tolerant Handler

struct FaultTolerantHandler<H: Handler> {
    inner: H,
    fallback_strategy: FallbackMode,
    error_count: u32,
    max_errors: u32,
}

Testing Custom Handlers

Systematic testing ensures handler correctness:

    #[test]
    fn test_basic_custom_handler() {
        let mut rng = thread_rng();
        let handler = BasicHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };

        let model = sample(addr!("test"), Normal::new(0.0, 1.0).unwrap());
        let (result, trace) = runtime::handler::run(handler, model);

        assert!(trace.choices.contains_key(&addr!("test")));
        assert!(trace.total_log_weight().is_finite());
        assert!(result.is_finite());
    }

    #[test]
    fn test_logging_handler() {
        let mut rng = thread_rng();
        let base_handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };
        let logging_handler = LoggingHandler::new(base_handler, false);

        let model = prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            observe(addr!("obs"), Normal::new(x, 0.1).unwrap(), 1.0);
            pure(x)
        );

        let (result, trace) = runtime::handler::run(logging_handler, model);

        assert!(trace.choices.contains_key(&addr!("x")));
        assert!(trace.log_likelihood.is_finite());
        assert!(result.is_finite());
    }

    #[test]
    fn test_statistics_handler() {
        let mut rng = thread_rng();
        let base_handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };
        let stats_handler = StatisticsHandler::new(base_handler);

        let model = prob!(
            let x <- sample(addr!("x"), Normal::new(0.0, 1.0).unwrap());
            let flag <- sample(addr!("flag"), Bernoulli::new(0.5).unwrap());
            pure((x, flag))
        );

        let (result, trace) = runtime::handler::run(stats_handler, model);

        assert_eq!(trace.choices.len(), 2);
        assert!(result.0.is_finite());
    }

    #[test]
    fn test_handler_composition() {
        let mut rng = thread_rng();
        let base_handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };

        // Compose multiple handlers
        let logged_handler = LoggingHandler::new(base_handler, false);
        let stats_handler = StatisticsHandler::new(logged_handler);

        let model = prob!(
            let x <- sample(addr!("param"), Normal::new(0.0, 1.0).unwrap());
            factor(-0.5);
            pure(x)
        );

        let (result, trace) = runtime::handler::run(stats_handler, model);

        assert!(trace.choices.contains_key(&addr!("param")));
        assert!(trace.log_factors.abs() > 0.0); // Factor was applied
        assert!(result.is_finite());
    }

Testing Strategy:

• Unit tests: Individual handler method behavior
• Integration tests: Handler with realistic models
• Property tests: Invariant verification across random inputs
• Composition tests: Multi-layer handler combinations

Production Considerations

Error Handling

impl Handler for ProductionHandler {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        match self.inner.on_sample_f64(addr, dist) {
            value if value.is_finite() => value,
            _ => {
                self.log_error(addr, "Non-finite sample");
                0.0 // Safe fallback
            }
        }
    }
}

Memory Management

struct MemoryEfficientHandler<H: Handler> {
    inner: H,
    choice_pool: Vec<Choice>, // Reusable allocations
    max_trace_size: usize,
}

Monitoring Integration

struct MetricsHandler<H: Handler> {
    inner: H,
    metrics_client: MetricsClient,
    model_name: String,
}

Common Patterns Summary

1. Decorator Pattern: Wrap handlers for additional functionality
2. State Accumulation: Track statistics and model behavior
3. Conditional Logic: Apply domain-specific constraints
4. Performance Monitoring: Identify bottlenecks and optimization opportunities
5. Custom Inference: Implement specialized algorithms
6. Composition: Chain multiple handlers for comprehensive capabilities
7. Error Handling: Graceful degradation and recovery
8. Resource Management: Efficient memory and computation usage

Best Practices

1. Single Responsibility: Each handler should have one clear purpose
2. Composability: Design handlers to work well in combination
3. Type Safety: Leverage Rust's type system for correctness
4. Performance: Minimize overhead in hot paths
5. Error Handling: Fail gracefully with meaningful diagnostics
6. Testing: Comprehensive unit and integration tests
7. Documentation: Clear API contracts and usage examples

Custom handlers transform Fugue from a probabilistic programming framework into a platform for building specialized inference systems, analytics tools, and production-ready probabilistic applications.

Production Deployment

Production deployment of probabilistic models requires reliability engineering, performance optimization, and operational excellence at scale. This guide establishes a mathematical framework for fault tolerance, service reliability, and system observability using Fugue's production-ready infrastructure patterns.

Error Handling and Graceful Degradation

Graceful degradation implements fault tolerance through systematic error recovery and service continuity. The mathematical foundation relies on Markov reliability models and circuit breaker theory:

stateDiagram-v2
    [*] --> Healthy
    Healthy --> Degraded : Error Rate > τ₁
    Degraded --> Failed : Error Rate > τ₂  
    Failed --> Recovery : Time > T_recovery
    Recovery --> Healthy : Success Rate > σ
    Degraded --> Healthy : Error Rate < τ₀
    
    note right of Healthy
        Error Rate: λ < τ₀
        SLA: 99.9%
        Full Functionality
    end note
    
    note right of Degraded  
        Error Rate: τ₀ < λ < τ₁
        SLA: 95%
        Limited Functionality
    end note
    
    note right of Failed
        Error Rate: λ > τ₂
        Circuit Open
        Fallback Mode
    end note

Circuit Breaker Mathematics: The failure rate follows a Poisson process with rate $\lambda (t)$. The circuit breaker transitions based on:

$$P(\text{trip})=1-e^{-{\int }_0^T\lambda (t)dt}$$

Error Budget Model: For SLA target $\alpha$, the error budget is: $$\text{Budget}(t)=(1-\alpha )\cdot t-{\int }_0^t{1}_{\text{error}}(s)ds$$

/// Production-ready handler that gracefully handles failures
struct RobustProductionHandler<H: Handler> {
    inner: H,
    error_count: u32,
    max_errors: u32,
    _fallback_values: HashMap<String, ChoiceValue>,
    circuit_breaker_open: bool,
}

impl<H: Handler> RobustProductionHandler<H> {
    fn new(inner: H, max_errors: u32) -> Self {
        let mut fallback_values = HashMap::new();
        fallback_values.insert("default_f64".to_string(), ChoiceValue::F64(0.0));
        fallback_values.insert("default_bool".to_string(), ChoiceValue::Bool(false));
        fallback_values.insert("default_u64".to_string(), ChoiceValue::U64(0));
        fallback_values.insert("default_usize".to_string(), ChoiceValue::Usize(0));

        Self {
            inner,
            error_count: 0,
            max_errors,
            _fallback_values: fallback_values,
            circuit_breaker_open: false,
        }
    }

    fn handle_error(&mut self, operation: &str, addr: &Address) -> bool {
        self.error_count += 1;
        eprintln!("PRODUCTION ERROR: {} failed at address {}", operation, addr);

        if self.error_count >= self.max_errors {
            self.circuit_breaker_open = true;
            eprintln!("CIRCUIT BREAKER: Too many errors, switching to fallback mode");
        }

        self.circuit_breaker_open
    }

    fn get_fallback_f64(&self, addr: &Address) -> f64 {
        // In production, this might come from a cache, configuration, or ML model
        match addr.as_str() {
            s if s.contains("temperature") => 20.0,
            s if s.contains("price") => 100.0,
            s if s.contains("probability") => 0.5,
            _ => 0.0,
        }
    }
}

impl<H: Handler> Handler for RobustProductionHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        if self.circuit_breaker_open {
            return self.get_fallback_f64(addr);
        }

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            self.inner.on_sample_f64(addr, dist)
        })) {
            Ok(value) if value.is_finite() => value,
            Ok(invalid_value) => {
                eprintln!("Invalid f64 sample: {} at {}", invalid_value, addr);
                self.handle_error("sample_f64", addr);
                self.get_fallback_f64(addr)
            }
            Err(_) => {
                self.handle_error("sample_f64_panic", addr);
                self.get_fallback_f64(addr)
            }
        }
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        if self.circuit_breaker_open {
            return false; // Safe fallback
        }

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            self.inner.on_sample_bool(addr, dist)
        })) {
            Ok(value) => value,
            Err(_) => {
                self.handle_error("sample_bool_panic", addr);
                false
            }
        }
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        if self.circuit_breaker_open {
            return 1; // Safe default
        }

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            self.inner.on_sample_u64(addr, dist)
        })) {
            Ok(value) => value,
            Err(_) => {
                self.handle_error("sample_u64_panic", addr);
                1
            }
        }
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        if self.circuit_breaker_open {
            return 0; // Safe array index
        }

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            self.inner.on_sample_usize(addr, dist)
        })) {
            Ok(value) => value,
            Err(_) => {
                self.handle_error("sample_usize_panic", addr);
                0
            }
        }
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        if !self.circuit_breaker_open
            && std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                self.inner.on_observe_f64(addr, dist, value)
            }))
            .is_err()
        {
            self.handle_error("observe_f64_panic", addr);
        }
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        if !self.circuit_breaker_open
            && std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                self.inner.on_observe_bool(addr, dist, value)
            }))
            .is_err()
        {
            self.handle_error("observe_bool_panic", addr);
        }
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        if !self.circuit_breaker_open
            && std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                self.inner.on_observe_u64(addr, dist, value)
            }))
            .is_err()
        {
            self.handle_error("observe_u64_panic", addr);
        }
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        if !self.circuit_breaker_open
            && std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                self.inner.on_observe_usize(addr, dist, value)
            }))
            .is_err()
        {
            self.handle_error("observe_usize_panic", addr);
        }
    }

    fn on_factor(&mut self, log_weight: f64) {
        if !log_weight.is_finite() {
            eprintln!("Invalid factor log-weight: {}", log_weight);
            self.error_count += 1;
            return; // Skip invalid factors
        }

        if !self.circuit_breaker_open {
            self.inner.on_factor(log_weight);
        }
    }

    fn finish(self) -> Trace {
        println!("✅ Production handler statistics:");
        println!(
            "   - Errors encountered: {}/{}",
            self.error_count, self.max_errors
        );
        println!(
            "   - Circuit breaker status: {}",
            if self.circuit_breaker_open {
                "OPEN (fallback mode)"
            } else {
                "CLOSED (normal)"
            }
        );

        if self.circuit_breaker_open {
            // Return minimal valid trace in fallback mode
            Trace::default()
        } else {
            self.inner.finish()
        }
    }
}

Robust Error Handling Features:

• Circuit Breaker Pattern: Prevents cascade failures by switching to fallback mode
• Panic Recovery: Catches panics and provides safe default values
• Input Validation: Ensures all inputs are finite and within expected ranges
• Fallback Values: Domain-specific defaults for different parameter types
• Error Counting: Tracks error rates to trigger circuit breaker activation

Configuration Management

Production models require flexible configuration for different environments:

#[derive(Debug, Clone)]
struct ModelConfig {
    // Model parameters
    temperature_prior_mean: f64,
    temperature_prior_std: f64,
    validity_probability: f64,
    sensor_noise_std: f64,

    // Runtime configuration
    max_inference_time_ms: u64,
    memory_pool_size: usize,
    enable_circuit_breaker: bool,
    error_threshold: u32,

    // Environment settings
    environment: String, // "development", "staging", "production"
    _log_level: String,
    enable_metrics: bool,
}

impl Default for ModelConfig {
    fn default() -> Self {
        Self {
            temperature_prior_mean: 20.0,
            temperature_prior_std: 5.0,
            validity_probability: 0.95,
            sensor_noise_std: 1.0,
            max_inference_time_ms: 1000,
            memory_pool_size: 100,
            enable_circuit_breaker: true,
            error_threshold: 10,
            environment: "production".to_string(),
            _log_level: "info".to_string(),
            enable_metrics: true,
        }
    }
}

struct ConfigurableModelRunner {
    config: ModelConfig,
    metrics: ProductionMetrics,
}

impl ConfigurableModelRunner {
    fn new(config: ModelConfig) -> Self {
        Self {
            metrics: ProductionMetrics::new(config.enable_metrics),
            config,
        }
    }

    fn create_model(&self) -> Model<(f64, bool)> {
        let config = self.config.clone();
        prob!(
            let temp <- sample(
                addr!("temperature"),
                Normal::new(config.temperature_prior_mean, config.temperature_prior_std).unwrap()
            );
            let valid <- sample(
                addr!("valid"),
                Bernoulli::new(config.validity_probability).unwrap()
            );
            // Simulate sensor reading with configured noise
            observe(
                addr!("sensor"),
                Normal::new(temp, config.sensor_noise_std).unwrap(),
                22.0
            );
            pure((temp, valid))
        )
    }

    fn run_inference(&mut self) -> Result<(f64, bool), String> {
        let start = Instant::now();

        // Configure handler based on environment
        let mut rng = thread_rng();
        let model = self.create_model(); // Create model before borrowing
        let result = if self.config.environment == "production" {
            // Use safe, fault-tolerant execution in production
            let base_handler = PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            };
            let robust_handler =
                RobustProductionHandler::new(base_handler, self.config.error_threshold);

            let (result, _trace) = runtime::handler::run(robust_handler, model);
            Ok(result)
        } else {
            // Use faster, less safe execution in development
            let handler = PriorHandler {
                rng: &mut rng,
                trace: Trace::default(),
            };
            let (result, _trace) = runtime::handler::run(handler, model);
            Ok(result)
        };

        let duration = start.elapsed();
        self.metrics.record_inference_time(duration);

        // Check timeout
        if duration.as_millis() > self.config.max_inference_time_ms as u128 {
            self.metrics.increment_timeout_count();
            return Err(format!(
                "Inference timeout: {}ms > {}ms",
                duration.as_millis(),
                self.config.max_inference_time_ms
            ));
        }

        result
    }
}

Configuration Best Practices:

• Environment-Specific Settings: Different behavior for development/staging/production
• Model Parameter Configuration: Tunable priors, noise levels, and thresholds
• Runtime Configuration: Timeout limits, error thresholds, retry budgets
• Deployment Configuration: Circuit breaker settings, logging levels, metrics enablement
• Type-Safe Defaults: Sensible fallbacks for all configuration parameters

Production Metrics and Observability

Observability requires systematic metric collection with statistical analysis and anomaly detection. The metric taxonomy follows the USE method (Utilization, Saturation, Errors) and RED method (Rate, Errors, Duration):

graph TD
    subgraph "Observability Architecture"
        A[Model Execution] --> B[Metric Collection]
        B --> C{Metric Type}
        C -->|USE| D[Resource Metrics]
        C -->|RED| E[Service Metrics]  
        C -->|Business| F[Domain Metrics]
        
        D --> G[Utilization: ρ = λ/μ]
        D --> H[Saturation: Queue Length]
        D --> I[Error Rate: λₑ]
        
        E --> J[Request Rate: λᵣ]
        E --> K[Error Rate: εᵣ] 
        E --> L[Duration: T₉₉]
        
        F --> M[Inference Accuracy]
        F --> N[Model Drift]
        F --> O[Business KPIs]
        
        G --> P[(Time Series DB)]
        H --> P
        I --> P
        J --> P
        K --> P
        L --> P
        M --> Q[(Analytics DB)]
        N --> Q
        O --> Q
    end

Statistical Process Control: Metrics follow control chart theory with statistical control limits:

$$\begin{array}{rl}\text{UCL}& =\bar{X}+3\sigma /\sqrt{n}\\ \text{LCL}& =\bar{X}-3\sigma /\sqrt{n}\end{array}$$

Anomaly Detection: Using exponentially weighted moving averages: $${\text{EWMA}}_t=\alpha X_t+(1-\alpha ){\text{EWMA}}_{t-1}$$

#[derive(Debug, Clone)]
struct ProductionMetrics {
    enabled: bool,
    inference_count: u64,
    error_count: u64,
    timeout_count: u64,
    total_inference_time: Duration,
    start_time: SystemTime,
}

impl ProductionMetrics {
    fn new(enabled: bool) -> Self {
        Self {
            enabled,
            inference_count: 0,
            error_count: 0,
            timeout_count: 0,
            total_inference_time: Duration::ZERO,
            start_time: SystemTime::now(),
        }
    }

    fn record_inference_time(&mut self, duration: Duration) {
        if self.enabled {
            self.inference_count += 1;
            self.total_inference_time += duration;
        }
    }

    fn _increment_error_count(&mut self) {
        if self.enabled {
            self.error_count += 1;
        }
    }

    fn increment_timeout_count(&mut self) {
        if self.enabled {
            self.timeout_count += 1;
        }
    }

    fn get_stats(&self) -> HashMap<String, f64> {
        let mut stats = HashMap::new();
        if self.enabled {
            let uptime = self.start_time.elapsed().unwrap_or(Duration::ZERO);
            let avg_inference_time = if self.inference_count > 0 {
                self.total_inference_time.as_millis() as f64 / self.inference_count as f64
            } else {
                0.0
            };

            stats.insert("inference_count".to_string(), self.inference_count as f64);
            stats.insert("error_count".to_string(), self.error_count as f64);
            stats.insert("timeout_count".to_string(), self.timeout_count as f64);
            stats.insert(
                "error_rate".to_string(),
                if self.inference_count > 0 {
                    self.error_count as f64 / self.inference_count as f64
                } else {
                    0.0
                },
            );
            stats.insert("avg_inference_time_ms".to_string(), avg_inference_time);
            stats.insert("uptime_seconds".to_string(), uptime.as_secs() as f64);
            stats.insert(
                "throughput_per_second".to_string(),
                if uptime.as_secs() > 0 {
                    self.inference_count as f64 / uptime.as_secs() as f64
                } else {
                    0.0
                },
            );
        }
        stats
    }

    fn export_prometheus_metrics(&self) -> String {
        let mut metrics = String::new();
        let stats = self.get_stats();

        metrics.push_str("# HELP fugue_inference_total Total number of inference runs\n");
        metrics.push_str("# TYPE fugue_inference_total counter\n");
        metrics.push_str(&format!(
            "fugue_inference_total {}\n",
            stats.get("inference_count").unwrap_or(&0.0)
        ));

        metrics.push_str("# HELP fugue_errors_total Total number of errors\n");
        metrics.push_str("# TYPE fugue_errors_total counter\n");
        metrics.push_str(&format!(
            "fugue_errors_total {}\n",
            stats.get("error_count").unwrap_or(&0.0)
        ));

        metrics.push_str(
            "# HELP fugue_inference_duration_ms Average inference duration in milliseconds\n",
        );
        metrics.push_str("# TYPE fugue_inference_duration_ms gauge\n");
        metrics.push_str(&format!(
            "fugue_inference_duration_ms {}\n",
            stats.get("avg_inference_time_ms").unwrap_or(&0.0)
        ));

        metrics.push_str("# HELP fugue_error_rate Error rate (errors/total inferences)\n");
        metrics.push_str("# TYPE fugue_error_rate gauge\n");
        metrics.push_str(&format!(
            "fugue_error_rate {}\n",
            stats.get("error_rate").unwrap_or(&0.0)
        ));

        metrics
    }
}

/// Production monitoring handler that integrates with metrics systems
struct MetricsHandler<H: Handler> {
    inner: H,
    metrics: Arc<std::sync::Mutex<ProductionMetrics>>,
    _model_name: String,
}

impl<H: Handler> MetricsHandler<H> {
    fn new(
        inner: H,
        metrics: Arc<std::sync::Mutex<ProductionMetrics>>,
        model_name: String,
    ) -> Self {
        Self {
            inner,
            metrics,
            _model_name: model_name,
        }
    }
}

impl<H: Handler> Handler for MetricsHandler<H> {
    fn on_sample_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>) -> f64 {
        let start = Instant::now();
        let result = self.inner.on_sample_f64(addr, dist);

        if let Ok(mut metrics) = self.metrics.lock() {
            metrics.record_inference_time(start.elapsed());
        }

        result
    }

    fn on_sample_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>) -> bool {
        self.inner.on_sample_bool(addr, dist)
    }

    fn on_sample_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>) -> u64 {
        self.inner.on_sample_u64(addr, dist)
    }

    fn on_sample_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>) -> usize {
        self.inner.on_sample_usize(addr, dist)
    }

    fn on_observe_f64(&mut self, addr: &Address, dist: &dyn Distribution<f64>, value: f64) {
        self.inner.on_observe_f64(addr, dist, value);
    }

    fn on_observe_bool(&mut self, addr: &Address, dist: &dyn Distribution<bool>, value: bool) {
        self.inner.on_observe_bool(addr, dist, value);
    }

    fn on_observe_u64(&mut self, addr: &Address, dist: &dyn Distribution<u64>, value: u64) {
        self.inner.on_observe_u64(addr, dist, value);
    }

    fn on_observe_usize(&mut self, addr: &Address, dist: &dyn Distribution<usize>, value: usize) {
        self.inner.on_observe_usize(addr, dist, value);
    }

    fn on_factor(&mut self, log_weight: f64) {
        self.inner.on_factor(log_weight);
    }

    fn finish(self) -> Trace {
        self.inner.finish()
    }
}

Metrics Collection:

• Performance Metrics: Inference time, throughput, operation counts
• Error Tracking: Error rates, timeout counts, failure categorization
• System Health: Uptime, resource utilization, allocation rate
• Prometheus Integration: Standard metrics format for monitoring systems
• Real-Time Dashboards: Live performance and health indicators

Health Checks and System Validation

Health monitoring implements continuous system validation through multi-level health checks with statistical thresholds and predictive alerting:

graph TD
    subgraph "Health Check Hierarchy"
        A[System Health H⁽ˢ⁾] --> B[Model Health H⁽ᵐ⁾]
        B --> C[Inference Health H⁽ⁱ⁾]
        C --> D[Resource Health H⁽ʳ⁾]

        A --> E{H⁽ˢ⁾ > θₛ?}
        B --> F{H⁽ᵐ⁾ > θₘ?}
        C --> G{H⁽ⁱ⁾ > θᵢ?}
        D --> H{H⁽ʳ⁾ > θʳ?}

        E -->|No| I[System Alert]
        F -->|No| J[Model Alert]
        G -->|No| K[Inference Alert]
        H -->|No| L[Resource Alert]

        E -->|Yes| M[System OK]
        F -->|Yes| M
        G -->|Yes| M
        H -->|Yes| M
    end

Health Score Calculation: Weighted combination of subsystem health: $$H^{(s)}=\sum _i{w}_i{H}^{(i)}$$

where $w_i$ are importance weights and $\sum w_i=1$.

Predictive Health Modeling: Using time series forecasting: $$H_{t+k}=\alpha H_t+\beta \frac{dH}{dt}{|}_t+\gamma \frac{d^{2}H}{d{t}^2}{|}_t$$

#[derive(Debug, Clone)]
struct HealthCheckResult {
    status: HealthStatus,
    message: String,
    details: HashMap<String, String>,
    timestamp: SystemTime,
}

#[derive(Debug, Clone, PartialEq)]
enum HealthStatus {
    Healthy,
    Degraded,
    Unhealthy,
}

struct ProductionHealthChecker {
    model_config: ModelConfig,
    metrics: Arc<std::sync::Mutex<ProductionMetrics>>,
}

impl ProductionHealthChecker {
    fn new(config: ModelConfig, metrics: Arc<std::sync::Mutex<ProductionMetrics>>) -> Self {
        Self {
            model_config: config,
            metrics,
        }
    }

    fn run_health_check(&self) -> HealthCheckResult {
        let mut details = HashMap::new();
        let mut overall_status = HealthStatus::Healthy;
        let mut messages = Vec::new();

        // Check 1: Model execution health
        match self.check_model_execution() {
            Ok(duration) => {
                details.insert("model_execution".to_string(), "healthy".to_string());
                details.insert(
                    "execution_time_ms".to_string(),
                    format!("{:.1}", duration.as_millis()),
                );
            }
            Err(e) => {
                overall_status = HealthStatus::Unhealthy;
                messages.push(format!("Model execution failed: {}", e));
                details.insert("model_execution".to_string(), "failed".to_string());
            }
        }

        // Check 2: Error rates
        if let Ok(metrics) = self.metrics.lock() {
            let stats = metrics.get_stats();
            let error_rate = stats.get("error_rate").unwrap_or(&0.0) * 100.0;
            details.insert(
                "error_rate_percent".to_string(),
                format!("{:.2}%", error_rate),
            );

            if error_rate > 5.0 {
                overall_status = HealthStatus::Degraded;
                messages.push(format!("High error rate: {:.1}%", error_rate));
            } else if error_rate > 20.0 {
                overall_status = HealthStatus::Unhealthy;
                messages.push(format!("Critical error rate: {:.1}%", error_rate));
            }

            let avg_time = stats.get("avg_inference_time_ms").unwrap_or(&0.0);
            details.insert(
                "avg_inference_time_ms".to_string(),
                format!("{:.1}", avg_time),
            );

            if *avg_time > self.model_config.max_inference_time_ms as f64 * 0.8 {
                overall_status = HealthStatus::Degraded;
                messages.push("Inference time approaching timeout threshold".to_string());
            }
        }

        // Check 4: System resources
        details.insert(
            "memory_pool_size".to_string(),
            self.model_config.memory_pool_size.to_string(),
        );
        details.insert(
            "circuit_breaker".to_string(),
            if self.model_config.enable_circuit_breaker {
                "enabled".to_string()
            } else {
                "disabled".to_string()
            },
        );

        let message = if messages.is_empty() {
            "All systems healthy".to_string()
        } else {
            messages.join("; ")
        };

        HealthCheckResult {
            status: overall_status,
            message,
            details,
            timestamp: SystemTime::now(),
        }
    }

    fn check_model_execution(&self) -> Result<Duration, String> {
        let start = Instant::now();
        let mut rng = thread_rng();

        // Run a simplified version of the model for health checking
        let health_model = || {
            prob!(
                let value <- sample(addr!("health_check"), Normal::new(0.0, 1.0).unwrap());
                pure(value)
            )
        };

        let handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            runtime::handler::run(handler, health_model())
        })) {
            Ok((result, trace)) => {
                if result.is_finite() && trace.total_log_weight().is_finite() {
                    Ok(start.elapsed())
                } else {
                    Err("Invalid model output".to_string())
                }
            }
            Err(_) => Err("Model execution panicked".to_string()),
        }
    }
}

Health Check Components:

• Model Execution Health: Verifies core functionality with simplified tests
• Memory Health: Monitors allocation rate and memory usage patterns
• Error Rate Analysis: Tracks and categorizes different failure modes
• Performance Monitoring: Identifies degradation before it impacts users
• Multi-Level Status: Healthy/Degraded/Unhealthy with detailed diagnostics

Input Validation and Security

Robust input validation prevents security vulnerabilities and system failures:

/// Secure input validator for production model parameters
struct InputValidator;

impl InputValidator {
    fn validate_temperature(temp: f64) -> Result<f64, String> {
        match temp {
            t if !t.is_finite() => Err("Temperature must be finite".to_string()),
            t if t < -50.0 => Err("Temperature too low (< -50°C)".to_string()),
            t if t > 100.0 => Err("Temperature too high (> 100°C)".to_string()),
            t => Ok(t),
        }
    }

    fn validate_probability(p: f64) -> Result<f64, String> {
        match p {
            p if !p.is_finite() => Err("Probability must be finite".to_string()),
            p if p < 0.0 => Err("Probability must be non-negative".to_string()),
            p if p > 1.0 => Err("Probability must not exceed 1.0".to_string()),
            p => Ok(p),
        }
    }

    fn _validate_sensor_reading(reading: f64) -> Result<f64, String> {
        match reading {
            r if !r.is_finite() => Err("Sensor reading must be finite".to_string()),
            r if r.abs() > 1000.0 => Err("Sensor reading out of reasonable range".to_string()),
            r => Ok(r),
        }
    }

    fn sanitize_address_component(component: &str) -> Result<String, String> {
        // Prevent injection attacks in address components
        if component
            .chars()
            .any(|c| !(c.is_alphanumeric() || c == '_' || c == '-'))
        {
            return Err("Address component contains invalid characters".to_string());
        }

        if component.len() > 50 {
            Err("Address component too long".to_string())
        } else if component.is_empty() {
            Err("Address component cannot be empty".to_string())
        } else {
            Ok(component.to_string())
        }
    }
}

/// Production model with comprehensive input validation
fn create_validated_model(
    temperature_reading: f64,
    sensor_id: &str,
    prior_prob: f64,
) -> Result<Model<(f64, bool)>, String> {
    // Validate all inputs before model creation
    let validated_temp = InputValidator::validate_temperature(temperature_reading)?;
    let validated_prob = InputValidator::validate_probability(prior_prob)?;
    let sanitized_sensor_id = InputValidator::sanitize_address_component(sensor_id)?;

    // Additional business logic validation
    if sanitized_sensor_id.starts_with("test_") && validated_prob > 0.5 {
        return Err("Test sensors cannot have high prior probability".to_string());
    }

    Ok(prob!(
        let true_temp <- sample(
            addr!("temperature"),
            Normal::new(validated_temp, 2.0).unwrap()
        );
        let is_working <- sample(
            addr!("sensor_working", sanitized_sensor_id.clone()),
            Bernoulli::new(validated_prob).unwrap()
        );

        // Safe observation with validated input
        observe(
            addr!("reading", sanitized_sensor_id),
            Normal::new(true_temp, if is_working { 0.5 } else { 5.0 }).unwrap(),
            validated_temp
        );

        pure((true_temp, is_working))
    ))
}

Security Measures:

• Range Validation: Ensure parameters are within physically meaningful bounds
• Type Safety: Validate all inputs before model construction
• Sanitization: Clean address components to prevent injection attacks
• Business Rule Enforcement: Domain-specific validation logic
• Error Messages: Informative feedback without revealing system internals

Deployment Strategies and Patterns

Deployment strategies implement risk management through controlled rollout and statistical validation. Each strategy provides different risk-latency tradeoffs:

graph TD
    subgraph "Deployment Strategy Matrix"  
        A[New Model Version] --> B{Strategy Selection}
        B -->|Low Risk| C[Blue-Green]
        B -->|Medium Risk| D[Canary]
        B -->|High Risk| E[A/B Test]

        C --> F[Instant Switch<br/>Risk: High<br/>Rollback: Fast]
        D --> G[Gradual Rollout<br/>Risk: Medium<br/>Validation: Statistical]
        E --> H[Statistical Test<br/>Risk: Low<br/>Duration: Long]

        F --> I{Success?}
        G --> J{Performance > Baseline?}
        E --> K{Significance Test?}

        I -->|No| L[Instant Rollback]
        J -->|No| M[Gradual Rollback]
        K -->|No| N[Maintain Status Quo]

        I -->|Yes| O[Full Deployment]
        J -->|Yes| P[Continue Rollout]
        K -->|Yes| Q[Gradual Migration]
    end

Canary Analysis: Statistical significance testing for canary deployments:

$$\text{Z-score}=\frac{({\bar{X}}_{\text{canary}}-{\bar{X}}_{\text{control}})}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}}$$

A/B Testing: Welch's t-test for unequal variances: $$t=\frac{{\bar{X}}_A-{\bar{X}}_B}{\sqrt{\frac{s_A^2}{n_A}+\frac{s_B^2}{n_B}}}$$

/// Production deployment manager with different strategies
#[derive(Debug, Clone)]
enum DeploymentStrategy {
    BlueGreen,
    CanaryRelease { percentage: f64 },
    RollingUpdate,
    _ImmediateSwitch,
}

struct ModelDeploymentManager {
    current_model_version: String,
    candidate_model_version: String,
    deployment_strategy: DeploymentStrategy,
    _rollback_threshold_error_rate: f64,
}

impl ModelDeploymentManager {
    fn new(strategy: DeploymentStrategy) -> Self {
        Self {
            current_model_version: "v1.0.0".to_string(),
            candidate_model_version: "v1.1.0".to_string(),
            deployment_strategy: strategy,
            _rollback_threshold_error_rate: 0.05, // 5% error rate triggers rollback
        }
    }

    fn should_use_candidate_model(&self, request_id: u64) -> bool {
        match &self.deployment_strategy {
            DeploymentStrategy::BlueGreen => {
                // In blue-green, we typically switch all traffic at once
                // For demo, we'll use request ID to simulate the switch
                request_id % 100 < 10 // 10% to candidate for testing
            }
            DeploymentStrategy::CanaryRelease { percentage } => {
                let hash = request_id % 100;
                (hash as f64) < (*percentage * 100.0)
            }
            DeploymentStrategy::RollingUpdate => {
                // Gradual rollout based on some criteria
                request_id % 10 < 3 // 30% rollout
            }
            DeploymentStrategy::_ImmediateSwitch => true,
        }
    }

    fn create_model(&self, use_candidate: bool) -> impl Fn() -> Model<f64> {
        let version = if use_candidate {
            self.candidate_model_version.clone()
        } else {
            self.current_model_version.clone()
        };

        move || {
            if version.starts_with("v1.1") {
                // Candidate model with improved parameters
                prob!(
                    let value <- sample(addr!("improved_param"), Normal::new(0.0, 0.8).unwrap());
                    factor(0.1); // Slight preference for this model
                    pure(value)
                )
            } else {
                // Current stable model
                prob!(
                    let value <- sample(addr!("stable_param"), Normal::new(0.0, 1.0).unwrap());
                    pure(value)
                )
            }
        }
    }

    fn process_request(&self, request_id: u64) -> Result<(f64, String), String> {
        let use_candidate = self.should_use_candidate_model(request_id);
        let version = if use_candidate {
            &self.candidate_model_version
        } else {
            &self.current_model_version
        };

        let model = self.create_model(use_candidate);

        // Execute with error handling
        let mut rng = thread_rng();
        let handler = PriorHandler {
            rng: &mut rng,
            trace: Trace::default(),
        };

        match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            runtime::handler::run(handler, model())
        })) {
            Ok((result, trace)) => {
                if result.is_finite() && trace.total_log_weight().is_finite() {
                    Ok((result, version.clone()))
                } else {
                    Err(format!("Invalid result from model {}", version))
                }
            }
            Err(_) => Err(format!("Model {} panicked", version)),
        }
    }
}

Deployment Patterns:

• Blue-Green Deployment: Instant traffic switching between model versions
• Canary Releases: Gradual rollout to percentage of traffic for risk mitigation
• Rolling Updates: Progressive deployment across infrastructure
• A/B Testing: Compare model performance with statistical significance
• Rollback Capability: Quick reversion to previous version on issues

Performance Optimization Patterns

Inference Execution

Each inference run uses a handler over a fresh Trace. PriorHandler is the standard forward-execution handler:

use fugue::runtime::interpreters::PriorHandler;

// Standard per-run execution
let handler = PriorHandler { rng: &mut rng, trace: Trace::default() };
let (result, trace) = runtime::handler::run(handler, model);

Batch Processing

// Process multiple inference requests efficiently
struct BatchProcessor {
    batch_size: usize,
}

impl BatchProcessor {
    fn process_batch(&mut self, requests: Vec<InferenceRequest>) -> Vec<InferenceResult> {
        requests.into_iter().map(|mut req| {
            let handler = PriorHandler { rng: &mut req.rng, trace: Trace::default() };
            self.run_single_inference(handler, req.model)
        }).collect()
    }
}

Connection Pooling

// Database connection management for model parameters
struct ModelParameterStore {
    connection_pool: Arc<ConnectionPool>,
    parameter_cache: LruCache<String, ModelParameters>,
}

Monitoring Integration

Prometheus Metrics

// Export metrics in Prometheus format
fn export_metrics(metrics: &ProductionMetrics) -> String {
    format!(
        "# HELP fugue_inference_total Total inference operations\n\
         TYPE fugue_inference_total counter\n\
         fugue_inference_total {}\n\
         HELP fugue_error_rate Current error rate\n\
         TYPE fugue_error_rate gauge\n\
         fugue_error_rate {}\n",
        metrics.inference_count,
        metrics.error_rate()
    )
}

Structured Logging

use serde_json::json;

// Structured logging for production debugging
fn log_inference_event(
    request_id: &str,
    model_version: &str,
    duration: Duration,
    result: &InferenceResult
) {
    let log_entry = json!({
        "event": "inference_completed",
        "request_id": request_id,
        "model_version": model_version,
        "duration_ms": duration.as_millis(),
        "success": result.is_success(),
        "timestamp": SystemTime::now(),
    });
    println!("{}", log_entry);
}

Alert Rules

// Define alerting thresholds
struct AlertRules {
    max_error_rate: f64,
    max_latency_ms: u64,
    min_throughput_per_sec: f64,
}

impl AlertRules {
    fn check_alerts(&self, metrics: &ProductionMetrics) -> Vec<Alert> {
        let mut alerts = Vec::new();

        if metrics.error_rate() > self.max_error_rate {
            alerts.push(Alert::HighErrorRate(metrics.error_rate()));
        }

        if metrics.avg_latency().as_millis() > self.max_latency_ms as u128 {
            alerts.push(Alert::HighLatency(metrics.avg_latency()));
        }

        alerts
    }
}

Testing in Production

Shadow Mode Testing

// Run new model versions in shadow mode
struct ShadowTester {
    primary_model: Box<dyn Fn() -> Model<f64>>,
    shadow_model: Box<dyn Fn() -> Model<f64>>,
    comparison_rate: f64,
}

impl ShadowTester {
    fn run_with_shadow(&mut self, input: &Input) -> (PrimaryResult, Option<ShadowResult>) {
        let primary = self.run_primary(input);

        let shadow = if rand::random::<f64>() < self.comparison_rate {
            Some(self.run_shadow(input))
        } else {
            None
        };

        (primary, shadow)
    }
}

Production Validation

// Continuous validation in production
fn validate_model_assumptions(trace: &Trace) -> ValidationResult {
    let mut issues = Vec::new();

    // Check log-weight stability
    if !trace.total_log_weight().is_finite() {
        issues.push("Non-finite log-weight detected".to_string());
    }

    // Check parameter ranges
    for (addr, choice) in &trace.choices {
        if let ChoiceValue::F64(value) = choice.value {
            if value.abs() > 1000.0 {
                issues.push(format!("Extreme value at {}: {}", addr, value));
            }
        }
    }

    ValidationResult { issues }
}

Operational Excellence

Infrastructure as Code

# Kubernetes deployment example
apiVersion: apps/v1
kind: Deployment
metadata:
  name: fugue-inference-service
spec:
  replicas: 3
  selector:
    matchLabels:
      app: fugue-inference
  template:
    metadata:
      labels:
        app: fugue-inference
    spec:
      containers:
      - name: inference-service
        image: fugue-inference:v1.2.0
        resources:
          requests:
            memory: "256Mi"
            cpu: "250m"
          limits:
            memory: "512Mi"
            cpu: "500m"
        livenessProbe:
          httpGet:
            path: /health
            port: 8080
          initialDelaySeconds: 30
          periodSeconds: 10

Service Level Objectives (SLOs)

// Define and monitor SLOs
struct ServiceLevelObjectives {
    availability_target: f64,    // 99.9%
    latency_p99_ms: u64,        // 100ms
    error_rate_threshold: f64,   // 0.1%
}

impl ServiceLevelObjectives {
    fn evaluate_slo_compliance(&self, metrics: &ProductionMetrics) -> SLOReport {
        SLOReport {
            availability: self.calculate_availability(metrics),
            latency_compliance: metrics.p99_latency() <= Duration::from_millis(self.latency_p99_ms),
            error_rate_compliance: metrics.error_rate() <= self.error_rate_threshold,
        }
    }
}

Capacity Planning

// Capacity planning and auto-scaling
struct CapacityPlanner {
    target_cpu_utilization: f64,
    target_memory_utilization: f64,
    scale_up_threshold: f64,
    scale_down_threshold: f64,
}

impl CapacityPlanner {
    fn recommend_scaling(&self, current_metrics: &SystemMetrics) -> ScalingRecommendation {
        if current_metrics.cpu_utilization > self.scale_up_threshold {
            ScalingRecommendation::ScaleUp(self.calculate_scale_factor(current_metrics))
        } else if current_metrics.cpu_utilization < self.scale_down_threshold {
            ScalingRecommendation::ScaleDown(0.5)
        } else {
            ScalingRecommendation::NoAction
        }
    }
}

Security Best Practices

Input Sanitization

Always validate and sanitize inputs before processing:

• Range checks for numerical parameters
• Character filtering for string inputs
• Business rule validation for domain constraints
• Rate limiting to prevent abuse
• Authentication and authorization for API access

Secret Management

// Secure configuration management
struct SecureConfig {
    database_url: SecretString,
    api_key: SecretString,
    model_parameters: ModelConfig,
}

impl SecureConfig {
    fn from_environment() -> Result<Self, ConfigError> {
        Ok(SecureConfig {
            database_url: env::var("DATABASE_URL")?.into(),
            api_key: env::var("API_KEY")?.into(),
            model_parameters: ModelConfig::from_file("model_config.toml")?,
        })
    }
}

Audit Logging

// Comprehensive audit trail
fn log_inference_request(
    user_id: &str,
    request: &InferenceRequest,
    response: &InferenceResponse,
) {
    let audit_log = AuditLogEntry {
        timestamp: SystemTime::now(),
        user_id: user_id.to_string(),
        action: "inference_request".to_string(),
        input_hash: hash_sensitive_data(&request.input),
        output_hash: hash_sensitive_data(&response.output),
        model_version: response.model_version.clone(),
        success: response.success,
    };

    audit_logger::log(audit_log);
}

Common Production Pitfalls

Redundant Per-Request Work

// Avoid: rebuilding immutable configuration on every request
// for request in requests {
//     let config = load_model_config(); // re-parsed every time!
//     process_request(config, request);
// }

// Do: build shared, immutable state once and reuse it. Each run still gets its
// own fresh Trace via the handler.
let config = load_model_config();
for mut request in requests {
    let handler = PriorHandler { rng: &mut request.rng, trace: Trace::default() };
    process_request(&config, handler, request.model);
}

Blocking Operations

// Avoid: Synchronous database calls in request handlers
// let result = database.query_sync(query); // Blocks event loop

// Do: Use async operations with proper timeouts
async fn process_request(request: Request) -> Result<Response, Error> {
    let timeout = Duration::from_millis(100);
    let result = tokio::time::timeout(timeout, database.query(query)).await??;
    Ok(result)
}

Error Propagation

// Avoid: Panicking on errors
// let value = risky_operation().unwrap(); // May crash service

// Do: Graceful error handling with fallbacks
let value = match risky_operation() {
    Ok(v) => v,
    Err(e) => {
        metrics.increment_error_count();
        log::warn!("Operation failed: {}, using fallback", e);
        fallback_value()
    }
};

Production deployment represents the culmination of probabilistic programming excellence, where theoretical foundations meet operational reality. Fugue's comprehensive tooling transforms academic probabilistic models into production-grade systems that deliver reliable, scalable, and maintainable probabilistic computing solutions.

Tutorials

• Foundation Tutorials - Bayesian inference, type safety, and traces
• Statistical Modeling - Regression, classification, mixtures, hierarchies
• Advanced Inference - SMC, ABC, and Variational Inference on the same worked example

Foundation Tutorials

Welcome to Fugue's Foundation Tutorials — your comprehensive introduction to the core concepts and unique features that make Fugue a revolutionary approach to probabilistic programming.

What You'll Learn

These tutorials build upon each other to give you a complete understanding of Fugue's foundational principles:

graph TB
    A["🪙 Bayesian Coin Flip<br/>Basic Probabilistic Modeling"] --> B["🔒 Type Safety Features<br/>Fugue's Type System Advantages"]
    B --> C["🔍 Trace Manipulation<br/>Runtime System & Custom Inference"]
    
    A --> D["Statistical Foundations"]
    B --> E["Type-Safe Programming"] 
    C --> F["Advanced Inference"]
    
    D --> G["Production Ready<br/>Probabilistic Models"]
    E --> G
    F --> G
    
    style A fill:#e1f5fe
    style B fill:#f3e5f5  
    style C fill:#e8f5e8
    style G fill:#fff3e0

Learning Path

🎯 Recommended Order

1. Bayesian Coin Flip (~45 minutes)

   • Start here for statistical foundations
   • Learn Bayesian inference principles
   • Understand model specification and analysis

2. Type Safety Features (~30 minutes)

   • Discover Fugue's unique advantages
   • Master type-safe probabilistic programming
   • Eliminate runtime errors with compile-time guarantees

3. Trace Manipulation (~60 minutes)

   • Deep dive into Fugue's runtime system
   • Learn custom inference and debugging techniques
   • Build production-ready probabilistic applications

Tutorial Overview

🪙 Bayesian Coin Flip

Foundation: Statistical inference and model analysis

Your introduction to Bayesian reasoning through the classic coin flipping problem. This tutorial demonstrates how prior beliefs are updated with evidence to form posterior distributions.

Key Concepts:

• Prior, likelihood, and posterior distributions
• Bayesian updating with Beta-Binomial conjugacy
• Model validation and parameter estimation
• Analytical vs computational solutions

What You'll Build:

• Complete Bayesian coin bias estimation model
• Prior sensitivity analysis framework
• Model validation with synthetic data

🔒 Type Safety Features

Foundation: Type-safe probabilistic programming

Explore Fugue's revolutionary type system that eliminates runtime errors while preserving full statistical expressiveness. Learn how dependent types make probabilistic programs both safer and faster.

Key Concepts:

• Natural return types for distributions (bool, u64, f64, usize, i64)
• Compile-time safety guarantees
• Safe array indexing with categorical distributions
• Parameter validation at construction time
• Performance benefits through zero-cost abstractions

What You'll Build:

• Type-safe hierarchical models
• Safe array indexing examples
• Performance comparison with traditional PPLs

🔍 Trace Manipulation

Foundation: Runtime system and advanced inference

Master Fugue's execution trace system — the foundation that enables sophisticated inference algorithms. Learn how traces record, replay, and analyze probabilistic model executions.

Key Concepts:

• Trace system architecture and execution history
• Handler system for flexible model interpretation
• Replay mechanics for MCMC and inference algorithms
• Custom handlers for specialized inference strategies
• Memory optimization for production deployment
• Diagnostic tools for convergence assessment

What You'll Build:

• Custom MCMC algorithm using trace replay
• Specialized handlers for debugging models
• Production inference pipeline with memory optimization
• Comprehensive diagnostic system for model validation

Learning Outcomes

After completing these foundation tutorials, you will:

📊 Statistical Mastery

• ✅ Understand Bayesian inference from first principles
• ✅ Build and validate probabilistic models confidently
• ✅ Interpret posterior distributions and uncertainty quantification
• ✅ Apply conjugate analysis and computational methods

🛡️ Type-Safe Programming

• ✅ Write probabilistic programs that catch errors at compile time
• ✅ Leverage natural return types for cleaner, safer code
• ✅ Understand performance benefits of zero-cost abstractions
• ✅ Build complex models with guaranteed type safety

⚙️ Advanced Inference

• ✅ Manipulate execution traces for custom inference algorithms
• ✅ Implement specialized handlers for unique requirements
• ✅ Debug and optimize problematic models systematically
• ✅ Deploy production-ready probabilistic systems

🔧 Production Skills

• ✅ Memory optimization techniques for high-throughput scenarios
• ✅ Convergence diagnostics and model validation workflows
• ✅ Custom inference algorithms tailored to specific problems
• ✅ Systematic debugging of numerical issues

Code Examples

All tutorials include comprehensive, tested code examples:

• 📁 examples/bayesian_coin_flip.rs - Complete Bayesian analysis
• 📁 examples/type_safety.rs - Type system demonstrations
• 📁 examples/trace_manipulation.rs - Runtime system examples

Each example includes:

• ✅ Comprehensive tests ensuring correctness
• ✅ Detailed comments explaining every concept
• ✅ Runnable code you can execute immediately
• ✅ Performance benchmarks where applicable

# Run any example to see concepts in action
cargo run --example bayesian_coin_flip
cargo run --example type_safety  
cargo run --example trace_manipulation

# Run tests to verify your understanding
cargo test --example bayesian_coin_flip
cargo test --example type_safety
cargo test --example trace_manipulation

Next Steps

After mastering these foundations, you're ready for:

📈 Statistical Modeling Tutorials

Apply your knowledge to real-world problems:

• Linear and logistic regression
• Hierarchical models and mixed effects
• Mixture models and clustering
• Time series and forecasting

🏗️ How-To Guides

Practical guidance for specific tasks:

• Building Complex Models
• Custom Handlers
• Optimizing Performance
• Debugging Models

🚀 Advanced Inference

Inference methods beyond adaptive MCMC, all demonstrated on one worked example so you can compare them directly:

• Sequential Monte Carlo (particle filtering with a log-evidence estimate)
• Approximate Bayesian Computation (likelihood-free inference)
• Variational Inference (fast, optimization-based approximation)

Getting Help

📚 Documentation

• Getting Started Guide - Fugue basics
• API Reference - Complete function documentation
• How-To Guides - Task-specific instructions

💡 Tips for Success

🔧 Troubleshooting

If you encounter issues:

1. Check prerequisites - Ensure you have the required mathematical background
2. Run examples step-by-step - Isolate where confusion arises
3. Review error messages - Fugue's type system provides helpful compile-time feedback
4. Consult diagnostics - Use trace analysis to debug model behavior

Ready to Begin?

Start your journey with Bayesian Coin Flip — the gateway to mastering probabilistic programming with Fugue.

Bayesian Coin Flip