Random — 乱数・確率分布

概要

sangi の乱数モジュールは 3 つの層からなる。

乱数エンジン — 再現可能なシード管理付き RandomEngine と、高速 PRNG (xoshiro256++) / 正規乱数生成器 (Ziggurat 等)。
確率分布 — 30 種超の連続・離散分布。各分布は確率密度 (pdf) / 累積分布 (cdf) / 分位点 (quantile) / モーメント (mean, variance) とサンプリング (sample) を統一インターフェースで提供する。
準乱数 — Halton / Sobol の低食い違い列。モンテカルロ積分の収束を $O(1/\sqrt{N})$ から準 $O(1/N)$ に改善する。

係数型 $T$ は double / float のほか多倍長 Float でも使える。サンプリングは標準 std::mt19937& を受け取る。

ビルド

#include <math/random/random.hpp>          // RandomEngine, Uniform/Normal/MultivariateNormal
#include <math/random/distributions.hpp>   // 連続・離散分布 (Gamma, Beta, Poisson, ...)
#include <math/random/quasi_random.hpp>    // Halton / Sobol 列

ヘッダオンリー。リンクするライブラリは不要 (ただし random.hpp は線形代数を使うため sangi_rational 系の static ライブラリと一緒にリンクすると確実)。

乱数エンジン

型	説明
`RandomEngine`	`std::mt19937` ラッパー。`uniform()`・`uniform(a,b)`・`normal(μ,σ)`・`uniformInt(a,b)`・`shuffle(c)`・`choice(c)` と再現可能シード
`random::Xoshiro256pp`	256-bit 高速 PRNG (xoshiro256++, 周期 $2^{256}{-}1$)
`random::NormalGenerator<T>`	Ziggurat 法の高速正規乱数生成器 (95% 超がテーブル参照のみ)
`random::MarsagliaPolarGenerator<T>`	Marsaglia polar 法 (sin/cos 不要)
`random::BoxMullerGenerator<T>`	Box-Muller 法 (参照・検証用)

RandomEngine eng(42);              // シード 42
double u = eng.uniform();          // [0,1)
double x = eng.uniform(2.0, 5.0);  // [2,5)
int d = eng.uniformInt(1, 6);      // サイコロ {1..6}
// 同じシードからは同じ列が再現される (RandomEngine eng2(42); eng2.uniform() == u)

分布の共通インターフェース

各分布クラスはパラメータをコンストラクタに取り、次の共通メンバを持つ (離散分布は pdf の代わりに確率質量 pmf(int k))。

メンバ	説明
`pdf(x)` / `pmf(k)`	確率密度 / 確率質量
`cdf(x)`	累積分布関数
`quantile(p)`	分位点 (逆 CDF)。一部の分布で提供
`mean()` / `variance()`	平均 / 分散
`sample(std::mt19937& gen)`	1 標本を生成
`sample(std::mt19937& gen, size_t n)`	$n$ 標本を `std::vector` で生成

NormalDistribution<double> nd(0.0, 1.0);   // 標準正規 N(0,1)
nd.mean();            // 実行結果: 0
nd.variance();        // 実行結果: 1
nd.pdf(0.0);          // 実行結果: 0.398942280401433   (= 1/sqrt(2π))
nd.cdf(0.0);          // 実行結果: 0.5
nd.quantile(0.975);   // 実行結果: 1.95996398612019    (両側 95% 点)

連続分布

クラス	パラメータ	説明
`UniformDistribution<T>`	`a, b`	一様分布 $U(a,b)$
`NormalDistribution<T>`	`μ, σ`	正規分布 $N(\mu,\sigma^2)$ (Ziggurat サンプリング)
`ExponentialDistribution<T>`	`λ`	指数分布 (率 $\lambda$)
`GammaDistribution<T>`	`α, β`	ガンマ分布 (Marsaglia-Tsang)
`BetaDistribution<T>`	`α, β`	ベータ分布 ($[0,1]$ 上)
`ChiSquaredDistribution<T>`	`k`	カイ二乗分布 $\chi^2(k)$
`FDistribution<T>`	`d1, d2`	F 分布
`StudentTDistribution<T>`	`ν`	スチューデント $t$ 分布
`CauchyDistribution<T>`	`x0, γ`	コーシー分布 (平均・分散は未定義)
`RayleighDistribution<T>`	`σ`	レイリー分布
`WeibullDistribution<T>`	`k, λ`	ワイブル分布
`LognormalDistribution<T>`	`μ, σ`	対数正規分布
`LaplaceDistribution<T>`	`μ, b`	ラプラス分布 (両側指数)
`ParetoDistribution<T>`	`α, xm`	パレート分布 (べき則)
`LogisticDistribution<T>`	`μ, s`	ロジスティック分布
`GumbelDistribution<T>`	`μ, β`	ガンベル分布 (極値)

