rspp 0.1.7

rust probolistic programming.
Documentation 
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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

use rand::distributions::Uniform;
extern crate nalgebra as na;
use na::*;
use nalgebra::{DMatrix, DVector};
use rand::prelude::Distribution;
use rand::seq::SliceRandom;
use rand::thread_rng;
use rand::Rng; // 0.8.5
use rand_distr::{Normal, NormalError};
use statrs::statistics::{MeanN, VarianceN};
use std::f64::consts::PI;
/*
先生成一个0到N间的随机整数i，代表选择第i列；

再生成一个0到1间的随机数，若其小于事件i占第i列矩形的面积的比例，则表示接受事件i，否则，接收第i列中不是事件i的另一个事件。

其实你可以算下这种方式事件i的概率，完全对应原来的概率分布。
*/

/// # Example
///
/// ```
/// use rspp::tools::resample;
/// let data = vec![1, 2, 3, 4, 5];
/// let weights = vec![0.4, 0.3, 0.2, 0.2, 0.6];
/// let ll = resample::resample_dup::<i32>(&data, &weights, 10);
/// dbg!("重采样===",ll);
/// ```
pub fn resample_dup<T: Clone>(data: &Vec<T>, weights: &Vec<f64>, size: usize) -> Vec<T> {
    let mut rng = rand::thread_rng();
    let n = data.len();
    let uni = Uniform::new(0.0, 1.0);
    let total_weight: f64 = weights.iter().sum();
    let norm_weight: Vec<f64> = weights.iter().map(|v| v / total_weight).collect();

    let mut res = Vec::new();
    while res.len() < size {
        let index = rng.gen_range(0..n);
        let threshold = uni.sample(&mut rng);
        if norm_weight[index] > threshold {
            let el = &data[index];
            res.push(el.clone());
        }
    }
    res
}
/// # Example
///
/// ```
/// use rspp::tools::resample::jacbin_one_row;
/// fn test_jacbin_f(x: &Vec<f64>) -> f64 {
///   x.iter().map(|i| i.powf(2.0)).sum()
/// }
/// let x: Vec<f64> = vec![1.0, 2.0];
/// let jac = jacbin_one_row(&x, test_jacbin_f);
/// dbg!("雅可比", jac);
/// ```
pub fn jacbin_one_row(x: &Vec<f64>, f: fn(&Vec<f64>) -> f64) -> Vec<f64> {
    let dt = 1e-8;
    let mut jac = Vec::new();
    for col in 0..x.len() {
        let djx = if x[col] != 0.0 { x[col] * dt.abs() } else { dt };
        let x1: Vec<f64> = x
            .iter()
            .enumerate()
            // 其他稳着不动,在col方向步进
            .map(|(i, v)| if i != col { *v } else { *v + djx })
            .collect();

        let res = (f(&x1) - f(x)) / djx;
        jac.push(res);
    }
    jac
}

pub fn random_from_stdnorm() -> f64 {
    let mut rng = rand::thread_rng();
    let normal = Normal::new(0.0, 1.0).unwrap();
    normal.sample(&mut rng)
}

pub fn gaussian_likehood(x: f64, std: f64) -> f64 {
    let n1 = -1. * x.powf(2.) / (2. * std.powf(2.));
    let n2 = (2. * PI).sqrt() * std;

    n1.exp() / n2
}

pub fn resample_cumsum(w: &Vec<f64>, np: usize) -> Vec<usize> {
    let wsum: f64 = w.iter().sum();
    let norm_w: Vec<f64> = w.iter().map(|v| v / wsum).collect();
    let w_cum: Vec<f64> = norm_w
        .into_iter()
        .scan(0.0, |acc, x| {
            *acc += x;
            Some(*acc)
        })
        .collect();
    let step = 1. / np as f64;
    let base: Vec<f64> = (0..np)
        .into_iter()
        .scan(0. - step, |acc, x| {
            *acc += step;
            Some(*acc)
        })
        .collect();

    let mut rng = rand::thread_rng();
    let uni = Uniform::new(0.0, step);
    let re_sample_id: Vec<f64> = base.iter().map(|v| v + uni.sample(&mut rng)).collect();
    let mut ind = 0;
    let mut index = Vec::new();
    for ip in 0..np {
        while re_sample_id[ip] > w_cum[ind] {
            ind += 1
        }
        index.push(ind);
    }
    index
}

pub fn rand_n_index(tlen: usize, n: usize) -> Vec<usize> {
    let l = 0..tlen;
    let mut numbers: Vec<usize> = l.into_iter().collect();
    numbers.shuffle(&mut rand::thread_rng());
    let res = numbers[0..n].to_vec();
    res
}
#[test]
fn test_resample() {
    let p = vec![1., 2., 3., 4.];
    resample_cumsum(&p, 100);
    let res = rand_n_index(6, 3);
    dbg!(res);
}