またまた、Rustの統計分布の関数を使ってみる、続編。

去年12月5日のtest04で、Beta分布の分布関数 (CDF) を、"statrs"に加えて"probability"ライブラリを用いて算出してみた。いずれもいい感じの数値を出していた。

（ていうか越年してたんかい。。）

今回は、Beta分布の分布関数の逆関数 (Inverse of CDF; InvCDF) を使ってみる。

環境はWindows 10。

自分的プロジェクト名はtest05。

ソースなどはこちら (Github)

Rustソース

Rustのソースファイルで、Beta分布の分布関数の逆関数 (InvCDF) を出し、さらに得られた値を分布関数に入れて、InvCDFとCDF(InvCDF)の絶対誤差も出す、そんな関数を書く。

この関数を、statrs版とprobability版に分けて書く。この関数の部分を示すと以下のとおり。beta_statrs_crate()はstatrs版、beta_prob_crate()はprobability版。

main.rs

fn beta_statrs_crate( alpha: f64, beta: f64, p: f64) -> ( f64, f64, f64) 
{

    use statrs::distribution::Beta;
    use statrs::distribution::ContinuousCDF;

    let dist = Beta::new( alpha, beta).unwrap();
    let z_hat = dist.inverse_cdf( p);
    let p_hat = dist.cdf( z_hat);
    let abserr = ( p - p_hat).abs();

    return ( z_hat, p_hat, abserr);

}

fn beta_prob_crate( alpha: f64, beta: f64, p: f64) -> ( f64, f64, f64)
{

    use probability::distribution::Beta;
    use probability::distribution::Inverse;
    use probability::distribution::Distribution;

    let dist = Beta::new( alpha, beta, 0.0, 1.0);
    let z_hat = dist.inverse( p);
    let p_hat = dist.distribution( z_hat);
    let abserr = ( p - p_hat).abs();

    return ( z_hat, p_hat, abserr);

}

Cargo.tomlファイル

Cargo.tomlの[dependencies]の部分はtest04と同様に次のような感じにする。

Cargo.toml

[dependencies]
statrs = "0.16"
probability = "0.20.1"

結果

実行結果として出力されたものの一部は次のような感じ。

Test for InvCDF() of Beta(alpha,beta)
 Assume p = CDF(z)
 First, set p, and then calculate z_hat using InvCDF()
 After that, calculate p_hat = CDF(z_hat) for validation

Test two crates; "1" means "statrs" crate, "2" means "probability" crate

alpha   beta    p   z_hat_1 p_hat_1 abserr_1    z_hat_2 p_hat_2 abserr_2
0.5 0.5 0.01    0.000274658203125   0.010551064856976648    0.000551064856976648    0.00024671981713422173  0.010000000000000004    0.000000000000000003469446951953614
0.5 0.5 0.02    0.001007080078125   0.020206218608167428    0.0002062186081674272   0.0009866357858642186   0.019999999999999997    0.000000000000000003469446951953614
0.5 0.5 0.03    0.002227783203125   0.030059238187395775    0.00005923818739577563  0.002219017698459991    0.029999999999999985    0.000000000000000013877787807814457
0.5 0.5 0.04    0.003936767578125   0.03997011295735626 0.00002988704264374198  0.003942649342761082    0.03999999999999999 0.000000000000000013877787807814457
0.5 0.5 0.05    0.006134033203125   0.04991121962081565 0.00008878037918435416  0.006155829702431135    0.04999999999999998 0.000000000000000020816681711721685

表示を簡略化して整形した表で示すとこんな感じ。

alpha	beta	p	z_hat_1	p_hat_1	abserr_1	z_hat_2	p_hat_2	abserr_2
0.5	0.5	0.01	0.000274658	0.010551065	0.000551065	0.00024672	0.01	3.46945E-18
0.5	0.5	0.02	0.00100708	0.020206219	0.000206219	0.000986636	0.02	3.46945E-18
0.5	0.5	0.03	0.002227783	0.030059238	5.92382E-05	0.002219018	0.03	1.38778E-17
0.5	0.5	0.04	0.003936768	0.039970113	2.9887E-05	0.003942649	0.04	1.38778E-17
0.5	0.5	0.05	0.006134033	0.04991122	8.87804E-05	0.00615583	0.05	2.08167E-17

ほう。。

いいんだ、"abserr_2"はとても小さい。つまり、"probability"版のInvCDFはCDFの逆関数としていい精度を出している。

他方で、"abserr_1"はまあまあ大きい。。つまり"statrs"版のInvCDFはCDFの逆関数としてちょっと誤差がある。。

確かにalpha < 1とかbeta < 1の範囲では分布が変則的な形状になるので、逆関数を数値的に出すなら困難がありそうだけど、alpha >= 1 and beta >= 1の範囲でも"statrs"版ではズレが大きいように見えたぜ。。

次回、もうちょっと詳しく調べようか。。

kojiynet on Hatena

Rustの統計分布の関数を使う～逆関数の精度は要チェケラ～

Rustソース

Cargo.tomlファイル

結果