341a2cadf0
Adds curiosity-driven exploration (ForwardDynamicsModel + ICM reward), configurable branching DQN fields (num_order_types, num_urgency_levels), regime-conditional importance sampling with ADX/CUSUM thresholds, and GPU curiosity training kernel support. Also fixes remaining 5 clippy errors from WIP merge: - gpu_smoketest: add 8 missing DQNConfig fields - curiosity.rs: replace needless_range_loop with slice fill - benchmark files: remove redundant #[cfg_attr] on unconditionally ignored tests 40 files changed, +1486/-1068 lines. 0 clippy errors, 0 warnings. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
477 lines
17 KiB
Rust
477 lines
17 KiB
Rust
//! Bayesian Online Changepoint Detection (BOCD)
|
||
//!
|
||
//! This module implements the Bayesian Online Changepoint Detection algorithm
|
||
//! for probabilistic regime change detection in financial time series.
|
||
//!
|
||
//! # Algorithm Overview
|
||
//!
|
||
//! BOCD maintains a distribution over the current "run length" (time since last changepoint)
|
||
//! and updates it online as new data arrives. The algorithm computes:
|
||
//!
|
||
//! 1. **Run-length distribution**: P(rₜ|x₁:ₜ) - probability that the current run length is r
|
||
//! 2. **Hazard function**: H(r) = 1/λ - probability of changepoint given run length r
|
||
//! 3. **Predictive probability**: P(xₜ|x₁:ₜ₋₁, rₜ) - likelihood of observation given run length
|
||
//!
|
||
//! # Mathematical Formulation
|
||
//!
|
||
//! The core update equations are:
|
||
//!
|
||
//! ```text
|
||
//! P(rₜ|x₁:ₜ) ∝ P(xₜ|rₜ, x₁:ₜ₋₁) × [
|
||
//! P(rₜ₋₁ = rₜ - 1|x₁:ₜ₋₁) × (1 - H(rₜ-1)) if rₜ > 0
|
||
//! Σᵣ P(rₜ₋₁ = r|x₁:ₜ₋₁) × H(r) if rₜ = 0
|
||
//! ]
|
||
//! ```
|
||
//!
|
||
//! Where:
|
||
//! - H(r) = 1/λ is the constant hazard function (λ = expected run length)
|
||
//! - P(xₜ|rₜ, x₁:ₜ₋₁) is computed using a conjugate Gaussian model
|
||
//!
|
||
//! # Usage
|
||
//!
|
||
//! ```rust
|
||
//! use ml::regime::bayesian_changepoint::{BayesianChangepointDetector, ChangepointInfo};
|
||
//!
|
||
//! // Create detector with hazard rate λ=100 (expect changepoint every 100 bars)
|
||
//! let mut detector = BayesianChangepointDetector::new(100.0, 0.3, 200);
|
||
//!
|
||
//! // Update with new values
|
||
//! for price in &[100.0, 101.0, 102.0, 150.0] { // Jump at 150.0
|
||
//! if let Some(info) = detector.update(*price) {
|
||
//! println!("Changepoint detected! Probability: {}", info.probability);
|
||
//! }
|
||
//! }
|
||
//!
|
||
//! // Query current state
|
||
//! let prob = detector.get_changepoint_probability();
|
||
//! let run_length = detector.get_expected_run_length();
|
||
//! ```
|
||
//!
|
||
//! # Performance
|
||
//!
|
||
//! - Target: <150μs per update (Bayesian computation intensive)
|
||
//! - Memory: O(max_run_length) for probability distribution
|
||
//! - Online: Constant time per update (no recomputation of history)
|
||
//!
|
||
//! # References
|
||
//!
|
||
//! - Adams & MacKay (2007): "Bayesian Online Changepoint Detection"
|
||
//! - Used in: Regime detection, structural break identification, adaptive strategies
|
||
|
||
use serde::{Deserialize, Serialize};
|
||
|
||
/// Information about a detected changepoint
|
||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||
pub struct ChangepointInfo {
|
||
/// Probability of changepoint (0.0 to 1.0)
|
||
pub probability: f64,
|
||
|
||
/// Expected run length since last changepoint
|
||
pub expected_run_length: f64,
|
||
|
||
/// Maximum a posteriori (MAP) run length
|
||
pub map_run_length: usize,
|
||
|
||
/// Current observation value
|
||
pub value: f64,
|
||
|
||
/// Time index of detection
|
||
pub time_index: usize,
|
||
}
|
||
|
||
/// Bayesian Online Changepoint Detector
|
||
///
|
||
/// Maintains a distribution over the current run length (time since last changepoint)
|
||
/// and updates it online using Bayesian inference.
|
||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||
pub struct BayesianChangepointDetector {
|
||
/// Hazard rate parameter (λ): Expected run length = 1/hazard_rate
|
||
/// Higher λ → more frequent changepoints
|
||
/// Lower λ → less frequent changepoints
|
||
hazard_rate: f64,
|
||
|
||
/// Probability threshold for changepoint detection (0.0 to 1.0)
|
||
/// Typical values: 0.2-0.4
|
||
changepoint_prob_threshold: f64,
|
||
|
||
/// Maximum run length to track (truncation for computational efficiency)
|
||
max_run_length: usize,
|
||
|
||
/// Current run-length probability distribution P(rₜ|x₁:ₜ)
|
||
/// Index r represents run length r (time since last changepoint)
|
||
run_length_probs: Vec<f64>,
|
||
|
||
/// Sufficient statistics for Gaussian model (online updates)
|
||
/// mean`[r]` = mean of observations for run length r
|
||
means: Vec<f64>,
|
||
|
||
/// Sufficient statistics for Gaussian model (online updates)
|
||
/// variance`[r]` = variance of observations for run length r
|
||
variances: Vec<f64>,
|
||
|
||
/// Number of observations for each run length
|
||
counts: Vec<f64>,
|
||
|
||
/// Current time index (number of observations processed)
|
||
time_index: usize,
|
||
|
||
/// Prior hyperparameters for Gaussian model
|
||
/// μ₀: Prior mean
|
||
prior_mean: f64,
|
||
|
||
/// Prior hyperparameters for Gaussian model
|
||
/// κ₀: Prior precision (pseudo-count)
|
||
prior_precision: f64,
|
||
|
||
/// Prior hyperparameters for Gaussian model
|
||
/// α₀: Prior degrees of freedom
|
||
prior_alpha: f64,
|
||
|
||
/// Prior hyperparameters for Gaussian model
|
||
/// β₀: Prior variance scale
|
||
prior_beta: f64,
|
||
}
|
||
|
||
impl BayesianChangepointDetector {
|
||
/// Create a new Bayesian changepoint detector
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `hazard_rate` - Expected run length = 1/hazard_rate (e.g., 100.0 → expect changepoint every 100 bars)
|
||
/// * `threshold` - Probability threshold for changepoint detection (0.0 to 1.0, typical: 0.2-0.4)
|
||
/// * `max_run_length` - Maximum run length to track (truncation, typical: 200-500)
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// New detector instance with default Gaussian priors
|
||
///
|
||
/// # Examples
|
||
///
|
||
/// ```rust
|
||
/// use ml::regime::bayesian_changepoint::BayesianChangepointDetector;
|
||
///
|
||
/// // Expect changepoint every 100 bars, detect at 30% probability, track up to 200 bars
|
||
/// let detector = BayesianChangepointDetector::new(100.0, 0.3, 200);
|
||
/// ```
|
||
pub fn new(hazard_rate: f64, threshold: f64, max_run_length: usize) -> Self {
|
||
// Initialize run-length distribution: P(r₀ = 0) = 1.0 (start with run length 0)
|
||
let mut run_length_probs = vec![0.0; max_run_length + 1];
|
||
run_length_probs[0] = 1.0;
|
||
|
||
// Non-informative priors for Gaussian model
|
||
// prior_beta = 1e4 gives prior variance = β/α = 1e4, broad enough for any
|
||
// data scale. This prevents the changepoint hypothesis from being starved
|
||
// by a narrow prior when the first observations set the mean far from 0.
|
||
let prior_mean = 0.0; // μ₀: No prior knowledge of mean
|
||
let prior_precision = 0.01; // κ₀: Low confidence in prior mean
|
||
let prior_alpha = 1.0; // α₀: Minimal degrees of freedom
|
||
let prior_beta = 1e4; // β₀: Broad variance scale (supports any data scale)
|
||
|
||
Self {
|
||
hazard_rate,
|
||
changepoint_prob_threshold: threshold,
|
||
max_run_length,
|
||
run_length_probs,
|
||
means: vec![prior_mean; max_run_length + 1],
|
||
variances: vec![1.0; max_run_length + 1],
|
||
counts: vec![0.0; max_run_length + 1],
|
||
time_index: 0,
|
||
prior_mean,
|
||
prior_precision,
|
||
prior_alpha,
|
||
prior_beta,
|
||
}
|
||
}
|
||
|
||
/// Update detector with a new observation
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `value` - New observation value
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// Some(ChangepointInfo) if changepoint detected (P(r=0) > threshold), None otherwise
|
||
///
|
||
/// # Algorithm
|
||
///
|
||
/// 1. Compute predictive probabilities P(xₜ|rₜ, x₁:ₜ₋₁) for all run lengths
|
||
/// 2. Update run-length distribution using Bayes' rule and hazard function
|
||
/// 3. Update sufficient statistics for Gaussian model
|
||
/// 4. Check if P(r=0) exceeds threshold
|
||
///
|
||
/// # Examples
|
||
///
|
||
/// ```rust
|
||
/// use ml::regime::bayesian_changepoint::BayesianChangepointDetector;
|
||
///
|
||
/// let mut detector = BayesianChangepointDetector::new(100.0, 0.3, 200);
|
||
///
|
||
/// // Stable regime
|
||
/// for i in 0..50 {
|
||
/// detector.update(100.0 + (i as f64 * 0.1));
|
||
/// }
|
||
///
|
||
/// // Regime change (sudden jump)
|
||
/// if let Some(info) = detector.update(150.0) {
|
||
/// println!("Changepoint detected at probability {}", info.probability);
|
||
/// }
|
||
/// ```
|
||
pub fn update(&mut self, value: f64) -> Option<ChangepointInfo> {
|
||
self.time_index += 1;
|
||
|
||
// Step 1: Compute predictive probabilities P(xₜ|rₜ, x₁:ₜ₋₁) for all run lengths
|
||
let mut predictive_probs = vec![0.0; self.max_run_length + 1];
|
||
for r in 0..=self.max_run_length {
|
||
if self.run_length_probs[r] > 1e-10 {
|
||
// Skip negligible probabilities
|
||
predictive_probs[r] = self.compute_predictive_probability(value, r);
|
||
}
|
||
}
|
||
|
||
// Step 2: Update run-length distribution using Bayes' rule and hazard function
|
||
let mut new_probs = vec![0.0; self.max_run_length + 1];
|
||
|
||
// Growth probabilities: P(rₜ = r+1|x₁:ₜ) from P(rₜ₋₁ = r|x₁:ₜ₋₁)
|
||
for r in 0..self.max_run_length {
|
||
if self.run_length_probs[r] > 1e-10 {
|
||
let survival_prob = 1.0 - self.hazard_function(r);
|
||
new_probs[r + 1] += self.run_length_probs[r] * predictive_probs[r] * survival_prob;
|
||
}
|
||
}
|
||
|
||
// Changepoint probability: P(rₜ = 0|x₁:ₜ) = Σᵣ P(rₜ₋₁ = r|x₁:ₜ₋₁) × H(r) × P(xₜ|rₜ=0)
|
||
let mut changepoint_prob = 0.0;
|
||
for r in 0..=self.max_run_length {
|
||
if self.run_length_probs[r] > 1e-10 {
|
||
changepoint_prob += self.run_length_probs[r] * self.hazard_function(r);
|
||
}
|
||
}
|
||
new_probs[0] = changepoint_prob * predictive_probs[0];
|
||
|
||
// Step 3: Normalize probabilities
|
||
let total: f64 = new_probs.iter().sum();
|
||
if total > 1e-10 {
|
||
for p in &mut new_probs {
|
||
*p /= total;
|
||
}
|
||
} else {
|
||
// Numerical underflow: Reset to initial state
|
||
new_probs = vec![0.0; self.max_run_length + 1];
|
||
new_probs[0] = 1.0;
|
||
}
|
||
|
||
// Step 3b: If a changepoint is detected, reset sufficient statistics
|
||
// so the new regime starts with fresh priors. Without this, the shared
|
||
// Welford statistics carry contamination from the old regime, preventing
|
||
// the algorithm from properly tracking subsequent regime changes.
|
||
if new_probs[0] > self.changepoint_prob_threshold {
|
||
self.counts = vec![0.0; self.max_run_length + 1];
|
||
self.means = vec![self.prior_mean; self.max_run_length + 1];
|
||
self.variances = vec![1.0; self.max_run_length + 1];
|
||
}
|
||
|
||
// Step 4: Update sufficient statistics for Gaussian model
|
||
self.update_sufficient_statistics(value, &new_probs);
|
||
|
||
// Update run-length distribution
|
||
self.run_length_probs = new_probs;
|
||
|
||
// Step 5: Check for changepoint detection
|
||
let cp_prob = self.run_length_probs[0];
|
||
(cp_prob > self.changepoint_prob_threshold).then(|| ChangepointInfo {
|
||
probability: cp_prob,
|
||
expected_run_length: self.get_expected_run_length(),
|
||
map_run_length: self.get_map_run_length(),
|
||
value,
|
||
time_index: self.time_index,
|
||
})
|
||
}
|
||
|
||
/// Get current changepoint probability P(r=0|x₁:ₜ)
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// Probability that a changepoint just occurred (0.0 to 1.0)
|
||
pub fn get_changepoint_probability(&self) -> f64 {
|
||
self.run_length_probs[0]
|
||
}
|
||
|
||
/// Get expected run length E[r|x₁:ₜ] = Σᵣ r × P(r|x₁:ₜ)
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// Expected number of bars since last changepoint
|
||
pub fn get_expected_run_length(&self) -> f64 {
|
||
self.run_length_probs
|
||
.iter()
|
||
.enumerate()
|
||
.map(|(r, &prob)| r as f64 * prob)
|
||
.sum()
|
||
}
|
||
|
||
/// Get maximum a posteriori (MAP) run length
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// Most likely run length (arg max P(r|x₁:ₜ))
|
||
pub fn get_map_run_length(&self) -> usize {
|
||
self.run_length_probs
|
||
.iter()
|
||
.enumerate()
|
||
.max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
|
||
.map(|(r, _)| r)
|
||
.unwrap_or(0)
|
||
}
|
||
|
||
/// Compute hazard function H(r) = 1/λ (constant hazard)
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `run_length` - Current run length r
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// Probability of changepoint given run length r
|
||
fn hazard_function(&self, _run_length: usize) -> f64 {
|
||
// Constant hazard: H(r) = 1/λ for all r
|
||
1.0 / self.hazard_rate
|
||
}
|
||
|
||
/// Compute predictive probability P(xₜ|rₜ, x₁:ₜ₋₁) using Student's t-distribution
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `value` - Observation value xₜ
|
||
/// * `run_length` - Run length r
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// Predictive probability (likelihood of observation)
|
||
///
|
||
/// # Algorithm
|
||
///
|
||
/// Uses conjugate Gaussian model with Normal-Inverse-Gamma priors.
|
||
/// Predictive distribution is Student's t with parameters:
|
||
/// - Location: μᵣ (posterior mean)
|
||
/// - Scale: √(βᵣ(κᵣ+1)/(αᵣκᵣ)) (posterior variance)
|
||
/// - Degrees of freedom: 2αᵣ
|
||
fn compute_predictive_probability(&self, value: f64, run_length: usize) -> f64 {
|
||
// r=0 (changepoint hypothesis): always use the prior predictive distribution.
|
||
// In Adams & MacKay (2007), a changepoint means a NEW regime starts with
|
||
// no evidence — the prior parameters describe the new regime's distribution.
|
||
// Without this, the shared sufficient statistics at r=0 get contaminated
|
||
// by observations from the current regime, preventing detection.
|
||
if run_length == 0 {
|
||
let variance = self.prior_beta / self.prior_alpha;
|
||
return self.gaussian_pdf(value, self.prior_mean, variance);
|
||
}
|
||
|
||
let n = self.counts[run_length];
|
||
|
||
if n < 1e-10 {
|
||
let variance = self.prior_beta / self.prior_alpha;
|
||
return self.gaussian_pdf(value, self.prior_mean, variance);
|
||
}
|
||
|
||
// Compute posterior parameters
|
||
let mean = self.means[run_length];
|
||
let variance = self.variances[run_length];
|
||
|
||
// Posterior precision: κᵣ = κ₀ + n (accumulates confidence with observations)
|
||
let location = mean;
|
||
let kappa_r = self.prior_precision + n;
|
||
let scale_squared = variance * (kappa_r + 1.0) / kappa_r;
|
||
|
||
// Cap at 1.0: gaussian_pdf for narrow distributions can exceed 1.0 at the
|
||
// mode (it's a density, not a probability). In BOCD, growth factor =
|
||
// pred × (1−H). When pred > 1, growth > 1−H per step, causing exponential
|
||
// amplification of hazard-redistributed probability mass across ALL run
|
||
// lengths, drowning the changepoint signal. Capping at 1.0 ensures growth
|
||
// ≤ (1−H) < 1, so hazard-redistributed mass decays geometrically.
|
||
let pred = if n > 10.0 {
|
||
self.gaussian_pdf(value, location, scale_squared)
|
||
} else {
|
||
let df = 2.0 * self.prior_alpha + n;
|
||
self.student_t_pdf(value, location, scale_squared, df)
|
||
};
|
||
pred.min(1.0)
|
||
}
|
||
|
||
/// Gaussian PDF: N(x|μ, σ²)
|
||
///
|
||
/// Uses 1e-300 floor (not 1e-10) to preserve relative probability ordering
|
||
/// between run lengths. A flat 1e-10 floor masks the discrimination between
|
||
/// established regimes (very low prob for outliers) and the prior (moderate
|
||
/// prob for any value).
|
||
fn gaussian_pdf(&self, x: f64, mean: f64, variance: f64) -> f64 {
|
||
if variance < 1e-10 {
|
||
return if (x - mean).abs() < 1e-6 { 1.0 } else { 1e-300 };
|
||
}
|
||
|
||
let diff = x - mean;
|
||
let exponent = -0.5 * diff * diff / variance;
|
||
let normalization = 1.0 / (2.0 * std::f64::consts::PI * variance).sqrt();
|
||
|
||
(normalization * exponent.exp()).max(1e-300)
|
||
}
|
||
|
||
/// Student's t-distribution PDF (approximation)
|
||
fn student_t_pdf(&self, x: f64, location: f64, scale: f64, df: f64) -> f64 {
|
||
if scale < 1e-10 {
|
||
return if (x - location).abs() < 1e-6 {
|
||
1.0
|
||
} else {
|
||
1e-300
|
||
};
|
||
}
|
||
|
||
// Simplified Student's t approximation (sufficient for BOCD)
|
||
let diff = (x - location) / scale.sqrt();
|
||
let factor = 1.0 + diff * diff / df;
|
||
let exponent = -(df + 1.0) / 2.0;
|
||
|
||
(factor.powf(exponent) / scale.sqrt()).max(1e-300)
|
||
}
|
||
|
||
/// Update sufficient statistics for Gaussian model (online updates)
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `value` - New observation
|
||
/// * `probs` - New run-length distribution
|
||
fn update_sufficient_statistics(&mut self, value: f64, probs: &[f64]) {
|
||
// Update statistics for each run length using online formulas
|
||
for r in 0..=self.max_run_length {
|
||
if probs[r] > 1e-10 {
|
||
let prev_count = self.counts[r];
|
||
let prev_mean = self.means[r];
|
||
|
||
// Update count (weighted by probability)
|
||
self.counts[r] = prev_count + probs[r];
|
||
|
||
// Update mean (online Welford's algorithm)
|
||
let delta = value - prev_mean;
|
||
self.means[r] = prev_mean + delta * probs[r] / self.counts[r];
|
||
|
||
// Update variance (online Welford's algorithm)
|
||
let delta2 = value - self.means[r];
|
||
self.variances[r] =
|
||
(prev_count * self.variances[r] + probs[r] * delta * delta2) / self.counts[r];
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Reset detector to initial state
|
||
pub fn reset(&mut self) {
|
||
self.run_length_probs = vec![0.0; self.max_run_length + 1];
|
||
self.run_length_probs[0] = 1.0;
|
||
self.means = vec![self.prior_mean; self.max_run_length + 1];
|
||
self.variances = vec![1.0; self.max_run_length + 1];
|
||
self.counts = vec![0.0; self.max_run_length + 1];
|
||
self.time_index = 0;
|
||
}
|
||
}
|
||
|
||
// Tests are in ml/tests/bayesian_changepoint_test.rs
|