initial commit
This commit is contained in:
7 files changed
+431
No files matched your search
@@ -0,0 +1 @@
|
||||
/target
|
||||
Generated
+214
@@ -0,0 +1,214 @@
|
||||
# This file is automatically @generated by Cargo.
|
||||
# It is not intended for manual editing.
|
||||
version = 3
|
||||
|
||||
[[package]]
|
||||
name = "anyhow"
|
||||
version = "1.0.52"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "84450d0b4a8bd1ba4144ce8ce718fbc5d071358b1e5384bace6536b3d1f2d5b3"
|
||||
|
||||
[[package]]
|
||||
name = "autocfg"
|
||||
version = "1.0.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "cdb031dd78e28731d87d56cc8ffef4a8f36ca26c38fe2de700543e627f8a464a"
|
||||
|
||||
[[package]]
|
||||
name = "cfg-if"
|
||||
version = "1.0.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd"
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-channel"
|
||||
version = "0.5.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e54ea8bc3fb1ee042f5aace6e3c6e025d3874866da222930f70ce62aceba0bfa"
|
||||
dependencies = [
|
||||
"cfg-if",
|
||||
"crossbeam-utils",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-deque"
|
||||
version = "0.8.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "6455c0ca19f0d2fbf751b908d5c55c1f5cbc65e03c4225427254b46890bdde1e"
|
||||
dependencies = [
|
||||
"cfg-if",
|
||||
"crossbeam-epoch",
|
||||
"crossbeam-utils",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-epoch"
|
||||
version = "0.9.6"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "97242a70df9b89a65d0b6df3c4bf5b9ce03c5b7309019777fbde37e7537f8762"
|
||||
dependencies = [
|
||||
"cfg-if",
|
||||
"crossbeam-utils",
|
||||
"lazy_static",
|
||||
"memoffset",
|
||||
"scopeguard",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-utils"
|
||||
version = "0.8.6"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "cfcae03edb34f947e64acdb1c33ec169824e20657e9ecb61cef6c8c74dcb8120"
|
||||
dependencies = [
|
||||
"cfg-if",
|
||||
"lazy_static",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "either"
|
||||
version = "1.6.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e78d4f1cc4ae33bbfc157ed5d5a5ef3bc29227303d595861deb238fcec4e9457"
|
||||
|
||||
[[package]]
|
||||
name = "getrandom"
|
||||
version = "0.2.3"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "7fcd999463524c52659517fe2cea98493cfe485d10565e7b0fb07dbba7ad2753"
|
||||
dependencies = [
|
||||
"cfg-if",
|
||||
"libc",
|
||||
"wasi",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "hermit-abi"
|
||||
version = "0.1.19"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33"
|
||||
dependencies = [
|
||||
"libc",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "lazy_static"
|
||||
version = "1.4.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
|
||||
|
||||
[[package]]
|
||||
name = "libc"
|
||||
version = "0.2.112"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "1b03d17f364a3a042d5e5d46b053bbbf82c92c9430c592dd4c064dc6ee997125"
|
||||
|
||||
[[package]]
|
||||
name = "megapermute"
|
||||
version = "0.1.0"
|
||||
dependencies = [
|
||||
"anyhow",
|
||||
"rand",
|
||||
"rayon",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "memoffset"
|
||||
version = "0.6.5"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "5aa361d4faea93603064a027415f07bd8e1d5c88c9fbf68bf56a285428fd79ce"
|
||||
dependencies = [
|
||||
"autocfg",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "num_cpus"
|
||||
version = "1.13.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "19e64526ebdee182341572e50e9ad03965aa510cd94427a4549448f285e957a1"
|
||||
dependencies = [
|
||||
"hermit-abi",
|
||||
"libc",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "ppv-lite86"
|
||||
version = "0.2.16"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "eb9f9e6e233e5c4a35559a617bf40a4ec447db2e84c20b55a6f83167b7e57872"
|
||||
|
||||
[[package]]
|
||||
name = "rand"
|
||||
version = "0.8.4"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "2e7573632e6454cf6b99d7aac4ccca54be06da05aca2ef7423d22d27d4d4bcd8"
|
||||
dependencies = [
|
||||
"libc",
|
||||
"rand_chacha",
|
||||
"rand_core",
|
||||
"rand_hc",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rand_chacha"
|
||||
version = "0.3.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88"
|
||||
dependencies = [
|
||||
"ppv-lite86",
|
||||
"rand_core",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rand_core"
|
||||
version = "0.6.3"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d34f1408f55294453790c48b2f1ebbb1c5b4b7563eb1f418bcfcfdbb06ebb4e7"
|
||||
dependencies = [
|
||||
"getrandom",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rand_hc"
|
||||
version = "0.3.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d51e9f596de227fda2ea6c84607f5558e196eeaf43c986b724ba4fb8fdf497e7"
|
||||
dependencies = [
|
||||
"rand_core",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rayon"
|
||||
version = "1.5.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "c06aca804d41dbc8ba42dfd964f0d01334eceb64314b9ecf7c5fad5188a06d90"
|
||||
dependencies = [
|
||||
"autocfg",
|
||||
"crossbeam-deque",
|
||||
"either",
|
||||
"rayon-core",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rayon-core"
|
||||
version = "1.9.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d78120e2c850279833f1dd3582f730c4ab53ed95aeaaaa862a2a5c71b1656d8e"
|
||||
dependencies = [
|
||||
"crossbeam-channel",
|
||||
"crossbeam-deque",
|
||||
"crossbeam-utils",
|
||||
"lazy_static",
|
||||
"num_cpus",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "scopeguard"
|
||||
version = "1.1.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd"
|
||||
|
||||
[[package]]
|
||||
name = "wasi"
|
||||
version = "0.10.2+wasi-snapshot-preview1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "fd6fbd9a79829dd1ad0cc20627bf1ed606756a7f77edff7b66b7064f9cb327c6"
|
||||
+11
@@ -0,0 +1,11 @@
|
||||
[package]
|
||||
name = "megapermute"
|
||||
version = "0.1.0"
|
||||
edition = "2021"
|
||||
|
||||
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
|
||||
|
||||
[dependencies]
|
||||
anyhow = "1.0.52"
|
||||
rand = "0.8.4"
|
||||
rayon = "1.5.1"
|
||||
+21
@@ -0,0 +1,21 @@
|
||||
MIT License
|
||||
|
||||
Copyright (c) 2022 Stephen Waits <steve@waits.net>
|
||||
|
||||
Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
of this software and associated documentation files (the "Software"), to deal
|
||||
in the Software without restriction, including without limitation the rights
|
||||
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
copies of the Software, and to permit persons to whom the Software is
|
||||
furnished to do so, subject to the following conditions:
|
||||
|
||||
The above copyright notice and this permission notice shall be included in all
|
||||
copies or substantial portions of the Software.
|
||||
|
||||
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
||||
SOFTWARE.
|
||||
@@ -0,0 +1,15 @@
|
||||
# megapermute
|
||||
|
||||
A million permutation hypothesis test utility.
|
||||
|
||||
## Usage
|
||||
|
||||
`megapermute` requires two input files, `control.dat` and `treatment.dat`. Each is a text
|
||||
file containing a list of numbers, treated internally as `f64`.
|
||||
|
||||
- `control.dat` the control sample
|
||||
- `treatment.dat` the treatment sample
|
||||
|
||||
With these two files in the current directly, running `megapermute` computes
|
||||
the empircal difference in means and p-value for the null hypothesis (H_0 ==
|
||||
H_a).
|
||||
+166
@@ -0,0 +1,166 @@
|
||||
use std::{
|
||||
fs::File,
|
||||
io::{BufRead, BufReader},
|
||||
iter,
|
||||
};
|
||||
|
||||
use anyhow::Result;
|
||||
use rand::{prelude::SliceRandom, thread_rng};
|
||||
use rayon::prelude::*;
|
||||
|
||||
// the number of threads to pass to rayon (which smartly splits it into actual number of threads)
|
||||
const N_THREADS: usize = 1_000;
|
||||
|
||||
// the number of permutation tests to run per "thread" (above)
|
||||
const N_PERMUTATIONS_PER_THREAD: usize = 1_000;
|
||||
|
||||
fn load_f64s(filename: &str) -> Result<Vec<f64>> {
|
||||
let f = BufReader::new(File::open(filename)?);
|
||||
Ok(f.lines()
|
||||
.map(|l| {
|
||||
l.expect("input error: unable to read line")
|
||||
.parse::<f64>()
|
||||
.expect("parse error: unable to parse f64")
|
||||
})
|
||||
.collect())
|
||||
}
|
||||
|
||||
// accepts an iterator of f64's and computes mean using Welford's online algorithm
|
||||
fn mean<'a>(iter: impl Iterator<Item = &'a f64>) -> f64 {
|
||||
iter.enumerate()
|
||||
.fold(0.0, |mu, (i, x)| mu + ((x - mu) / (i + 1) as f64))
|
||||
}
|
||||
|
||||
#[derive(Clone, PartialEq, Eq)]
|
||||
enum Group {
|
||||
Control,
|
||||
Treatment,
|
||||
}
|
||||
|
||||
fn permute(control: &[f64], treatment: &[f64], mu_diff: f64) -> f64 {
|
||||
let count: f64 = (0..N_THREADS)
|
||||
.into_par_iter()
|
||||
.map(|_| {
|
||||
// create an rng so we can shuffle later
|
||||
let mut rng = thread_rng();
|
||||
|
||||
// create an index array of [Control, Treatment] which we'll shuffle repeatedly to make our
|
||||
// selections during each permutation
|
||||
let mut index: Vec<Group> = iter::repeat(Group::Control)
|
||||
.take(control.iter().len())
|
||||
.chain(iter::repeat(Group::Treatment).take(treatment.iter().len()))
|
||||
.collect();
|
||||
|
||||
// do the actual permutations
|
||||
(0..N_PERMUTATIONS_PER_THREAD)
|
||||
.filter(|_| {
|
||||
// shuffle our group selections
|
||||
index.shuffle(&mut rng);
|
||||
|
||||
// comput emean of this "control permutation"
|
||||
let permuted_mean_control = mean(
|
||||
control
|
||||
.iter()
|
||||
.chain(treatment.iter())
|
||||
.enumerate()
|
||||
.filter_map(|(i, x)| {
|
||||
if index[i] == Group::Control {
|
||||
Some(x)
|
||||
} else {
|
||||
None
|
||||
}
|
||||
}),
|
||||
);
|
||||
|
||||
// comput emean of this "treatment permutation"
|
||||
let permuted_mean_treatment = mean(
|
||||
control
|
||||
.iter()
|
||||
.chain(treatment.iter())
|
||||
.enumerate()
|
||||
.filter_map(|(i, x)| {
|
||||
if index[i] == Group::Treatment {
|
||||
Some(x)
|
||||
} else {
|
||||
None
|
||||
}
|
||||
}),
|
||||
);
|
||||
|
||||
// select this permutation if the diff in these permuted means is more than the
|
||||
// empirical diff of means
|
||||
(permuted_mean_treatment - permuted_mean_control) > mu_diff
|
||||
})
|
||||
.count() as f64
|
||||
})
|
||||
.sum();
|
||||
|
||||
// p-value is the ratio of permutations where delta(mean) exceeded empircal delta(mean) to
|
||||
// total permutations
|
||||
let p_value = count / (N_THREADS * N_PERMUTATIONS_PER_THREAD) as f64;
|
||||
|
||||
// adjust for left or right tail
|
||||
if mu_diff < 0.0 {
|
||||
1.0 - p_value
|
||||
} else {
|
||||
p_value
|
||||
}
|
||||
}
|
||||
|
||||
fn main() -> Result<()> {
|
||||
// read data
|
||||
let control = load_f64s("control.dat")?;
|
||||
let treatment = load_f64s("treatment.dat")?;
|
||||
|
||||
// compute empircal means
|
||||
let mean_control = mean(control.iter());
|
||||
let mean_treatment = mean(treatment.iter());
|
||||
println!(" mu_control = {}", mean_control);
|
||||
println!(" mu_treatment = {}", mean_treatment);
|
||||
println!(
|
||||
"(mu_treatment - mu_control) = {}",
|
||||
mean_treatment - mean_control
|
||||
);
|
||||
|
||||
// run permutation test to compute p-value
|
||||
let pvalue = permute(&control, &treatment, mean_treatment - mean_control);
|
||||
println!(" p-value = {}", pvalue);
|
||||
|
||||
Ok(())
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
// import names from main scope
|
||||
use super::*;
|
||||
|
||||
const MEAN_EPSILON: f64 = 0.000001;
|
||||
const PVALUE_EPSILON: f64 = 0.001;
|
||||
|
||||
// for float comparison
|
||||
fn approx_eq(a: f64, b: f64, epsilon: f64) -> bool {
|
||||
(a - b).abs() < epsilon
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_mouse_data() {
|
||||
// test data from Table 2.1 in "An Introduction to the Bootstrap" (book)
|
||||
let control = vec![52.0, 104.0, 146.0, 10.0, 51.0, 30.0, 40.0, 27.0, 46.0];
|
||||
let treatment = vec![94.0, 197.0, 16.0, 38.0, 99.0, 141.0, 23.0];
|
||||
|
||||
// compute empircal means
|
||||
let mean_control = mean(control.iter());
|
||||
assert!(approx_eq(mean_control, 56.22222222222222, MEAN_EPSILON));
|
||||
let mean_treatment = mean(treatment.iter());
|
||||
assert!(approx_eq(mean_treatment, 86.85714285714286, MEAN_EPSILON));
|
||||
assert!(approx_eq(
|
||||
mean_treatment - mean_control,
|
||||
30.63492063492064,
|
||||
MEAN_EPSILON
|
||||
));
|
||||
|
||||
// run permutation test to compute p-value
|
||||
let pvalue = permute(&control, &treatment, mean_treatment - mean_control);
|
||||
assert!(approx_eq(pvalue, 0.13896357, PVALUE_EPSILON));
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user