Add MOM estimator options (#30)

* Add MOM estimator options

* Fix repeated SPOT fit tail reset

Author: shenxiangzhuang

crates/libspot-rs/README.md

@@ -36,6 +36,32 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
 
 ## Features
 
+### Estimator and threshold options
+
+`SpotDetector::new(config)` keeps the default behavior: it uses `SpotEstimator::Best`, the P2 initial threshold, and the historical `>=` excess update rule.
+
+To reproduce algorithms that explicitly use FluxEV-style MOM-SPOT, configure all options explicitly:
+
+```rust,ignore
+use libspot_rs::{
+    SpotConfig, SpotDetector, SpotEstimator, SpotExcessUpdate, SpotInitialThreshold,
+};
+
+let config = SpotConfig {
+    q: 0.001,
+    level: 0.98,
+    max_excess: 10_000,
+    ..SpotConfig::default()
+};
+
+let mut detector = SpotDetector::new_with_full_options(
+    config,
+    SpotEstimator::Mom,
+    SpotInitialThreshold::Empirical,
+    SpotExcessUpdate::Greater,
+)?;
+```
+
 ### Serialization (Model Persistence)
 
 Serialization support is **enabled by default**. SPOT detectors can be serialized and deserialized for model deployment:

crates/libspot-rs/src/config.rs

@@ -1,5 +1,45 @@
 //! Configuration types for SPOT detector
 
+/// GPD parameter estimator used by SPOT.
+///
+/// The default keeps the historical libspot-rs behavior: try the supported
+/// estimators and keep the fit with the best log-likelihood.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
+pub enum SpotEstimator {
+    /// Try the available estimators and keep the best log-likelihood.
+    #[default]
+    Best,
+    /// Force Method of Moments estimation.
+    Mom,
+}
+
+/// Initial excess-threshold selection strategy.
+///
+/// The default keeps the historical libspot-rs behavior.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
+pub enum SpotInitialThreshold {
+    /// Use the P2 streaming quantile estimator.
+    #[default]
+    P2,
+    /// Use the empirical sorted quantile from the initial batch.
+    Empirical,
+}
+
+/// Streaming excess update condition.
+///
+/// The default keeps the historical libspot-rs behavior.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
+pub enum SpotExcessUpdate {
+    /// Update the tail when `value >= excess_threshold`.
+    #[default]
+    GreaterOrEqual,
+    /// Update the tail only when `value > excess_threshold`.
+    Greater,
+}
+
 /// Configuration parameters for SPOT detector
 ///
 /// # Serialization
@@ -72,4 +112,22 @@ mod tests {
         assert_relative_eq!(config1.level, config2.level);
         assert_eq!(config1.max_excess, config2.max_excess);
     }
+
+    #[test]
+    fn test_spot_estimator_default_keeps_existing_behavior() {
+        assert_eq!(SpotEstimator::default(), SpotEstimator::Best);
+    }
+
+    #[test]
+    fn test_spot_initial_threshold_default_keeps_existing_behavior() {
+        assert_eq!(SpotInitialThreshold::default(), SpotInitialThreshold::P2);
+    }
+
+    #[test]
+    fn test_spot_excess_update_default_keeps_existing_behavior() {
+        assert_eq!(
+            SpotExcessUpdate::default(),
+            SpotExcessUpdate::GreaterOrEqual
+        );
+    }
 }

crates/libspot-rs/src/estimator.rs

@@ -30,6 +30,34 @@ pub fn mom_estimator(peaks: &Peaks) -> (f64, f64, f64) {
     (gamma, sigma, log_likelihood)
 }
 
+/// Method of Moments estimator using sample variance.
+///
+/// This matches the MOM estimator used by FluxEV's SPOT integration.
+pub fn mom_sample_variance_estimator(peaks: &Peaks) -> (f64, f64, f64) {
+    let nt_local = peaks.size();
+    if nt_local < 2 {
+        return (f64::NAN, f64::NAN, f64::NAN);
+    }
+
+    let e = peaks.mean();
+    let mut sum_sq_delta = 0.0;
+    for &value in peaks.container().raw_data().iter().take(nt_local) {
+        let delta = value - e;
+        sum_sq_delta += delta * delta;
+    }
+    let v = sum_sq_delta / (nt_local as f64 - 1.0);
+
+    if e.is_nan() || v.is_nan() || v <= 0.0 {
+        return (f64::NAN, f64::NAN, f64::NAN);
+    }
+
+    let r = e * e / v;
+    let gamma = 0.5 * (1.0 - r);
+    let sigma = 0.5 * e * (1.0 + r);
+
+    (gamma, sigma, 0.0)
+}
+
 /// Grimshaw estimator for GPD parameters
 pub fn grimshaw_estimator(peaks: &Peaks) -> (f64, f64, f64) {
     let mini = peaks.min();
@@ -299,6 +327,21 @@ mod tests {
         assert!(sigma > 0.0); // Sigma should be positive
     }
 
+    #[test]
+    fn test_mom_sample_variance_estimator_normal_case() {
+        let mut peaks = Peaks::new(10).unwrap();
+        for value in [1.0, 2.0, 3.0, 4.0] {
+            peaks.push(value);
+        }
+
+        let (gamma, sigma, llhood) = mom_sample_variance_estimator(&peaks);
+
+        assert!(gamma.is_finite());
+        assert!(sigma.is_finite());
+        assert!(sigma > 0.0);
+        assert_relative_eq!(llhood, 0.0);
+    }
+
     #[test]
     fn test_log_likelihood_gamma_zero() {
         let mut peaks = Peaks::new(10).unwrap();

crates/libspot-rs/src/lib.rs

@@ -59,7 +59,7 @@ mod tail;
 mod ubend;
 
 // Re-export public types
-pub use config::SpotConfig;
+pub use config::{SpotConfig, SpotEstimator, SpotExcessUpdate, SpotInitialThreshold};
 pub use error::{SpotError, SpotResult};
 pub use peaks::Peaks;
 pub use spot::SpotDetector;