#![no_std]
use core::marker::PhantomData;
use core::ops::{AddAssign, Mul};
use num_traits::{AsPrimitive, ConstZero, Float};
#[derive(Debug, Clone, Default, Copy)]
pub struct Pearson<T, O> {
x: T,
y: T,
xy: T,
x2: T,
y2: T,
len: usize,
output: PhantomData<O>,
}
impl<T, O> Pearson<T, O>
where
T: ConstZero,
{
pub const fn new() -> Self {
Self {
x: T::ZERO,
y: T::ZERO,
xy: T::ZERO,
x2: T::ZERO,
y2: T::ZERO,
len: 0,
output: PhantomData,
}
}
}
impl<T, O> Pearson<T, O>
where
T: AddAssign + Mul<Output = T> + Copy,
{
pub fn record(&mut self, x: T, y: T) {
self.x += x;
self.y += y;
self.xy += x * y;
self.x2 += x * x;
self.y2 += y * y;
self.len += 1;
}
}
impl<T, O> Pearson<T, O>
where
O: Float + 'static,
T: AsPrimitive<O> + ConstZero,
usize: AsPrimitive<O>,
{
pub fn correlate(&mut self) -> O {
let x_prod_y = self.x.as_() * self.y.as_();
let x_prod_x = self.x.as_() * self.x.as_();
let y_prod_y = self.y.as_() * self.y.as_();
let n = self.len.as_();
let numerator = self.xy.as_() - (x_prod_y / n);
let denominator =
((self.x2.as_() - (x_prod_x / n)) * (self.y2.as_() - (y_prod_y / n))).sqrt();
*self = Self::new();
numerator / denominator
}
}
pub fn correlate<T, O>(x: &[T], y: &[T]) -> O
where
O: Float + 'static,
T: AsPrimitive<O> + ConstZero + AddAssign + Mul<Output = T> + Copy,
usize: AsPrimitive<O>,
{
assert_eq!(x.len(), y.len());
let mut p = Pearson::new();
x.iter().zip(y.iter()).for_each(|(&x, &y)| p.record(x, y));
p.correlate()
}
#[test]
fn simple() {
let x1 = [1, 2, 3, 4, 5];
let x2 = [5, 4, 3, 2, 1];
let y = [1, 2, 3, 4, 5];
let c1: f32 = correlate(&x1, &y);
let c2: f32 = correlate(&x2, &y);
assert_eq!(c1, 1.0);
assert_eq!(c2, -1.0);
}