Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- extern crate rand;
- use self::rand::distributions::Uniform;
- use self::rand::Rng;
- fn repartition<T: Clone>(xs: &Vec<T>, rng: &mut impl Rng) -> (Vec<T>, Vec<T>) {
- let n = xs.len();
- let mut xs_t = Vec::with_capacity(n / 2);
- let mut xs_f = Vec::with_capacity(n / 2);
- let u = Uniform::new(0.0, 1.0);
- xs.iter().for_each(|x| {
- if rng.sample(u) < 0.5 {
- xs_t.push(x.clone())
- } else {
- xs_f.push(x.clone())
- }
- });
- (xs_t, xs_f)
- }
- fn perm_test<T, F>(
- xs: &Vec<T>,
- ys: &Vec<T>,
- func: F, // Could be very slow to execute
- n_perms: usize,
- mut rng: &mut impl Rng,
- ) -> f64
- where
- F: Fn(&Vec<T>, &Vec<T>) -> f64,
- T: Clone,
- {
- let f0 = func(&xs, &ys);
- let mut xy = xs.clone();
- xy.append(&mut ys.clone());
- let incr = 1.0 / n_perms as f64;
- (0..n_perms).fold(0.0, |acc, _| {
- let (x, y) = repartition(&xy, &mut rng);
- if func(&x, &y) > f0 {
- acc + incr
- } else {
- acc
- }
- })
- }
- fn mean_diff(xs: &Vec<f64>, ys: &Vec<f64>) -> f64 {
- let mean_xs: f64 = xs.iter().fold(0.0, |acc, x| acc + x) / xs.len() as f64;
- let mean_ys: f64 = ys.iter().fold(0.0, |acc, y| acc + y) / ys.len() as f64;
- (mean_xs - mean_ys).abs()
- }
- fn main() {
- let mut rng = rand::thread_rng();
- let ux = Uniform::new(0.0, 1.0);
- let uy = Uniform::new(0.2, 1.2);
- let xs: Vec<f64> = (0..100).map(|_| rng.sample(&ux)).collect();
- let ys: Vec<f64> = (0..100).map(|_| rng.sample(&uy)).collect();
- println!("X-Y mean diff: {}", mean_diff(&xs, &ys));
- let xy_p_val = perm_test(&xs, &ys, mean_diff, 1_000, &mut rng);
- let xx_p_val = perm_test(&xs, &ys, mean_diff, 1_000, &mut rng);
- println!("X-Y permutation test: {}", xy_p_val);
- println!("X-X permutation test: {}", xx_p_val);
- }
Add Comment
Please, Sign In to add comment