[Rust] 正規分布に従う乱数を生成する [Box-Mullerの実装]

正規分布を実装したので, 正規乱数も作りたいなぁと思い実装しました. これ使ってガウスノイズ付きのデータ生成 & 線形回帰とかやっていこうかと思います.

そんなことより, ボックス=ミュラー法, 今日知りました.

うん, 普通に感動したわ.(小並感)

理論の参考はwikipediaと,

Rustで実装する.

実装は脳死で数式を打ち込めばオッケーですね.
とりあえず正規乱数を生成する関数の部分だけ. nは個数指定のためです. 一様分布からの乱数はrand使ってんのかい！ってツッコミはなしで.

    // ボックスミュラー法で標準正規分布に従った乱数を生成する.
    fn norm_random(&mut self, n: i32) -> Vec<f64> {
        let mut y = vec![];
        let mut rng = rand::thread_rng();
        for i in 0..n {
            let u1: f64 = rng.gen_range(0f64,1f64);
            let u2: f64 = rng.gen_range(0f64,1f64);
            let a = (- 2f64 * u1.ln()).sqrt();
            let b = (2f64*f64::consts::PI*u2).cos();
            y.push( self.mu + self.sigma*a*b )
        }
        y
    }

実際に100個の乱数生成して, ヒストグラムにしてみました.

f:id:ket-30:20170202043150p:plain:w400:h400

histgramはgnuplotのboxesで書いてるんですけど, そこに渡すベクトルの作り方がくそ汚い. しかもx軸にはboxの中心になる数を与えるのにミスってるし. スマートな方法わかるまで, とりあえずこれで. (match初めてだったし, 多少はね…)

以下全コード.

fn main() {

    //平均 0, 分散 1のNormalDistribution型が生成される.
    let mut norm = NormalDistributionBilder::new().set();

    // 標準正規分布に従う乱数の生成.
    let y: Vec<f64> = norm.norm_random(100);

    // histgram生成のための下準備
    // クソコード
    let mut n: Vec<i32> = vec![0; 6];
    let m: Vec<i32> = vec![−3, -2, -1, 0, 1, 2];

    for e in &y {
        match () {
            _ if (e < &-2f64) => n[0] += 1,
            _ if (&-2f64 <= e && e < &-1f64) => n[1] += 1,
            _ if (&-1f64 <= e && e < &0f64) => n[2] += 1,
            _ if ( &0f64 <= e && e < &1f64) => n[3] += 1,
            _ if ( &1f64 <= e && e < &2f64) => n[4] += 1,
            _ if ( &2f64 <= e) => n[5] += 1,
            _ => println!("error"),
        };
    }

    //プロット
    let mut fg = Figure::new();
    fg.axes2d().boxes(&m, &n, &[Color("blue")]);
    fg.set_terminal("png", "norm_random_hist.png");
    fg.show();
}

struct NormalDistribution {
    mu : f64,
    sigma : f64,
}

impl NormalDistribution {
    fn new(mu: f64, sigma: f64) -> NormalDistribution {
        NormalDistribution{mu: mu, sigma: sigma,}
    }

    fn calc(&mut self, x: &Vec<f64>) -> Vec<f64> {
        let mut y: Vec<f64> = vec![];
        for e in x{
            let a = - (e - self.mu).powi(2) / (2.0_f64 * self.sigma.powi(2));
            let b = (2.0_f64 * f64::consts::PI).sqrt() * self.sigma;
            y.push( a.exp() / b );
        }

        y
    }

    fn show_info(&mut self){
        println!("------- Information -------");
        println!("Average: {}, Variance: {}", self.mu, self.sigma);
        println!("---------------------------");
    }

    // ボックスミュラー法で標準正規分布に従った乱数を生成する.
    fn norm_random(&mut self, n: i32) -> Vec<f64> {
        let mut y = vec![];
        let mut rng = rand::thread_rng();
        for i in 0..n {
            let u1: f64 = rng.gen_range(0f64,1f64);
            let u2: f64 = rng.gen_range(0f64,1f64);
            let a = (- 2f64 * u1.ln()).sqrt();
            let b = (2f64*f64::consts::PI*u2).cos();
            y.push( self.mu + self.sigma*a*b )
        }
        y
    }

    fn plot(&mut self, x: &Vec<f64>) {
        let mut fg = Figure::new();
        let y: Vec<f64> = self.calc(&x);
        fg.axes2d()
        .lines(x, &y, &[Color("blue")]); 
        fg.set_terminal("png", "NormalDistribution.png");
        fg.show();
    }
}

struct NormalDistributionBilder {
    mu: f64,
    sigma: f64,
}

impl NormalDistributionBilder {
    fn new() -> NormalDistributionBilder {
        NormalDistributionBilder{mu: 0f64, sigma: 1f64,}
    }

    fn mu(&mut self, mu: f64) -> &mut NormalDistributionBilder {
        self.mu = mu;
        self
    }

    fn sigma(&mut self, sigma: f64) -> &mut NormalDistributionBilder {
        self.sigma = sigma;
        self
    }

    fn set(&mut self) -> NormalDistribution {
        NormalDistribution{ mu: self.mu, sigma: self.sigma}
    }
}