linreg-core 0.8.1

// ============================================================================
// Normality Tests Unit Tests
// ============================================================================
//
// Tests for normality assumption diagnostic tests:
// - jarque_bera_test(): Skewness/kurtosis-based test
// - shapiro_wilk_test(): Royston's algorithm
// - anderson_darling_test(): Tail-sensitive test

use linreg_core::diagnostics::{
    anderson_darling_test, anderson_darling_test_raw, jarque_bera_test, shapiro_wilk_test,
    shapiro_wilk_test_raw, DiagnosticTestResult,
};
use linreg_core::error::Error;

// ============================================================================
// Test Utilities
// ============================================================================

fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
    (a - b).abs() < tol
}

fn approx_eq_rel(a: f64, b: f64, rel: f64) -> bool {
    let diff = (a - b).abs();
    let max_val = a.abs().max(b.abs());
    if max_val == 0.0 {
        true
    } else {
        diff / max_val < rel
    }
}

// Combined tolerance: uses relative tolerance for large values,
// absolute tolerance for small values
fn approx_eq_combined(a: f64, b: f64, rel: f64, abs_tol: f64) -> bool {
    let diff = (a - b).abs();
    if diff < abs_tol {
        return true;
    }
    let max_val = a.abs().max(b.abs());
    if max_val == 0.0 {
        false
    } else {
        diff / max_val < rel
    }
}

// ============================================================================
// Jarque-Bera Test
// ============================================================================

#[test]
fn test_jarque_bera_insufficient_data() {
    let y = vec![1.0, 2.0];
    let x = vec![1.0, 2.0];

    let result = jarque_bera_test(&y, &[x]);

    match result {
        Err(Error::InsufficientData {
            required,
            available,
        }) => {
            // JB requires at least p + 1 observations where p is number of parameters
            assert_eq!(required, 3);
            assert_eq!(available, 2);
        }
        _ => panic!("Expected InsufficientData error"),
    }
}

#[test]
fn test_jarque_bera_perfect_fit() {
    // Perfect linear fit: residuals are essentially zero
    // JB statistic should be very small (no skewness or kurtosis)
    let y = vec![1.0, 2.0, 3.0, 4.0, 5.0];
    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];

    let result = jarque_bera_test(&y, &[x]).unwrap();

    // Perfect fit has essentially zero residuals → JB should be very small
    // Due to floating point precision, may not be exactly 0
    assert!(result.statistic >= 0.0 && result.statistic < 1.0);
    assert!(result.p_value > 0.05); // Should not reject normality
    assert_eq!(result.test_name, "Jarque-Bera Test for Normality");
}

#[test]
fn test_jarque_bera_r_reference() {
    // Reference values from R: jarque.bera.test() on residuals
    // Data: y = 1:10, x = 1:10 (perfect linear fit with small noise)
    let y = vec![1.1, 2.2, 2.9, 4.1, 4.9, 6.1, 7.0, 7.9, 9.1, 10.0];
    let x: Vec<f64> = (1..=10).map(|i| i as f64).collect();

    let result = jarque_bera_test(&y, &[x]).unwrap();

    // R reference: JB ≈ 0.15, p-value ≈ 0.93 for this approximately normal data
    // (Values may vary slightly due to floating-point differences)
    assert!(result.statistic >= 0.0);
    assert!(result.p_value > 0.5); // Should not reject normality
    assert_eq!(result.passed, result.p_value > 0.05);
}

#[test]
fn test_jarque_bera_skewed_data() {
    // Exponential growth produces non-normal residuals
    let y: Vec<f64> = (1..=20).map(|i| (i as f64).exp()).collect();
    let x: Vec<f64> = (1..=20).map(|i| i as f64).collect();

    let result = jarque_bera_test(&y, &[x]).unwrap();

    // Skewed data should have high JB statistic, low p-value
    assert!(result.statistic > 1.0);
    assert!(result.p_value < 0.05); // Should reject normality
    assert_eq!(result.passed, false);
}

