feat: rework clayton copula

rust-dd · Jan 9, 2025 · c9c0abe · c9c0abe
1 parent f02d4af
commit c9c0abe
Show file tree

Hide file tree

Showing 4 changed files with 149 additions and 51 deletions.
diff --git a/Cargo.toml b/Cargo.toml
@@ -43,11 +43,11 @@ num-complex = { version = "0.4.6", features = ["rand"] }
 plotly = { version = "0.10.0", features = ["plotly_ndarray"] }
 polars = { version = "0.43.1", features = ["lazy"] }
 prettytable-rs = "0.10.0"
-# pyo3 = { version = "0.22.3", features = ["extension-module", "abi3-py38"] }
 quadrature = "0.1.2"
 rand = "0.8.5"
 rand_distr = "0.4.3"
 rayon = "1.10.0"
+roots = "0.0.8"
 sci-rs = "0.4.1"
 scilib = "1.0.0"
 statrs = "0.17.1"

diff --git a/pyproject.toml b/pyproject.toml
diff --git a/src/stats/copulas/bivariate.rs b/src/stats/copulas/bivariate.rs
@@ -1,9 +1,11 @@
+use core::f64;
 use std::{cmp::Ordering, error::Error};
 
 use ndarray::{stack, Array1, Array2, Axis};
 use ndarray_rand::RandomExt;
 use ndarray_stats::QuantileExt;
 use rand_distr::Uniform;
+use roots::{find_root_brent, SimpleConvergency};
 
 pub mod clayton;
 pub mod frank;
@@ -18,6 +20,8 @@ pub enum CopulaType {
   Independence,
 }
 
