stochastic-rs-stats 0.1.3

A Rust library for the statistical analysis of stochastic processes
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

use std::f64::consts::LN_2;

use linreg::linear_regression;
use ndarray::Array1;

/// Fractal dimension.
pub struct FractalDim {
  pub x: Array1<f64>,
}

impl FractalDim {
  #[must_use]
  pub fn new(x: Array1<f64>) -> Self {
    Self { x }
  }

  /// Calculate the variogram of the path.
  pub fn variogram(&self, p: Option<f64>) -> f64 {
    if self.x.len() < 3 {
      panic!("A path hossza legalább 3 kell, hogy legyen.");
    }

    let p = p.unwrap_or(1.0);
    let sum1: f64 = (1..self.x.len())
      .map(|i| (self.x[i] - self.x[i - 1]).abs().powf(p))
      .sum();
    let sum2: f64 = (2..self.x.len())
      .map(|i| (self.x[i] - self.x[i - 2]).abs().powf(p))
      .sum();

    let vp = |increments: f64, l: usize, x_len: usize| -> f64 {
      1.0 / (2.0 * (x_len - l) as f64) * increments
    };

    let v1 = vp(sum1, 1, self.x.len());
    let v2 = vp(sum2, 2, self.x.len());

    2.0 - (1.0 / p) * ((v2.ln() - v1.ln()) / LN_2)
  }

  /// Calculate the Higuchi fractal dimension of the path.
  pub fn higuchi_fd(&self, kmax: usize) -> f64 {
    let n_times = self.x.len();

    let mut lk = Array1::<f64>::zeros(kmax);
    let mut x_reg = Array1::<f64>::zeros(kmax);
    let mut y_reg = Array1::<f64>::zeros(kmax);

    for k in 1..=kmax {
      let mut lm = Array1::zeros(k);

      for m in 0..k {
        let mut ll = 0.0;
        let n_max = ((n_times - m - 1) as f64 / k as f64).floor() as usize;

        for j in 1..n_max {
          ll += (self.x[m + j * k] - self.x[m + (j - 1) * k]).abs();
        }

        ll /= k as f64;
        ll *= (n_times - 1) as f64 / (k * n_max) as f64;
        lm[m] = ll;
      }

      lk[k - 1] = lm.iter().sum::<f64>() / k as f64;
      x_reg[k - 1] = (1.0 / k as f64).ln();
      y_reg[k - 1] = lk[k - 1].ln();
    }

    let (slope, _) =
      linear_regression(x_reg.as_slice().unwrap(), y_reg.as_slice().unwrap()).unwrap();
    slope
  }
}

#[cfg(test)]
mod tests {
  use super::*;
  use approx::assert_relative_eq;
  use stochastic_rs::{process::fbm::Fbm, Sampling};

  #[test]
  fn test_variogram() {
    let hurst = 0.75;
    let x = Fbm::new(&Fbm {
      hurst,
      n: 10_000,
      t: None,
      m: None,
      ..Default::default()
    });
    let fd = FractalDim::new(x.sample());
    let result = fd.variogram(None);
    assert_relative_eq!(2.0 - result, hurst, epsilon = 1e-1);
  }

  #[test]
  fn test_higuchi_fd() {
    let hurst = 0.75;
    let x = Fbm::new(&Fbm {
      hurst,
      n: 10_000,
      t: None,
      m: None,
      ..Default::default()
    });
    let fd = FractalDim::new(x.sample());
    let result = fd.higuchi_fd(10);
    assert_relative_eq!(2.0 - result, hurst, epsilon = 1e-2);
  }
}