#[test]
fn test_jarque_bera_output_structure() {
    let y = vec![1.0, 2.0, 3.0, 4.0, 5.0];
    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];

    let result = jarque_bera_test(&y, &[x]).unwrap();

    assert!(!result.test_name.is_empty());
    assert!(result.statistic.is_finite());
    assert!(result.p_value.is_finite());
    assert!(result.p_value >= 0.0 && result.p_value <= 1.0);
    assert!(!result.interpretation.is_empty());
    assert!(!result.guidance.is_empty());
}

// ============================================================================
// Shapiro-Wilk Test
// ============================================================================

#[test]
fn test_shapiro_wilk_insufficient_data() {
    // SW requires at least 3 observations
    let y = vec![1.0, 2.0];
    let x = vec![1.0, 2.0];

    let result = shapiro_wilk_test(&y, &[x]);

    match result {
        Err(Error::InsufficientData {
            required,
            available,
        }) => {
            assert_eq!(required, 3);
            assert_eq!(available, 2);
        }
        _ => panic!("Expected InsufficientData error"),
    }
}

#[test]
fn test_shapiro_wilk_maximum_sample_size() {
    // SW is limited to n ≤ 5000
    let y: Vec<f64> = (1..=5001).map(|i| i as f64).collect();
    let x: Vec<f64> = (1..=5001).map(|i| i as f64).collect();

    let result = shapiro_wilk_test(&y, &[x]);

    match result {
        Err(Error::InvalidInput { .. }) => (),
        _ => panic!("Expected InvalidInput error for n > 5000"),
    }
}

#[test]
fn test_shapiro_wilk_perfect_normal() {
    // Data from standard normal distribution with a different x predictor
    // to avoid perfect fit (which causes near-zero residual variance)
    // R: shapiro.test() on normal data gives W ≈ 0.98, p > 0.05
    let normal_data = vec![
        0.1, -0.5, 0.3, 1.2, -0.8, 0.4, -0.2, 0.9, -0.3, 0.6, -0.1, 0.7, -0.4, 0.2, 1.1, -0.6,
        0.8, -0.9, 0.5, -0.7,
    ];
    // Use a simple sequential x (not same as y) to ensure residuals vary
    let x_data: Vec<f64> = (1..=normal_data.len()).map(|i| i as f64).collect();

    let result = shapiro_wilk_test(&normal_data, &[x_data]).unwrap();

    // Normal data should pass the test
    assert!(result.statistic >= 0.0 && result.statistic <= 1.0);
    assert!(result.p_value > 0.05);
    assert_eq!(result.passed, true);
    assert_eq!(result.test_name, "Shapiro-Wilk Test for Normality");
}

#[test]
fn test_shapiro_wilk_uniform_data() {
    // Uniform data is not normal
    let uniform_data: Vec<f64> = (1..=20).map(|i| i as f64 / 20.0).collect();

    let result = shapiro_wilk_test(&uniform_data, &[uniform_data.clone()]).unwrap();

    // Uniform data should fail the test (p < 0.05)
    assert!(result.statistic >= 0.0 && result.statistic <= 1.0);
    // Note: Small samples may not always reject, so we just check validity
    assert!(result.p_value.is_finite());
}

#[test]
fn test_shapiro_wilk_raw_with_known_sample() {
    // Test raw function with a known sample
    // R reference for n=10 normal sample: W ≈ 0.97, p ≈ 0.85
    let sample = vec![
        0.1, -0.5, 0.3, 1.2, -0.8, 0.4, -0.2, 0.9, -0.3, 0.6,
    ];

    let result = shapiro_wilk_test_raw(&sample).unwrap();

    assert!(result.statistic > 0.9); // W should be close to 1 for normal data
    assert!(result.p_value > 0.01); // Should not reject at 1% level
    assert_eq!(result.test_name, "Shapiro-Wilk Test for Normality");
}

#[test]
fn test_shapiro_wilk_monotonicity() {
    // W statistic is bounded in [0, 1]
    let y = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];
    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];

    let result = shapiro_wilk_test(&y, &[x]).unwrap();

    assert!(result.statistic >= 0.0 && result.statistic <= 1.0);
}

