//! Bessel's Correction Test
//!
//! Validates that variance calculations use the unbiased estimator (N-1 denominator)
//! instead of the biased estimator (N denominator).
//!
//! Background:
//! - Biased variance: σ² = Σ(x - μ)² / N
//! - Unbiased variance: s² = Σ(x - μ)² / (N-1)  [Bessel's correction]
//!
//! Bessel's correction compensates for using the sample mean instead of the
//! population mean, providing an unbiased estimate of the population variance.

use approx::assert_relative_eq;

/// Manual calculation of variance with Bessel's correction
fn calculate_variance_unbiased(data: &[f32]) -> f32 {
    if data.len() < 2 {
        return 0.0; // Variance undefined for N=1
    }

    let mean: f32 = data.iter().sum::<f32>() / data.len() as f32;
    let sum_sq: f32 = data.iter().map(|&x| (x - mean).powi(2)).sum();
    sum_sq / (data.len() - 1) as f32 // N-1 (unbiased)
}

/// Manual calculation of variance without Bessel's correction (biased)
fn calculate_variance_biased(data: &[f32]) -> f32 {
    let mean: f32 = data.iter().sum::<f32>() / data.len() as f32;
    let sum_sq: f32 = data.iter().map(|&x| (x - mean).powi(2)).sum();
    sum_sq / data.len() as f32 // N (biased)
}

#[test]
fn test_windowed_variance_uses_bessel_correction() {
    // Given: Small window with known values
    let window = vec![1.0f32, 2.0, 3.0, 4.0, 5.0]; // N=5

    // When: Calculate variance (manual computation)
    let variance_unbiased = calculate_variance_unbiased(&window);
    let variance_biased = calculate_variance_biased(&window);

    // Then: Verify mathematical correctness
    // Mean = 3.0
    // Sum of squared deviations = (1-3)^2 + (2-3)^2 + (3-3)^2 + (4-3)^2 + (5-3)^2 = 10
    // Unbiased variance = 10 / (5-1) = 2.5
    // Biased variance = 10 / 5 = 2.0 (WRONG)

    assert_relative_eq!(variance_unbiased, 2.5, epsilon = 1e-5);
    assert_relative_eq!(variance_biased, 2.0, epsilon = 1e-5);

    // Verify they are different
    assert!((variance_unbiased - variance_biased).abs() > 0.4);
}

#[test]
fn test_bessel_correction_impact() {
    // Given: Realistic price window
    let prices = vec![5000.0f32, 5010.0, 5020.0, 5030.0, 5040.0];

    // When: Calculate std with and without Bessel's correction
    let std_biased = calculate_variance_biased(&prices).sqrt();
    let std_unbiased = calculate_variance_unbiased(&prices).sqrt();

    // Then: Unbiased should be ~sqrt(N/(N-1)) ≈ 1.118x larger
    let ratio = std_unbiased / std_biased;
    assert_relative_eq!(ratio, 1.118, epsilon = 0.01); // sqrt(5/4) = 1.118

    // Verify magnitudes make sense
    assert!(std_unbiased > std_biased);
    assert!(std_unbiased > 15.0); // Should be ~15.81
    assert!(std_biased > 14.0); // Should be ~14.14
}

#[test]
fn test_bessel_correction_edge_case_n1() {
    // Given: Single element window
    let single = vec![42.0f32];

    // When: Calculate variance
    let variance = calculate_variance_unbiased(&single);

    // Then: Should return 0.0 (variance undefined for N=1)
    assert_eq!(variance, 0.0);
}

#[test]
fn test_bessel_correction_edge_case_n2() {
    // Given: Two element window
    let pair = vec![10.0f32, 20.0];

    // When: Calculate variance
    let variance_unbiased = calculate_variance_unbiased(&pair);
    let variance_biased = calculate_variance_biased(&pair);

    // Then: Unbiased uses N-1=1, biased uses N=2
    // Mean = 15.0
    // Sum of squared deviations = (10-15)^2 + (20-15)^2 = 50
    // Unbiased variance = 50 / 1 = 50.0
    // Biased variance = 50 / 2 = 25.0

    assert_relative_eq!(variance_unbiased, 50.0, epsilon = 1e-5);
    assert_relative_eq!(variance_biased, 25.0, epsilon = 1e-5);

    // N=2 shows maximum impact: 2x difference
    let ratio = variance_unbiased / variance_biased;
    assert_relative_eq!(ratio, 2.0, epsilon = 1e-5);
}

#[test]
fn test_bessel_correction_converges_large_n() {
    // Given: Large window (N=1000)
    let large_window: Vec<f32> = (0..1000).map(|i| i as f32).collect();

    // When: Calculate variance
    let variance_unbiased = calculate_variance_unbiased(&large_window);
    let variance_biased = calculate_variance_biased(&large_window);

    // Then: For large N, difference should be small (~0.1%)
    let ratio = variance_unbiased / variance_biased;
    assert_relative_eq!(ratio, 1.001, epsilon = 0.001); // 1000/999 ≈ 1.001

    // But they should still be different
    assert!(variance_unbiased > variance_biased);
}

#[test]
fn test_bessel_correction_constant_values() {
    // Given: Constant window (zero variance)
    let constant = vec![42.0f32, 42.0, 42.0, 42.0];

    // When: Calculate variance
    let variance_unbiased = calculate_variance_unbiased(&constant);
    let variance_biased = calculate_variance_biased(&constant);

    // Then: Both should be 0.0 (no variation)
    assert_eq!(variance_unbiased, 0.0);
    assert_eq!(variance_biased, 0.0);
}