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

use crate::{Distribution, DistributionError, RandomVariable};

/// Sample b from posterior p(b|a) with likelihood p(a|b) and prior p(b)

pub struct ImportanceSampler<T, D, PD>
where
    T: RandomVariable,
    D: Distribution<Value = T, Condition = ()>,
    PD: Distribution<Value = T, Condition = ()>,
{
    distribution: D,
    proposal: PD,
}

#[derive(thiserror::Error, Debug)]
pub enum ImportanceSamplingError {
    #[error("out of range")]
    OutOfRange,
    #[error("Unknown error")]
    Unknown,
}

impl<T, D, PD> ImportanceSampler<T, D, PD>
where
    T: RandomVariable,
    D: Distribution<Value = T, Condition = ()>,
    PD: Distribution<Value = T, Condition = ()>,
{
    pub fn new(distribution: D, proposal: PD) -> Result<Self, DistributionError> {
        Ok(Self {
            distribution,
            proposal,
        })
    }

    pub fn expectation(&self, f: impl Fn(&T) -> f64, x: &[T]) -> Result<f64, DistributionError> {
        let wi_fxi = x
            .iter()
            .map(|xi| -> Result<_, DistributionError> {
                let wi =
                    self.distribution.p_kernel(&xi, &())? / self.proposal.p_kernel(&xi, &())?;
                let fxi = f(xi);
                Ok((wi, fxi))
            })
            .collect::<Result<Vec<_>, _>>()?;
        let sum_w = wi_fxi.iter().map(|&(wi, _)| wi).sum::<f64>();
        let sum_w_fx = wi_fxi.iter().map(|&(wi, fxi)| wi * fxi).sum::<f64>();

        Ok(sum_w_fx / sum_w)
    }
}