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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
// use std::f64;
//
// use ordered_float::NotNan;
// use bio::stats::{LogProb, Prob};
// use bio::stats::combinatorics::combinations;
//
// use model::{Variant, DiscreteAlleleFreqs, AlleleFreq};
// use priors::{Model, InfiniteSitesNeutralVariationModel, PairModel};
/*
pub struct NormalNormalModel {
inner: InfiniteSitesNeutralVariationModel,
ploidy: u32,
allele_freqs: DiscreteAlleleFreqs
}
impl NormalNormalModel {
pub fn new(ploidy: u32, heterozygosity: Prob) -> Self {
NormalNormalModel {
inner: InfiniteSitesNeutralVariationModel::new(2, ploidy, heterozygosity),
ploidy: ploidy,
allele_freqs: allele_freqs(ploidy)
}
}
}
impl PairModel<DiscreteAlleleFreqs, DiscreteAlleleFreqs> for NormalNormalModel {
fn prior_prob(&self, af_first: AlleleFreq, af_second: AlleleFreq, _: &Variant) -> LogProb {
let sample_prior = |af: AlleleFreq| {
let m = *af * self.ploidy as f64;
let valid_genotypes = combinations(self.ploidy, m);
let all_genotypes = combinations()
};
}
fn joint_prob<L, O>(
&self,
af_first: &DiscreteAlleleFreqs,
af_second: &DiscreteAlleleFreqs,
likelihood_first: &L,
likelihood_second: &O,
_: &Variant
) -> LogProb where
L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
O: Fn(AlleleFreq, AlleleFreq) -> LogProb
{
let p_second = LogProb::ln_sum_exp(&af_second.iter().map(|&af_second| {
likelihood_second(af_second, AlleleFreq(0.0))
}).collect_vec());
let prob = LogProb::ln_sum_exp(&af_first.iter().map(|&af_first| {
let p_first = likelihood_first(af_first, AlleleFreq(0.0));
let prob = p_first + p_second;
prob
}).collect_vec());
prob
}
fn marginal_prob<L, O>(
&self,
likelihood_first: &L,
likelihood_second: &O,
variant: &Variant
) -> LogProb where
L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
O: Fn(AlleleFreq, AlleleFreq) -> LogProb
{
let p = self.joint_prob(
self.allele_freqs().0,
self.allele_freqs().1,
likelihood_first,
likelihood_second,
variant
);
p
}
fn map<L, O>(
&self,
likelihood_first: &L,
likelihood_second: &O,
_: &Variant
) -> (AlleleFreq, AlleleFreq) where
L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
O: Fn(AlleleFreq, AlleleFreq) -> LogProb
{
fn calc_map<L: Fn(AlleleFreq, AlleleFreq) -> LogProb>(likelihood: &L, afs: &DiscreteAlleleFreqs) -> AlleleFreq {
let (_, map) = afs.iter().minmax_by_key(|&af| {
let p = likelihood(*af, AlleleFreq(0.0));
NotNan::new(*p).expect("probability is NaN")
}).into_option().expect("prior has empty allele frequency spectrum");
*map
}
let map_first = calc_map(likelihood_first, self.allele_freqs().0);
let map_second = calc_map(likelihood_second, self.allele_freqs().1);
(map_first, map_second)
}
fn allele_freqs(&self) -> (&DiscreteAlleleFreqs, &DiscreteAlleleFreqs) {
(&self.allele_freqs, &self.allele_freqs)
}
}*/