Guest User

Untitled

a guest
Jun 22nd, 2018
105
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.81 KB | None | 0 0
  1. extern crate rand;
  2.  
  3. use self::rand::distributions::Uniform;
  4. use self::rand::Rng;
  5.  
  6. fn repartition<T: Clone>(xs: &Vec<T>, rng: &mut impl Rng) -> (Vec<T>, Vec<T>) {
  7. let n = xs.len();
  8. let mut xs_t = Vec::with_capacity(n / 2);
  9. let mut xs_f = Vec::with_capacity(n / 2);
  10.  
  11. let u = Uniform::new(0.0, 1.0);
  12. xs.iter().for_each(|x| {
  13. if rng.sample(u) < 0.5 {
  14. xs_t.push(x.clone())
  15. } else {
  16. xs_f.push(x.clone())
  17. }
  18. });
  19. (xs_t, xs_f)
  20. }
  21.  
  22. fn perm_test<T, F>(
  23. xs: &Vec<T>,
  24. ys: &Vec<T>,
  25. func: F, // Could be very slow to execute
  26. n_perms: usize,
  27. mut rng: &mut impl Rng,
  28. ) -> f64
  29. where
  30. F: Fn(&Vec<T>, &Vec<T>) -> f64,
  31. T: Clone,
  32. {
  33. let f0 = func(&xs, &ys);
  34.  
  35. let mut xy = xs.clone();
  36. xy.append(&mut ys.clone());
  37.  
  38. let incr = 1.0 / n_perms as f64;
  39. (0..n_perms).fold(0.0, |acc, _| {
  40. let (x, y) = repartition(&xy, &mut rng);
  41. if func(&x, &y) > f0 {
  42. acc + incr
  43. } else {
  44. acc
  45. }
  46. })
  47. }
  48.  
  49. fn mean_diff(xs: &Vec<f64>, ys: &Vec<f64>) -> f64 {
  50. let mean_xs: f64 = xs.iter().fold(0.0, |acc, x| acc + x) / xs.len() as f64;
  51. let mean_ys: f64 = ys.iter().fold(0.0, |acc, y| acc + y) / ys.len() as f64;
  52. (mean_xs - mean_ys).abs()
  53. }
  54.  
  55.  
  56. fn main() {
  57. let mut rng = rand::thread_rng();
  58. let ux = Uniform::new(0.0, 1.0);
  59. let uy = Uniform::new(0.2, 1.2);
  60. let xs: Vec<f64> = (0..100).map(|_| rng.sample(&ux)).collect();
  61. let ys: Vec<f64> = (0..100).map(|_| rng.sample(&uy)).collect();
  62.  
  63. println!("X-Y mean diff: {}", mean_diff(&xs, &ys));
  64.  
  65. let xy_p_val = perm_test(&xs, &ys, mean_diff, 1_000, &mut rng);
  66. let xx_p_val = perm_test(&xs, &ys, mean_diff, 1_000, &mut rng);
  67.  
  68. println!("X-Y permutation test: {}", xy_p_val);
  69. println!("X-X permutation test: {}", xx_p_val);
  70. }
Add Comment
Please, Sign In to add comment