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

// Copyright 2016 Johannes Köster.
// Licensed under the MIT license (http://opensource.org/licenses/MIT)
// This file may not be copied, modified, or distributed
// except according to those terms.

//! Utilities for Bayesian statistics.

pub mod bayes_factors;
pub mod model;
pub use self::bayes_factors::BayesFactor;
pub use self::model::Model;

use itertools::Itertools;
use ordered_float::OrderedFloat;

use crate::stats::LogProb;

/// For each of the hypothesis tests given as posterior error probabilities
/// (PEPs, i.e. the posterior probability of the null hypothesis), estimate the FDR
/// for the case that all null hypotheses with at most this PEP are rejected.
/// FDR is calculated as presented by Müller, Parmigiani, and Rice,
/// "FDR and Bayesian Multiple Comparisons Rules" (July 2006).
/// Johns Hopkin's University, Dept. of Biostatistics Working Papers. Working Paper 115.
///
/// # Returns
///
/// A vector of expected FDRs in the same order as the given PEPs.
pub fn expected_fdr(peps: &[LogProb]) -> Vec<LogProb> {
    // sort indices
    let sorted_idx =
        (0..peps.len()).sorted_by(|&i, &j| OrderedFloat(*peps[i]).cmp(&OrderedFloat(*peps[j])));
    // estimate FDR
    let mut expected_fdr = vec![LogProb::ln_zero(); peps.len()];
    for (j, (expected_fp, i)) in LogProb::ln_cumsum_exp(sorted_idx.clone().map(|i| peps[i]))
        .zip(sorted_idx)
        .enumerate()
    {
        let fdr = LogProb(*expected_fp - ((j + 1) as f64).ln());
        expected_fdr[i] = if fdr <= LogProb::ln_one() {
            fdr
        } else {
            LogProb::ln_one()
        };
    }

    expected_fdr
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::stats::LogProb;

    #[test]
    fn test_expected_fdr() {
        let peps = [
            LogProb(0.1f64.ln()),
            LogProb::ln_zero(),
            LogProb(0.25f64.ln()),
        ];
        let fdrs = expected_fdr(&peps);
        println!("{:?}", fdrs);

        assert_relative_eq!(*fdrs[1], *LogProb::ln_zero());
        assert_relative_eq!(*fdrs[0], *LogProb(0.05f64.ln()));
        assert_relative_eq!(*fdrs[2], *LogProb((0.35 / 3.0f64).ln()), epsilon = 0.000001);
    }
}