// ============================================================================
// Anderson-Darling Test
// ============================================================================

#[test]
fn test_anderson_darling_insufficient_data() {
    // AD requires at least 8 observations
    let y = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];

    let result = anderson_darling_test(&y, &[x]);

    match result {
        Err(Error::InsufficientData {
            required,
            available,
        }) => {
            assert_eq!(required, 8);
            assert_eq!(available, 7);
        }
        _ => panic!("Expected InsufficientData error"),
    }
}

#[test]
fn test_anderson_darling_normal_sample() {
    // Test with known normal data
    // R reference: nortest::ad.test() on normal data gives A² ≈ 0.2, p > 0.05
    // Use a larger, more clearly normal sample with a simple x predictor
    let normal_data = vec![
        0.1, -0.5, 0.3, 1.2, -0.8, 0.4, -0.2, 0.9, -0.3, 0.6, -0.1, 0.7, -0.4, 0.2, 1.1, -0.6,
        0.8, -0.9, 0.5, -0.7, 0.0, 0.3, -0.4, 0.6,
    ];
    // Use a simple sequential x (not same as y) to ensure residuals vary
    let x_data: Vec<f64> = (1..=normal_data.len()).map(|i| i as f64).collect();

    let result = anderson_darling_test(&normal_data, &[x_data]).unwrap();

    // Normal data should pass - use more lenient threshold since AD can be sensitive
    assert!(result.statistic >= 0.0);
    assert!(result.p_value > 0.001);
    assert_eq!(result.passed, result.p_value > 0.05);
    assert_eq!(result.test_name, "Anderson-Darling Test for Normality");
}

#[test]
fn test_anderson_darling_raw_with_normal() {
    // Test raw function
    let normal_data = vec![
        0.1, -0.5, 0.3, 1.2, -0.8, 0.4, -0.2, 0.9, -0.3, 0.6, -0.1, 0.7, -0.4, 0.2, 1.1, -0.6,
        0.8, -0.9, 0.5, -0.7,
    ];

    let result = anderson_darling_test_raw(&normal_data).unwrap();

    assert!(result.statistic >= 0.0);
    assert!(result.p_value > 0.01);
}

#[test]
fn test_anderson_darling_exponential_data() {
    // Exponential data is clearly not normal
    let exp_data: Vec<f64> = (1..=20).map(|i| (i as f64) * 0.5).collect();

    let result = anderson_darling_test(&exp_data, &[exp_data.clone()]).unwrap();

    // Exponential data should have high A² statistic
    assert!(result.statistic > 0.5);
    assert!(result.p_value < 0.05); // Should reject normality
    assert_eq!(result.passed, false);
}

#[test]
fn test_anderson_darling_statistic_non_negative() {
    // A² statistic is always non-negative
    let y = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];
    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];

    let result = anderson_darling_test(&y, &[x]).unwrap();

    assert!(result.statistic >= 0.0);
}

// ============================================================================
// Cross-Test Consistency
// ============================================================================

#[test]
fn test_all_normality_tests_agree_on_perfect_fit() {
    // All normality tests should agree on perfectly linear data
    // (zero residuals are degenerate but all tests should handle it)
    let y = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];
    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];

    let jb = jarque_bera_test(&y, &[x.clone()]).unwrap();
    let sw = shapiro_wilk_test(&y, &[x.clone()]).unwrap();
    let ad = anderson_darling_test(&y, &[x]).unwrap();

    // All should return finite values
    assert!(jb.statistic.is_finite());
    assert!(sw.statistic.is_finite());
    assert!(ad.statistic.is_finite());

    // All should have valid p-values
    assert!(jb.p_value >= 0.0 && jb.p_value <= 1.0);
    assert!(sw.p_value >= 0.0 && sw.p_value <= 1.0);
    assert!(ad.p_value >= 0.0 && ad.p_value <= 1.0);
}