r2rs-stats 0.1.1

Statistics programming for Rust based on R's stats package
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

// "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 nalgebra::DMatrix;

use crate::funcs::total_sum_of_squares;

pub fn linear_regression_standardized(x: &DMatrix<f64>, y: &DMatrix<f64>) -> DMatrix<f64> {
    let mean_y = y.iter().sum::<f64>() / y.shape().0 as f64;
    let column_means = x
        .column_iter()
        .map(|column| column.iter().sum::<f64>() / column.len() as f64)
        .collect::<Vec<_>>();

    let corrected_sum_of_squares = x
        .column_iter()
        .enumerate()
        .map(|(j, column)| {
            column
                .iter()
                .map(|x_ij| (x_ij - column_means[j]).powi(2))
                .sum::<f64>()
        })
        .collect::<Vec<_>>();

    let sst = total_sum_of_squares(&x.to_owned().insert_column(0, 1.0), y);

    let mut w = x.clone();
    w.column_iter_mut().enumerate().for_each(|(j, mut column)| {
        column.iter_mut().enumerate().for_each(|(i, w_ij)| {
            *w_ij = (x[(i, j)] - column_means[j]) / corrected_sum_of_squares[j].sqrt()
        });
    });

    let mut y_0 = y.clone();
    y_0.iter_mut()
        .enumerate()
        .for_each(|(i, y_0_i)| *y_0_i = (y[(i, 0)] - mean_y) / sst.sqrt());

    let w_r = w.clone().qr().r();

    let standardized_betas = w_r.clone().pseudo_inverse(f64::EPSILON).unwrap()
        * w_r.transpose().pseudo_inverse(f64::EPSILON).unwrap()
        * w.transpose()
        * y_0;

    let mut betas = (0..x.shape().1)
        .map(|j| standardized_betas[j] * (sst / corrected_sum_of_squares[j]).sqrt())
        .collect::<Vec<_>>();

    betas.insert(
        0,
        mean_y
            - (0..x.shape().1)
                .map(|j| betas[j] * column_means[j])
                .sum::<f64>(),
    );
    DMatrix::from_iterator(betas.len(), 1, betas.into_iter())
}