パラメータ (主なもの):

記号	型	定義域・意味
`a, b`	`T`	区間の下限・上限 ($a < b$)
`μ, σ`	`T`	平均・標準偏差 ($\sigma > 0$)
`λ`	`T`	率パラメータ ($\lambda > 0$)
`α, β`	`T`	形状・率/尺度パラメータ ($> 0$)
`k` / `ν`	`T`	自由度 ($> 0$)

UniformDistribution<double> ud(0.0, 1.0);
ud.mean();        // 実行結果: 0.5
ud.variance();    // 実行結果: 0.0833333333333333   (= 1/12)
ud.pdf(0.5);      // 実行結果: 1

ExponentialDistribution<double> ed(2.0);  // 率 λ = 2
ed.mean();        // 実行結果: 0.5
ed.variance();    // 実行結果: 0.25

離散分布

クラス	パラメータ	説明
`BernoulliDistribution<T>`	`p`	ベルヌーイ分布 (1 試行 {0,1})
`PoissonDistribution<T>`	`λ`	ポアソン分布 (計数)
`BinomialDistribution<T>`	`n, p`	二項分布 ($n$ 試行の成功数)
`GeometricDistribution<T>`	`p`	幾何分布 (初成功までの失敗数)
`NegativeBinomialDistribution<T>`	`r, p`	負の二項分布 ($r$ 回成功までの失敗数)

PoissonDistribution<double> pd(3.0);
pd.mean();        // 実行結果: 3
pd.variance();    // 実行結果: 3

BinomialDistribution<double> bd(10, 0.5);
bd.mean();        // 実行結果: 5
bd.variance();    // 実行結果: 2.5
bd.pmf(5);        // 実行結果: 0.24609375   (= C(10,5)/2^10)

多変量・特殊分布

クラス	パラメータ	説明
`MultivariateNormal<T>`	`mean, covariance`	多変量正規分布 (Cholesky 前計算で `pdf`/`sample`)
`DirichletDistribution<T>`	`alpha[]`	ディリクレ分布 ($K$ 次元単体上)
`ArcsineDistribution<T>`	`a, b`	逆正弦分布 (U 字)
`InverseGammaDistribution<T>`	`α, β`	逆ガンマ分布
`InverseGaussianDistribution<T>`	`μ, λ`	逆ガウス (Wald) 分布
`SkewNormalDistribution<T>`	`ξ, ω, α`	歪正規分布
`TriangularDistribution<T>`	`a, b, c`	三角分布 (最頻値 $c$)
`NoncentralChiSquaredDistribution<T>`	`k, λ`	非心カイ二乗分布
`NoncentralTDistribution<T>`	`ν, δ`	非心 $t$ 分布
`NoncentralFDistribution<T>`	`d1, d2, λ`	非心 F 分布

準乱数 (低食い違い列)

準乱数列は一様性が高く、モンテカルロ積分の収束を擬似乱数の $O(1/\sqrt{N})$ から準 $O((\log N)^D / N)$ に改善する。

クラス	説明
`HaltonSequence<Dim>`	Halton 列 (素数基底の van der Corput, $1 \le \text{Dim} \le 50$)
`SobolSequence<Dim>`	Sobol 列 (Gray code, Joe-Kuo 方向数, $1 \le \text{Dim} \le 8$)

メンバ: next() (次の点 std::array<double, Dim>)、generate(n) ($n$ 点)、 index()、reset()。

HaltonSequence<2> h;     // 基底 (2, 3)
auto p0 = h.next();      // 実行結果: (0, 0)
auto p1 = h.next();      // 実行結果: (0.5, 0.333333333333333)
auto p2 = h.next();      // 実行結果: (0.25, 0.666666666666667)
auto p3 = h.next();      // 実行結果: (0.75, 0.111111111111111)

使用例

#include <math/random/random.hpp>
#include <math/random/distributions.hpp>
#include <random>
#include <iostream>
using namespace sangi;

int main() {
    // 閉形式の値 (決定的)
    NormalDistribution<double> nd(0.0, 1.0);
    std::cout << nd.pdf(0.0) << "\n";       // 0.398942280401433
    std::cout << nd.quantile(0.975) << "\n"; // 1.95996398612019

    // モンテカルロ: N(2, 0.5) から 10 万標本の標本平均 ≈ 2.0
    std::mt19937 gen(12345);
    NormalDistribution<double> g(2.0, 0.5);
    double s = 0.0;
    for (int i = 0; i < 100000; ++i) s += g.sample(gen);
    std::cout << s / 100000 << "\n";        // ≈ 2.0 (シード依存)
}