+const EPSILON: f64 = 1e-12;
+
 pub trait Bivariate {
   fn r#type(&self) -> CopulaType;
 
@@ -121,41 +125,65 @@ pub trait Bivariate {
     Ok(())
   }
 
-  fn probability_density(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>>;
-
-  fn pdf(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>> {
-    self.probability_density(X)
-  }
+  fn pdf(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>>;
 
-  fn log_probability_density(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>> {
-    let pdf = self.probability_density(X)?;
+  fn log_pdf(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>> {
+    let pdf = self.pdf(X)?;
     let log_pdf = pdf.mapv(|val| (val + 1e-32).ln());
     Ok(log_pdf)
   }
 
-  fn cumulative_distribution(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>>;
-
-  fn cdf(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>> {
-    self.cumulative_distribution(X)
-  }
+  fn cdf(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>>;
 
   fn percent_point(&self, y: &Array1<f64>, V: &Array1<f64>) -> Result<Array1<f64>, Box<dyn Error>> {
-    todo!()
+    let n = y.len();
+    let mut results = Array1::zeros(n);
+
+    for i in 0..n {
+      let y_i = y[i];
+      let v_i = V[i];
+
+      let f = |u| self.partial_derivative_scalar(u, v_i).unwrap() - y_i;
+      let mut convergency = SimpleConvergency {
+        eps: f64::EPSILON,
+        max_iter: 50,
+      };
+      let min = find_root_brent(f64::EPSILON, 1.0, f, &mut convergency);
+      results[i] = min.unwrap_or(f64::EPSILON);
+    }
+
+    Ok(results)
   }
 
   fn ppf(&self, y: &Array1<f64>, V: &Array1<f64>) -> Result<Array1<f64>, Box<dyn Error>> {
     self.percent_point(y, V)
   }
 
-  fn partial_derivative(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn Error>> {
-    todo!()
+  fn partial_derivative(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
+    let n = X.nrows();
+    let mut X_prime = X.clone();
+    let mut delta = Array1::zeros(n);
+    for i in 0..n {
+      delta[i] = if X[[i, 1]] > 0.5 { -0.0001 } else { 0.0001 };
+      X_prime[[i, 1]] = X[[i, 1]] + delta[i];
+    }
+
+    let f = self.cdf(X).unwrap();
+    let f_prime = self.cdf(&X_prime).unwrap();
+
+    let mut deriv = Array1::zeros(n);
+    for i in 0..n {
+      deriv[i] = (f_prime[i] - f[i]) / delta[i];
+    }
+
+    Ok(deriv)
   }
 
-  fn partial_derivative_scalar(
-    &self,
-    u: Array1<f64>,
-    v: Array1<f64>,
-  ) -> Result<f64, Box<dyn Error>> {
-    todo!()
+  fn partial_derivative_scalar(&self, U: f64, V: f64) -> Result<f64, Box<dyn Error>> {
+    self.check_fit()?;
+    let X = stack![Axis(1), Array1::from(vec![U]), Array1::from(vec![V])];
+    let out = self.partial_derivative(&X);
+
+    Ok(out?.get(0).unwrap().clone())
   }
 }
diff --git a/src/stats/copulas/bivariate/clayton.rs b/src/stats/copulas/bivariate/clayton.rs
@@ -55,33 +55,97 @@ impl Bivariate for Clayton {
   }
 
   fn generator(&self, t: &Array1<f64>) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
-    todo!()
+    self.check_fit()?;
+
+    let theta = self.theta.unwrap();
+    Ok((1.0 / theta) * (t.powf(-theta) - 1.0))
   }
 
-  fn probability_density(
-    &self,
-    X: &Array2<f64>,
-  ) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
-    todo!()
+  fn pdf(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
+    self.check_fit()?;
+
+    let U = X.column(0);
+    let V = X.column(1);
+
+    let theta = self.theta.unwrap();
+    let a = (theta + 1.0) * (&U * &V).powf(-theta - 1.0);
+    let b = U.powf(-theta) + V.powf(-theta) - 1.0;
+    let c = -(2.0 * theta + 1.0) / theta;
+    Ok(a * b.powf(c))
   }
 
-  fn cumulative_distribution(
-    &self,
-    X: &Array2<f64>,
-  ) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
-    todo!()
+  fn cdf(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
+    self.check_fit()?;
+
+    let U = X.column(0);
+    let V = X.column(1);
+
+    let V_all_zeros = V.iter().all(|&v| v == 0.0);
+    let U_all_zeros = U.iter().all(|&u| u == 0.0);
+
+    if V_all_zeros || U_all_zeros {
+      let shape = V.shape();
+      return Ok(Array1::zeros(shape[0]));
+    }
+
+    let theta = self.theta.unwrap();
+    let mut cdfs = Array1::<f64>::zeros(U.len());
+
+    for i in 0..U.len() {
+      let u = U[i];
+      let v = V[i];
+
+      if u > 0.0 && v > 0.0 {
+        cdfs[i] = (u.powf(-theta) + v.powf(-theta) - 1.0).powf(-1.0 / theta);
+      } else {
+        cdfs[i] = 0.0;
+      }
+    }
+
+    Ok(cdfs)
   }
 
   fn percent_point(
     &self,
     y: &Array1<f64>,
     V: &Array1<f64>,
   ) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
-    todo!()
+    self.check_fit()?;
+
+    let theta = self.theta.unwrap();
+
+    if theta == 0.0 {
+      return Ok(V.clone());
+    }
+
+    let a = y.powf(theta / (-1.0 - theta));
+    let b = V.powf(theta);
+
+    let b_all_zeros = b.iter().all(|&v| v == 0.0);
+
+    if b_all_zeros {
+      return Ok(Array1::ones(V.len()));
+    }
+
+    Ok(((a + &b - 1.0) / b).powf(-1.0 / theta))
   }
 
   fn partial_derivative(&self, X: &Array2<f64>) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
-    todo!()
+    self.check_fit()?;
+
+    let U = X.column(0);
+    let V = X.column(1);
+
+    let theta = self.theta.unwrap();
+    let A = V.powf(-theta - 1.0);
+
+    if A.is_all_infinite() {
+      return Ok(Array1::zeros(V.len()));
+    }
+
+    let B = V.powf(-theta) + U.powf(-theta) - 1.0;
+    let h = B.powf((-1.0 - theta) / theta);
+    Ok(A * h)
   }
 
   fn compute_theta(&self) -> f64 {
@@ -94,3 +158,25 @@ impl Bivariate for Clayton {
     2.0 * tau / (1.0 - tau)
   }
 }
+
+#[cfg(test)]
+mod tests {
+  use ndarray::array;
+
+  use crate::stats::copulas::bivariate::Bivariate;
+
+  use super::Clayton;
+
+  #[test]
+  fn test_fit() {
+    let data = array![[0.5, 0.5], [0.5, 0.4], [0.5, 0.3], [0.2, 0.2], [0.5, 0.1]];
+    let mut clayton = Clayton::new();
+    clayton.fit(&data).unwrap();
+
+    let data1 = array![0.5, 0.4, 0.3, 0.2, 0.1];
+    let data2 = array![0.5, 0.4, 0.3, 0.2, 0.1];
+
+    let percentile_value = clayton.percent_point(&data1, &data2).unwrap();
+    println!("{:?}", percentile_value);
+  }
+}