r2rs-stats 0.1.1

// "Whatever you do, work at it with all your heart, as working for the Lord,
// not for human masters, since you know that you will receive an inheritance
// from the Lord as a reward. It is the Lord Christ you are serving."
// (Col 3:23-24)

use std::str::FromStr;

use nalgebra::DMatrix;
use num_traits::ToPrimitive;
use strafe_testing::{
    r::{RGenerator, RString, RTester},
    r_assert_relative_equal, r_assert_relative_equal_result,
};
use strafe_trait::{Model, ModelBuilder, Statistic, StatisticalTest};
use strafe_type::ModelMatrix;

use crate::regression::lm::LeastSquaresRegressionBuilder;

pub fn linear_regression_test_inner(seed: u16, num_features: f64, num_rows: f64, fuzz_amount: f64) {
    let num_features = num_features.to_usize().unwrap();
    let num_rows = num_rows.to_usize().unwrap();
    let fuzz_amount = fuzz_amount.to_usize().unwrap();

    let xs = RGenerator::new().set_seed(seed).generate_uniform_matrix(
        num_rows,
        num_features + 1,
        (-10, 10),
    );

    let features = RGenerator::new()
        .set_seed(seed + 1)
        .generate_uniform(num_features + 1, (-100, 100));

    let errors =
        RGenerator::new()
            .set_seed(seed + 2)
            .generate_errors(num_rows, (0, 1), fuzz_amount);

    let ys = (0..num_rows)
        .map(|i| {
            (0..num_features).fold(features[0] + errors[i], |ret, j| {
                ret + (features[j + 1] * xs[(i, j)])
            })
        })
        .collect::<Vec<_>>();

    let y = DMatrix::from_iterator(ys.len(), 1, ys.iter().copied());

    let mut model = LeastSquaresRegressionBuilder::new()
        .with_x(&ModelMatrix::from(&xs))
        .with_y(&ModelMatrix::from(&y))
        .build();
    let model_test = model.test(&()).unwrap();

    let r_regression = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};",
            RString::from_f64_slice(&ys),
            RString::from_f64_matrix(&xs)
        )))
        .set_display(&RString::from_str("summary(lm(y~x))").unwrap())
        .run()
        .expect("R error")
        .as_lm_regression();

    r_assert_relative_equal_result!(
        Ok(r_regression.f),
        &model_test.significance_of_regression.statistic(),
        1e-4
    );
    if r_regression.p <= 0.05 {
        if model_test.significance_of_regression.probability_value() > 0.05 {
            r_assert_relative_equal_result!(
                Ok(r_regression.p),
                &model_test.significance_of_regression.probability_value(),
                1e-4
            );
        }
    } else {
        if model_test.significance_of_regression.probability_value() <= 0.05 {
            r_assert_relative_equal_result!(
                Ok(r_regression.p),
                &model_test.significance_of_regression.probability_value(),
                1e-4
            );
        }
    }

    for i in 0..num_features {
        assert!(
            r_regression.coefs[i].estimate
                >= model.parameters().unwrap()[i].confidence_interval().0
        );
        assert!(
            r_regression.coefs[i].estimate
                <= model.parameters().unwrap()[i].confidence_interval().1
        );
        r_assert_relative_equal_result!(
            Ok(r_regression.coefs[i].estimate),
            &model.parameters().unwrap()[i].estimate(),
            1e-4
        );
        r_assert_relative_equal_result!(
            Ok(r_regression.coefs[i].estimate),
            &model.parameters().unwrap()[i].estimate(),
            1e-4
        );
        r_assert_relative_equal_result!(
            Ok(r_regression.coefs[i].t),
            &model_test.significance_of_coefficients[i].statistic(),
            1e-4
        );
        if r_regression.coefs[i].p <= 0.05 {
            if model_test.significance_of_coefficients[i].probability_value() > 0.05 {
                r_assert_relative_equal_result!(
                    Ok(r_regression.coefs[i].p),
                    &model_test.significance_of_coefficients[i].probability_value(),
                    1e-4
                );
            }
        } else {
            if model_test.significance_of_coefficients[i].probability_value() <= 0.05 {
                r_assert_relative_equal_result!(
                    Ok(r_regression.coefs[i].p),
                    &model_test.significance_of_coefficients[i].probability_value(),
                    1e-4
                );
            }
        }
    }

    r_assert_relative_equal_result!(
        Ok(r_regression.r_squared),
        &model.determination().unwrap().statistic(),
        1e-3
    );

    let r_conf = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};",
            RString::from_f64_slice(&ys),
            RString::from_f64_matrix(&xs)
        )))
        .set_display(&RString::from_str("predict(lm(y~x), interval='confidence')").unwrap())
        .run()
        .expect("R error")
        .as_f64_matrix()
        .expect("Cannot get R matrix");

    for (&r_lower_conf, &rust_lower_conf) in r_conf
        .column(0)
        .iter()
        .zip(model.predictions().unwrap().estimate().matrix().iter())
    {
        r_assert_relative_equal!(r_lower_conf, rust_lower_conf, 1e-4);
    }
    for (&r_lower_conf, &rust_lower_conf) in r_conf.column(1).iter().zip(
        model
            .predictions()
            .unwrap()
            .confidence_interval()
            .matrix()
            .column(0)
            .iter(),
    ) {
        r_assert_relative_equal!(r_lower_conf, rust_lower_conf, 1e-4);
    }
    for (&r_upper_conf, &rust_upper_conf) in r_conf.column(2).iter().zip(
        model
            .predictions()
            .unwrap()
            .confidence_interval()
            .matrix()
            .column(1)
            .iter(),
    ) {
        r_assert_relative_equal!(r_upper_conf, rust_upper_conf, 1e-4);
    }

    let r_pred = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};",
            RString::from_f64_slice(&ys),
            RString::from_f64_matrix(&xs)
        )))
        .set_display(&RString::from_str("predict(lm(y~x), interval='prediction')").unwrap())
        .run()
        .expect("R error")
        .as_f64_matrix()
        .expect("Cannot get R matrix");

    for (&r_lower_conf, &rust_lower_conf) in r_pred.column(1).iter().zip(
        model
            .manipulator()
            .unwrap()
            .confidence_interval()
            .matrix()
            .column(0)
            .iter(),
    ) {
        r_assert_relative_equal!(r_lower_conf, rust_lower_conf, 1e-4);
    }
    for (&r_upper_conf, &rust_upper_conf) in r_pred.column(2).iter().zip(
        model
            .manipulator()
            .unwrap()
            .confidence_interval()
            .matrix()
            .column(1)
            .iter(),
    ) {
        r_assert_relative_equal!(r_upper_conf, rust_upper_conf, 1e-4);
    }

    let r_res = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};",
            RString::from_f64_slice(&ys),
            RString::from_f64_matrix(&xs),
        )))
        .set_display(&RString::from_str("unname(residuals(lm(y~x)))").unwrap())
        .run()
        .expect("R error")
        .as_f64_vec()
        .expect("R error");

    for (&rust_ret, &r_ret) in model.residuals().unwrap().matrix().iter().zip(r_res.iter()) {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-4);
    }

    let r_vcov = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};",
            RString::from_f64_slice(&ys),
            RString::from_f64_matrix(&xs),
        )))
        .set_display(&RString::from_str("vcov(lm(y~x, TAU='R'))").unwrap())
        .run()
        .expect("R error")
        .as_f64_matrix()
        .expect("R error");

    for (&rust_ret, &r_ret) in model.variance().unwrap().iter().zip(r_vcov.iter()) {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-4);
    }

    let r_res = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};",
            RString::from_f64_slice(&ys),
            RString::from_f64_matrix(&xs),
        )))
        .set_display(&RString::from_str("unname(rstudent(lm(y~x)))").unwrap())
        .run()
        .expect("R error")
        .as_f64_vec()
        .expect("R error");

    for (&rust_ret, &r_ret) in model
        .studentized_residuals()
        .unwrap()
        .matrix()
        .iter()
        .zip(r_res.iter())
    {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-1);
    }

    let r_res = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};",
            RString::from_f64_slice(&ys),
            RString::from_f64_matrix(&xs),
        )))
        .set_display(&RString::from_str("unname(rstandard(lm(y~x)))").unwrap())
        .run()
        .expect("R error")
        .as_f64_vec()
        .expect("R error");

    for (&rust_ret, &r_ret) in model
        .standardized_residuals()
        .unwrap()
        .matrix()
        .iter()
        .zip(r_res.iter())
    {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-2);
    }
}