Skip to content

Latest commit

 

History

History
1129 lines (817 loc) · 43.2 KB

043-random.md

File metadata and controls

1129 lines (817 loc) · 43.2 KB
Raw

乱数

乱数はプログラミングにおいてよく使う。例えば6面ダイスをプログラムで実装するには、1, 2, 3, 4, 5, 6までのいずれかの目を出す。

$ ./dice
1
5
$ ./dice
3
5 1 6
$ ./dice
10
5 1 6 6 1 6 6 2 4 2

このプログラムdiceは標準入力から整数型の値nを取り、1, 2, 3, 4, 5, 6のいずれかをそれぞれ$\frac{1}{6}$の確率でn個出力する。

まずこのdiceプログラムを作ることを目標にC++の乱数アルゴリズムである<random>の使い方を学んでいく。

疑似乱数

コンピューターで使われる乱数のほとんどは疑似乱数と呼ばれる方法で生成されている。さまざまなアルゴリズムがあるが、とても簡単に理解できる疑似乱数のアルゴリズムに、線形合同法(Linear congruential generator)がある。

線形合同法ではいまの乱数を$X_n$、次の乱数を$X_{n+1}$とすると、$X_{n+1}$は以下のように求められる。

$$ X_{n+1} = (a \times X_{n} + c) \bmod m $$

たとえば$a = 3, c = 5, m = 2^{sizeof(std::uint32_t) \times 8}$の場合で、$X_0 = 0$のとき、

$$X_0 = 0$$

$$X_1 = 3 \times 0 + 5 \bmod 2^{32}-1 = 5$$

$$X_2 = 3 \times X_1 + 5 \bmod 2^{32}-1 = 20$$

$$X_3 = 3 \times X_2 + 5 \bmod 2^{32}-1 = 65$$

「これはぜんぜん乱数ではない。予測可能じゃないか」と考えるかもしれない。しかし中でどのように乱数が生成されているかわからなければ、外部からは乱数のように見える。これが擬似乱数の考え方だ。

乱数エンジン

乱数エンジンは生の乱数を生成するライブラリだ。クラスで実装されている。

乱数エンジンはメンバー関数min()で最小値を、メンバー関数max()で最大値を、operator()で最小値から最大値の間の乱数を返す。

template < typename Engine >
void f( Engine & e )
{
    // 最小値
    auto a = e.min() ;
    // 最大値
    auto b = e.max() ;
    // 乱数
    auto r1 = e() ;
    // 次の乱数
    auto r2 = e() ;
}

乱数エンジンのオブジェクトeはoperator ()を呼び出すたび、つまりe()をするたびに変更される。これは疑似乱数のための内部状態を更新するためだ。そのため、乱数エンジンはconstでは新しい乱数を作るのに使えない。

標準ライブラリはデフォルトの乱数エンジンとしてstd::default_random_engineを提供している。

以下のプログラムはデフォルトの乱数エンジンから乱数を10個出力する。

int main()
{
    // 乱数エンジン
    std::default_random_engine e ;
    for ( int i = 0 ; i != 10 ; ++i )
    {
        // 乱数を出力
        std::cout << e() << "\n"sv ;
    }
}

標準ライブラリの提供する乱数エンジンにはさまざまなものがあるが、本書ではもう1つ、メルセンヌツイスターというアルゴリズムを実装した乱数エンジンを紹介する。std::mt19937だ。

std::mt19937を使うには、st::default_random_engineを置き換えるだけでいい。

int main()
{
    std::mt19937 e ;
    for ( int i = 0 ; i != 10 ; ++i )
    {
        std::cout << e() << "\n"sv ;
    }
}

メルセンヌツイスターはとても優秀な乱数エンジンだ。乱数が必要な多くの場面では、メルセンヌツイスターを使っておけばまず問題はない。

では乱数エンジンを使って、生の乱数を標準入力で得た個数だけ出力するプログラムを書いてみよう。

int main()
{
    // 乱数エンジン
    std::mt19937 e ;

    // 標準入力からnを得る
    unsigned n {} ;
    std::cin >> n ;
    // n個出力
    for ( unsigned int i = 0 ; i != n ; ++i )
    {
        std::cout << e() << " "sv ;
    }
}

実行結果は以下のようになる。

$ dice
10
3499211612 581869302 3890346734 3586334585 545404204 4161255391 3922919429 949333985 2715962298 1323567403

乱数エンジンで生成されるのは生の乱数だ。これは通常、32bit符号なし整数とか64bit符号なし整数で表現できる全範囲の値として生成される。これは実際に必要な乱数とは値の範囲が違う。実際に必要な乱数とは、例えば6面ダイスの場合は、int型で1, 2, 3, 4, 5, 6のいずれかの値がそれぞれ$\frac{1}{6}$の確率で出てほしい。

乱数分布

乱数分布とは生の乱数を望みの範囲の乱数に加工するためのライブラリだ。クラスで実装されている。

乱数分布ライブラリにもさまざまなものがあるが、6面ダイスのプログラムを実装するのに使うのはstd::uniform_int_distribution<T>だ。

この乱数分布ライブラリは、Tにほしい乱数の整数型を指定する。コンストラクター引数を2つ取るので、1つ目の引数に最小値、2つ目の引数に最大値を指定する。

std::uniform_int_distribution<int> d(a, b) ;

この乱数分布クラスの変数dは、$a \leq r \leq b$までの範囲のint型の乱数rを作り出す。

6面ダイスを作るには、d(a, b)をd(1, 6)にすればよい。

std::uniform_int_distribution<int> d(1, 6) ;

乱数分布クラスのオブジェクトdを作ったならば、operator()に乱数エンジンのオブジェクトを引数に渡すことで乱数が作れる。乱数エンジンのオブジェクトをeとすると、d(e)だ。

template < typename Engine, typename Distribution >
void f( Engine & e, Distribution d)
{
    // 乱数
    auto r1 = d(e) ;
    // 次の乱数   
    auto r2 = d(e) ;
    // 次の乱数
    auto r3 = d(e) ; 
}

以上の知識を利用して、プログラムdiceを作ってみよう。

int main()
{
    // 乱数エンジン
    std::mt19937 e ;
    // 乱数分布
    std::uniform_int_distribution<int> d(1, 6) ;

    // 入力を処理
    unsigned n {} ;
    std::cin >> n ;
    for ( unsigned int i = 0 ; i != n ; ++i )
    {
        // 乱数出力
        std::cout << d(e) << " "sv ;
    }
}

さっそく実行してみよう。

$ ./dice
5
5 1 6 6 1
$ ./dice
10
5 1 6 6 1 6 6 2 4 2
$ ./dice
20
5 1 6 6 1 6 6 2 4 2 1 4 2 2 4 6 6 6 6 6

この場合、乱数値は正の整数しか生成しないので、型をintではなくunsigned intにすることもできる。

std::uniform_int_distribution<unsigned int> d( 1, 6 ) ;

ただし、乱数の結果の型をunsigned intにすると、生成した乱数を使うときに負数が出てくるような計算で問題になる。例えば6面ダイスを2回振り、1回目の出目から2回目の出目を引き算するコードを書いてみよう。

int main()
{
    std::mt19937 e ;
    std::uniform_int_distribution<unsigned int> d( 1, 6 ) ;
    auto a = d(e) ; // 1回目
    auto b = d(e) ; // 2回目

    auto result = a - b ; // 結果
    std::cout << result ;
}

もし2回目の出目の方が1回目の出目より大きかった場合、結果は負数になってしまうが、unsigned int型は負数を表現できない。

そのため、通常は符号付きの整数型を使った方が安全だ。

また、分布クラスのテンプレートパラメーターにはデフォルトテンプレート実引数が指定されているので、デフォルトでよければ省略することもできる。

// std::uniform_int_distribution<int>と同じ
std:uniform_int_distribution d( 1, 6 ) ;

ところで、上のコードは動くのだが、別のプログラムを実行しても毎回同じ出力になる。これでは実用的な6面ダイスプログラムとは言えない。プログラムの実行ごとに結果を買えたい場合、シードを設定する。

シード

線形合同法を思い出してみよう。線形合同法で次の乱数$X_{n+1}$を計算するには、いまの乱数$X_{n}$に対して$X_{n+1} = (a \times X_{n} + c) \bmod m$という計算をする。

線形合同法とは現在の乱数値を内部状態として持ち、そこに計算を加えることで次の乱数を作り出すのだ。

一般化すると、疑似乱数は内部状態$S_n$を持ち、計算を加える関数$f(x)$を適用することで、次の内部状態$S_{n+1}=f(S_n)$を作り出すのだ。単純な線形合同法の場合、内部状態がそのまま乱数の値になるが、複雑な疑似乱数アルゴリズムでは、内部状態から乱数を求めるのにさらに計算を加えるものもある。

乱数エンジンをデフォルト初期化すると、この内部状態もデフォルト初期化される。そのため、いままで使っていた乱数は、プログラムの実行ごとに同じ乱数列を作り出すのだ。

疑似乱数の内部状態の初期状態を設定するための値をシード(seed)という。シードを設定するには、std::seed_seqというクラスのオブジェクトを乱数エンジンのコンストラクターに渡す。

int main()
{
    std::seed_seq s{124} ;
    std::mt19937 e(s) ;
    std::cout << e() ;
}

std::seed_seq s({n})のnの値を変更し、異なるシード値が異なる生の乱数を生成しているのを確かめよう。

シード値は乱数エンジンのメンバー関数seed(s)でも渡すことができる。

int main()
{
    std::seed_seq s{123} ;
    // シード値を設定
    std::mt19937 e(s) ;
    // 乱数を生成
    // 内部状態が変わる
    auto r1 = e() ;
    // シード値を設定
    e.seed(s) ;
    // 乱数を生成
    auto r2 = e() ;
    // 同じシード値による乱数は同じ値になる
    // r1 == r2
}

r1 == r2になるのは、同じシード値を渡して内部状態を設定しているからだ。

std::seed_seqには複数の符号なし32bit整数を渡すことができる。= {n1, n2, n3,...}のように初期化することもできるし、イテレーターを使って(first, last)のように設定することもできる。

int main()
{
    std::seed_seq s1 = {1,2,3,4,5} ;

    std::vector<std::uint32_t> v = {1,2,3,4,5} ; 
    std::seed_seq s2( std::begin(v), std::end(v) ) ;
}

乱数エンジンをコピーすると、その内部状態もコピーされる。これを利用して、乱数を保存しておくこともできる。

int main()
{
    std::mt19937 e1 ;
    // 2回乱数を生成
    e1() ; e1() ;
    // コピー、内部状態もコピーされる
    std::mt19937 e2 = e1 ;
    // true
    bool b1 = e1() == e2() ;
    // true
    bool b2 = e1() == e2() ;
}

乱数エンジンe1, e2は同じ状態を持っているので、同じ回数乱数生成をすると、同じ乱数列が生成される。

予測不可能な乱数

シード値を設定すれば乱数エンジンに異なった乱数列を生成させることができる。しかし、シード値はどうやって生成すればいいのだろうか。シード値をデフォルト初期化した乱数エンジンで生成しても意味がない。なぜならば初期状態はプログラムの実行ごとに同じだからだ。

int main()
{
    std::mt19937 e1 ;
    std::seed_seq s = { e1(), e1(), e1() } ;
    // 意味がない
    std::mt19937 e2(s) ;
}

内部状態を更新するのではない、本当に予測不可能な乱数を生成するには、ハードウェアの支援が必要だ。例えば放射性同位体がいつ放射性崩壊を起こすかは予測不可能だ。したがって放射線量を計測するガイガーカウンターの値は予測不可能だ。コンピューターにガイガーカウンターが取り付けられていれば、その値を読むことによって予測不可能な値を得ることができる。ほかにもコンピューターにはさまざまな予測不可能な値を得る方法がある。std::random_deviceはそのような実装依存のコンピューターの支援を受け、予測不可能な乱数を生成する乱数エンジンだ。

int main()
{
    std::random_device rd ;
    for ( int i = 0 ; i != 10 ; ++i )
        std::cout << rd() << " "sv ;
}

std::random_deviceを使えば、std::seed_seqを予測不可能な値で初期化できる。

// 予測不可能な乱数エンジン
std::random_device rd ;
// シード値
std::seed_seq s = { rd(), rd(), rd() } ;
// シード値を指定して初期化
std::mt19937 e(s) ;

十分なシード値の量

std::seed_seqの初期化では32bit符号なし整数をいくつでも指定できる。

std::random_device rd ;
std::seed_seq s = { rd(), rd(), rd(), ...} ;

ではいくつの値を渡せばいいのだろうかということは、初期化する乱数エンジンの内部状態のサイズによって異なってくる。現在、C++標準規格には乱数エンジンを適切に初期化する簡単な方法がない。

1つの目安としては、乱数エンジンのオブジェクトサイズがある。

int main()
{
    std::cout << sizeof(std::mt19937) ;   
}

これを筆者の環境で実行すると、5000と出力された。つまり筆者の使っているC++実装のstd::mt19937のオブジェクトサイズは5000バイトだということだ。

std::random_deviceはunsigned int型の乱数を返す。筆者の環境ではsizeof(unsigned int) == 4になる。すると$5000 \div 4 = 1250$ となる。とすると安全のためには、std::seed_seqにはstd::random_deviceの乱数を1250個渡すべきだろう。

std::random_device rd ;
std::seed_seq s = { rd(), rd(), rd(), ... /*残り1247個のrd()*/ } ;
std::mt19937 e( s ) ;

筆者の環境ではsizeof(std::default_random_engine) == 8であった。すると2個でよいことになる。

std::random_device rd ;
std::seed_seq s = { rd(), rd() } ;
std::default_random_engine e(s) ;

C++標準規格にはいずれ、乱数エンジンを予測不可能なシード値で適切に初期化する簡単な方法が追加されるはずだ。

乱数分布ライブラリ

生の乱数は使いづらい。生の乱数というのは$n$ビットの整数値だ。それに対して、我々が使いたい実用的な乱数というのは以下のようなものだ。

  • コイントスの結果
  • 6面ダイスを振った結果
  • 当選確率1%のくじ引きの結果
  • 浮動小数点数0.0から1.0の範囲の値

コイントスの結果は表か裏かの2値になる。いま、$n$ビットの整数値rのすべてのビットが等しく乱数ビットであるならば、2値の乱数は単にr & 0b1で得られる。

では6面ダイスはどうか。6面ダイスは$1 \leq r \leq 6$までの6通りの状態を持つ乱数が必要だ。6通りの状態を表現するには、少なくとも3ビットの乱数が必要になる計算だ。しかし、3ビットの乱数は実際には多すぎる。3ビットの乱数で表現できるのは$2^3=8$通りの状態だ。したがって、r & 0x111というわけにはいかない。

巷には間違った乱数の分布方法として、$a \leq n \leq b$の範囲の乱数nを生の乱数rから得るために、以下のような計算式を用いる方法が蔓延している。

$$ n = r \bmod |b-a| + a $$

この間違った計算式を使うと、6面ダイスの乱数値nは生の乱数r($0 \leq r$)から以下のようにC++で計算できる。

template < typename Engine >
int dice( Engine & e )
{
    auto r = e() ;
    auto n = r % 6 + 1 ;
    return n ;
}

剰余を使って値を$0 \leq n \leq 5$までの範囲にし、そこに1を加えることで$1 \leq n \leq 6$にしようというものだ。残念ながら、この方法は偏った6面ダイスを作ってしまう。なぜか。

生の乱数rには3ビット以上の情報が必要だ。コンピューターは整数をビット列で表現するのですべてのビットが等しく乱数の場合、$n$ビットの乱数値は$2^n$個の状態を持つ。これを$0 \leq r \leq 2^n-1$に割り振った符号なし整数にしたとする。rが3ビットの場合、その値の範囲は$0 \leq r \leq 7$だ。

上のコードでは、0から5まではそのまま1から6になる。剰余のため、6と7はそれぞれ1と2になる。すべての取り得る乱数を書き出してみよう。

生の乱数 ダイス目

0 1 1 2 2 3 3 4 4 5 5 6 6 1 7 2

するとこの6面ダイスは1,2の出目の確率が$\frac{2}{8}$で、3,4,5,6の出目の確率が$\frac{1}{8}$になる。

よりビット数の大きな生の乱数を使ってもこの問題は解決しない。ビット数を増やせば増やすほど、偏りを減らすことはできるが、偏りは絶対になくならない。理由は、6の素因数3は2で割り切れないためだ。

では$1 \leq n \leq 6$までの乱数を得るにはどうするのかというと以下のようなアルゴリズムで分布を行う。

  1. 3bitの生の乱数rを得る
  2. rが$0 \leq r \leq 5$なら``r`+1'が分布された乱数
  3. それ以外の場合、1.に戻る

これを実装すると以下のようになる。

template < typename Engine >
auto dice( Engine & e )
{
    // ループを実行する
    while (true)
    {
        // 3bitの生の乱数を得る
        auto r = e() & 0b111 ;
        // 0-5なら乱数分布終わり
        if ( r <= 5 )
            return r + 1 ;
        // それ以外ならば振り直し
    }
}

この関数の実行時間は確率的に決まる。この実装はとても非効率的に見えるが、これ以外に公平に2で割り切れない素因数を含む状態数の乱数を生成する方法はない。

このコードは1回の乱数生成をキャッシュして複数回の3bitの乱数を切り出すなどの処理をすれば、乱数生成にコストがかかる場合、その分高速化できる。

分布クラス

分布クラスには共通の機能がある。本書ではすべてを解説しないが、重要な機能を解説する。

まず標準ライブラリの分布クラスに共通する機能を説明する。

分布クラスはオブジェクトを作り、そのオブジェクトを乱数エンジンと組み合わせて使う。その際、コンストラクターの引数で細かい設定を指定する。

distribution_type d( /*設定*/ ) ;

分布クラスはoperator ()を呼び出して乱数を分布させる。その際、引数には乱数エンジンへの非constなリファレンスを指定する

template < typename Engine, typename Distribution >
void f( Engine & e, Distribution & d )
{
    auto r = d(e) ;
}

rが分布された乱数。乱数エンジンeと乱数分布dは乱数を生成したので内部状態が変更される。

分布クラスはメンバー関数minとmaxで分布する乱数の最小値、最大値が得られる。

template < typename Distribution >
void f( Distribution & d )
{
    auto a = d.min() ;
    auto b = d.max() ;
}

分布クラスは構築時の実引数を同名のメンバー関数で取得することができる。

例えば、std::uniform_int_distribution( a, b )の場合、構築時に渡したa, bを引数の名前でメンバー関数として取得できる。

std::uniform_int_distribution d( 1, 6 ) ;
d.a() ; // 1
d.b() ; // 6

分布クラスは内部状態のリセットができる。

分布クラスは内部的に乱数値をキャッシュしている可能性がある。例えば乱数値が0か1である場合、1 bitの乱数しか必要ない。ここで渡した乱数エンジンが2 bit以上の乱数を生成できるのであれば、乱数値をキャッシュしておいて、1 bitずつ切り出して使うという最適化が考えられる。

しかしこの場合、同じ乱数エンジンを渡したのに、結果が違うということが起こり得る。

int main()
{
    std::uniform_int_distribution a( 1, 6 ) ;
    std::uniform_int_distribution b( 1, 6 ) ;

    std::mt19937 x ;

    // 乱数を生成
    // aは内部に乱数をキャッシュするかもしれない。
    a( x ) ;

    // yはxと同じ内部状態を持つ
    // つまり生成する生の乱数は同じ
    std::mt19937 y = x ;

    auto r1 = a( x ) ;
    auto r2 = b( y ) ;

    // r1 == r2 である保証はない
    
}

このような場合に、内部状態をリセットするメンバー関数resetを呼び出せば、同じ内部状態になることが保証される。

// 内部状態をリセット
a.reset() ;
// true
auto bool = ( a(x) == b(y) ) ;

また、この内部状態を取り出すこともできる。内部状態はネストされた型名param_typeで保持できる。内部状態を取り出すにはメンバー関数param()を呼び出す。分布クラスのコンストラクターにこのparam_typeの値を渡すと、同じ内部状態の分布クラスを作り出すことができる。またメンバー関数param(parm)でparam_typeの値を渡して内部状態を設定することも可能だ。

template < typename Distribution >
void f( Distribution & d )
{
    // 内部状態の取り出し
    // Distribution::param_type型
    auto p = d.param() ;

    // dと同じ内部状態を持つ変数
    Distribution same_d( p ) ;

    Distribution other ;
    // 既存の変数の内部状態を変更
    other.param( p ) ;
}

一様分布(Uniform Distribution)

一様分布とは乱数の取り得る状態がすべて等しい確率で出現する乱数のことだ。

整数の一様分布(std::uniform_int_distribution<IntType>)

uniform_int_distribution<IntType>は整数型の乱数$i$, $a \leq i \leq b$を以下の定数離散確率関数に従って分布させる。

$$ P(i,|,a,b) = 1 / (b - a + 1) \text{ .} $$

変数の宣言:

std::uniform_int_distribution<IntType> d( a, b ) ;

IntTypeは整数型でデフォルトはint、aは最小値、bは最大値。ただし$a \leq b$。

エンジンも含めた使い方は以下のとおり。

template < typename Engine >
void f( Engine & e )
{
    std::uniform_int_distribution d(1, 10) ;
    d.a() ; // 1
    d.b() ; // 10

    // 1から10までの範囲の乱数
    auto r = d(e) ;
}

値の範囲には負数も使える。

std::uniform_int_distribution d( -3, 3 ) ;

この分布は、$-3$, $-2$, $-1$, 0, 1, 2, 3のいずれかをそれぞれ$\frac{1}{7}$の等しい確率で返す。

浮動小数点数の一様分布(uniform_real_distribution<RealType>)

uniform_real_distribution<RealType>は浮動小数点数型の乱数$x$, $a \leq x &lt; b$を以下の定数確率密度関数に従って分布させる。

$$ p(x,|,a,b) = 1 / (b - a) \text{ .} $$

a == bのときは未定義となる。

変数の宣言:

std::uniform_real_distribution<RealType> d( a, b ) ;

RealTypeは浮動小数点数型でデフォルトはdouble、aは最小値、bは最大値。値の範囲は$a \leq b$ かつ $b - a \leq \text{RealType型の最大値}$。

エンジンも含めた使い方は以下のとおり。

template < typename Engine >
void f( Engine & e )
{
    std::uniform_real_distribution d(0.0, 1.0 ) ;
    d.a() ; // 0.0
    d.b() ; // 1.0

    // 0.0から1.0までの範囲の乱数
    auto r = d(e) ;
}

浮動小数点数の難しいところは、整数と違って値の範囲の状態が多いことだ。例えば0.0と1.0の間には0.5もあるし、0.01もあるし0.001もある。浮動小数点数の実装が表現できる状態はとても多い。uniform_real_distributionは指定された値の範囲で浮動小数点数が表現できるすべての値のうちから乱数を生成してくれる。そのため読者は難しいことを考える必要はない。

ベルヌーイ分布(Bernoulli distributions)

ベルヌーイ分布(bernoulli distribution)とは、ベルヌーイ試行(bernoulli trial)に関する分布だ。

ベルヌーイとは数学者ヤコブ・ベルヌーイ(Jacob Bernoulli, 1655-1705)に由来する。ヤコブ・ベルヌーイは西洋数学史上、ジェロラモ・カルダーノ(Gerolamo Cardano, 1501-1576)に続く二人目の、数学的に乱数をまともに文章に書き残した数学者で、現在の統計の基礎を切り開いた人物だ。西洋数学史において乱数と統計の研究は遅れた。この理由は主に宗教的なもので、運命とは神の決定したもうことであり、人の子の及ぶところではないとする考え方が一般的だった。そのため、まともな数学者は乱数を研究しなかった。ヤコブ・ベルヌーイの乱数に関する論文も、けっきょく本人は完成させることができず、論文完成に息子が着手しようとするも、これまた乱数はまともな数学者の取り組むべきところではないという周囲の圧力のために断念するなどの興味深い歴史がある。

ベルヌーイ試行

ベルヌーイ試行とは、独立した試行で結果が2種類のものだ。

「独立した試行」というのは、試行結果が前回の試行に影響されないことをいう。例えばコイントスの結果は表と裏だが、前回のコイントスの結果は今回のコイントスに影響しない。

結果が2種類というのは、試行をした結果、2種類の結果のうちのどちらか一方が出ることを言う。成功/失敗、表/裏、勝ち/負け、true/falseなど、なんでもいい。数学的には成功/失敗を使うが、C++ではtrue/falseで表現する。

ベルヌーイ試行において、確率$p$が成功する確率である場合、確率$q=1-p$が失敗の確率だ。

  • $p=1$の場合、必ず成功する。失敗しない
  • $p=0$の場合、必ず失敗する。成功しない
  • $p=0.5$の場合、成功と失敗は同じ確率になる
  • $p=0.4$の場合、$\frac{2}{5}$の確率で成功し、$\frac{3}{5}$の確率で失敗する

具体的なベルヌーイ試行の例を挙げると、

  • コイントスの結果、表か、裏か
  • 6面ダイスを振って6が出るか、6以外が出るか
  • 6面ダイスを振って5,6が出るか、1,2,3,4が出るか
  • 確率1%で当たるくじ引きの結果がアタリか、ハズレか
  • 赤玉と白玉が多数入ったツボの中身をよくかき混ぜ、玉を1つだけ取り出し、戻す。引いた玉の色が赤か、白か

このような結果を2種類に分けることができ、そのうちのどちらか一方だけが結果として出る、かつ1回1回が独立した試行をベルヌーイ試行と呼ぶ。

ベルヌーイ分布を使うと、一様分布ではない2値(true/false)の確率的な結果について乱数で得ることができる。例えば、ビデオゲームで宝箱を開けると32%の確率でアイテムが入っており、68%の確率で空っぽであるとする。これを一様分布で実装すると、以下のようになる。

// 宝箱にアイテムが入っている場合trueを返す
template < typename Engine >
bool open_chest( Engine & e )
{
    // 1から100までの整数の乱数を生成する
    std::uniform_int_distribution d(1, 100) ;
    // 32以下ならアイテムが入っている
    // 33以上ならば空っぽ
    return d(e) <= 32 ;
}

このようなコードを書くのは間違いの元だ。確率32%というのは32.0/100.0というdouble型の値で表現できる。この値だけ指定して、残りはライブラリに任せたい。そのようなときに使うのがベルヌーイ分布だ。

ベルヌーイ分布(std::bernoulli_distribution)

ベルヌーイ分布(bernoulli distribution)は1回のベルヌーイ試行の結果を乱数として返す。

std::bernoulli_distributionはbool型の乱数$b$を以下の離散確率関数に従って分布する。

$$ P(b,|,p) = \left{ \begin{array}{ll} p & \text{ if $b = true$, or} \\ 1 - p & \text{ if $b = false$.} \end{array}\right. $$

確率$p$でtrueが、確率$1-p$でfalseが返る。

変数の宣言:

std::bernoulli_distribution d( p ) ;

pはdouble型で、値の範囲は$0 \leq p \leq 1$。

使い方:

int main()
{
    std::bernoulli_distribution d( 0.5 ) ;
    d.p() ; // 0.5 ;

    std::mt19937 e ;
    // 乱数生成
    d(e) ;
}

bernoulli_distributionはテンプレートクラスではない。生成する乱数の型はboolだ。pはdouble型で確率$p$のことだ。値の範囲は$0 \leq p \leq 1$。

例えば前述の32%の確率でアイテムが入っている宝箱を実装するには以下のようになる。

template < typename Engine >
bool open_chest( Engine & e )
{
    std::bernoulli_distribution d( 32.0 / 100.0 )
    return d(e) ;
}

この関数open_chestは確率32%でtrueを、確率68%でfalseを返す。

本当にそうだろうか。確かめてみよう。

32%の確率でtrueになり、68%の確率でfalseになっているかどうかを確かめるには、大量の乱数を生成してtrue/falseをカウントし、それぞれ乱数を生成した数で割って割合を見ればよい。

int main()
{
    // 試行回数
    const int trial_count = 100 ;

    std::mt19937 e ;
    std::bernoulli_distribution d( 32.0 / 100.0 ) ;

    std::array<int, 2> result{} ;
    for ( int i = 0 ; i != trial_count ; ++i )
        // boolからintへの変換は
        // falseが0, trueが1
        ++result[ d(e) ] ;

    std::cout << "false: "sv << double(result[0]) / double(trial_count) * 100.0 << "%\n"sv
        << "true : "sv << double(result[1]) / double(trial_count) * 100.0 << "%\n"sv ;
}

これを実行してみると、筆者の環境では、

false: 72%
true : 28%

と出力された。少し違う。乱数なので試行回数が少なすぎる場合は、大きく偏ることもある。では試行回数を増やしてみよう。変数trial_countが試行回数だ。

200回試行すると、

false: 72.5%
true : 27.5%

まだダメだ。1000回試行してみよう。

false: 68.5%
true : 31.5%

だいぶ近くなった。

ちなみに、このままどんどん1万回、10万回と試行回数を増やしていっても、精度はそれほど上がらない。このことはヤコブ・ベルヌーイの研究と関わってくる。

std::bernoulli_distributionのコンストラクターに与えるdouble型の確率pの範囲は$0.0 \leq p \leq 1.0$だ。

$p=1.0$ならば常にtrue, $p=0.0$なら常にfalse、$p=0.5$ならばtrue/falseが一様分布する。

二項分布(std::binomial_distribution<IntType>)

二項分布(binomial distribution)は確率$p$で成功するベルヌーイ試行を$t$回行ったときに成功した回数$i$を乱数として返す。

具体的に例えると、100回コイントスをした結果出た表の数だ。コイントスは表と裏とそれぞれ 50%ずつの確率で出す。表を成功(true)とすると、つまり、$p=0.5$のベルヌーイ試行だ。100回コイントスをするというのは$t=100$だ。つまり、100回コイントスをした結果出た表の数というのは、100回ベルヌーイ試行した結果の成功数になる。この結果は、期待値としては50だが、ここで生成するのは乱数なので、50回出るわけではない。運が悪ければ1回も表が出ないこともあり得る。

6面ダイスを60回振った結果出た1の目の回数もそうだ。この場合、$p=\frac{1}{6}$のベルヌーイ試行を$t=60$回行うことになる。期待値は$10$だ。成功を2,3の目が出た回数と考えることもできる。この場合期待値は$20$だ。

確率1%で当たるくじを100回引いた場合もそうだ。この場合、$p=0.01$で$t=100$になる。期待値は$1$なので、1回当たることが平均的に期待できる。ちなみに、実際に100回くじ引きをして1回でも当たる確率は約63%だ。

std::binomial_distribution<IntType>はIntType型の乱数$i \geq 0$を以下の離散確率関数に従って分布する。

$$ P(i,|,t,p) = \binom{t}{i} \cdot p^i \cdot (1-p)^{t-i} \text{ .} $$

変数の宣言:

std::binomial_distribution<IntType> d( t, p ) ;

IntTypeは整数型でデフォルトはintだ。tはIntType型の整数値で、値の範囲は$0 \leq t$だ。pはdouble型の値で確率を指定する。pの値の範囲は$0 \leq p \leq 1$だ。

使い方:

int main()
{
    std::binomial_distribution d( 1, 0.5 ) ;
    d.t() ; // 1
    d.p() ; // 0.5

    std::mt19937 e ;
    // 乱数生成
    d(e) ;
}

100回コイントスをした結果、表が出た回数を乱数で得る関数coinflips100は以下のように書ける。

template < typename Engine >
auto coinflips100( Engine & e )
{
    // t == 100, p == 0.5
    std::binomial_distribution d( 100, 0.5 ) ;
    return d(e) ;
}

100回のベルヌーイ試行をするので$t=100$で、ベルヌーイ試行の成功確率は$p=\frac{1}{2}=0.5$になる。

これを10回ぐらい呼んでみよう。

int main()
{
    std::mt19937 e ;
    for ( int i = 0 ; i != 10 ; ++i )
        std::cout << coinflips100( e ) << ", "sv ;
}

筆者の環境では結果は以下のようになった。

53, 54, 43, 56, 51, 50, 45, 48, 49, 47,

期待値は50なので、50前後の乱数が出やすい。

6面ダイスを60回振った結果出た1の目の合計を乱数で返す関数roll_for_oneは以下のようになる。

template < typename Engine >
auto roll_for_one( Engine & e )
{
    // t == 60, p == 1.0/ 6.0
    std::binomial_distribution d( 60, 1.0 / 6.0 ) ;
    return d(e) ;
}

60回のベルヌーイ試行をするので$t=60$で、ベルヌーイ試行の確率は6面ダイスの1の目が出る確率なので、$p=\frac{1}{6}$になる。

確率1%で当たるくじを100回引いた結果アタリの回数を返す関数lootboxは以下のように実装できる。

template < typename Engine >
auto lootbox( Engine & e )
{
    // t == 100, p = 0.01
    std::binomial_distribution d( 100, 1.0 / 100.0 ) ;
    return d(e) ;
}

この関数を10回呼び出してみると結果は以下のようになる。

1, 0, 2, 1, 0, 0, 0, 0, 1, 3,

確率1%で当たるくじを100回引くと、複数回当たることもあれば、1回も当たらないこともある。期待値は1だが、期待値というのは平均的に期待できる結果でしかない。読者諸君もくじ引きをするときは確率に気を付けよう。たとえくじが毎回公平であったとしても、確率は無記憶性なのだ。「もう90回くじを引いたからあと10回引けば当たるはず」という考え方は通用しない。

幾何分布(std::geometric_distribution)

幾何分布(geometric distribution)とは、確率$p$で成功するベルヌーイ試行を初めて成功するまで行った回数を乱数として分布する。

具体的な例で例えると、

  • コイントスを表が出るまで行った回数
  • 6面ダイスを1の目が出るまで振った回数
  • 確率1%で当たるくじ引きをアタリが出るまで引いた回数

コイントスの例で考えよう。コイントス1回をベルヌーイ試行とし、成功を表とする。表が出るまでコイントスをしてみよう。コイントスを何回する必要があるだろうか。運がよければ1回で表が出るので1回だ。運が悪ければ、5回コイントスをしても全部裏なこともあるだろう。100回コイントスをして表が一度も出ないことは、確率的にはあり得る。ただしその確率は$\frac{1}{2^{100}}$なので、およそあり得ない確率ではある。

std::geometric_distribution<IntType>はIntType型の乱数$i$, $i \geq 0$を以下の離散確率関数に従って分布する。

$$ P(i,|,p) = p \cdot (1-p)^{i} \text{ .} $$

変数の宣言:

std::geometric_distribution<IntType> d( p ) ;

IntTypeは整数型でデフォルトはint、pは確率で値の範囲は$0 < p < 1$だ。pの値の範囲に注意すること。0と1であってはならない。幾何分布は成功するまでベルヌーイ試行した回数を返すので、$p=0$の場合、必ず失敗するベルヌーイ試行になり意味がない。$p=1$のときは必ず成功するベルヌーイ試行であり、やはり意味がない。

geometric_distributionの生成する乱数の範囲にも注意が必要だ。生成される乱数$i$の範囲は$i \geq 0$だ。0もあり得る。0ということは、最初のベルヌーイ試行が成功したということだ。1は2回目のベルヌーイ試行が成功したということだ。幾何分布はベルヌーイ試行が初めて成功するまでのベルヌーイ試行の回数を返すので、成功したベルヌーイ試行は回数に含めない。

使い方:

int main()
{
    // p == 0.5
    std::geometric_distribution d( 0.5 ) ;
    d.p() ; // 0.5 ;

    std::mt19937 e ;
    d(e) ;
}

コイントスを表が出るまで繰り返し、その合計回数を乱数で返す関数try_coinflipsを書いてみよう。

template < typename Engine >
auto try_coinflips( Engine & e )
{
    std::geometric_distribution d( 0.5 ) ;
    return d(e) + 1;
}

最後に+1しているのは、この文脈では表を出したときのコイントスも数えるからだ。つまり成功したベルヌーイ試行も回数に数えるので、幾何分布の生成する乱数より1多い数になる。

10回呼び出してみたところ、以下のような戻り値を返した。

1, 3, 6, 1, 1, 2, 1, 8, 9, 5,

運がよければ1回で表が出るが、運が悪ければ9回かかる。もちろんもっとかかる可能性もある。

6面ダイスを1の目が出るまで振り、その合計回数を返す関数try_rollsを書いてみよう。

template < typename Engine >
auto try_rolls( Engine & e )
{
    std::geometric_distribution d( 1.0 / 6.0 ) ;
    return d(e) + 1;
}

これも10回呼び出してみると筆者の環境では以下のようになった。

1, 10, 20, 2, 3, 5, 2, 28, 31, 19,

6面ダイスを振ると、運がよければ1回で1の目が出るが、運が悪いと何十回も振る必要がある。

確率1%のくじを初めて当たるまで引き続け、くじを引いた回数を返す関数try_lootboxesも書いてみよう。

template < typename Engine >
unsigned int try_lootboxes( Engine & e )
{
    std::geometric_distribution d( 1.0 / 100.0 ) ;
    return d(e) + 1;
}

10回呼び出してみよう。

15, 180, 346, 25, 37, 79, 21, 493, 562, 342,

確率1%のくじを当てるには、運が悪いと何百回も引かなければならない。

負の二項分布(std::negative_binomial_distribution)

負の二項分布(negative binomial distribution)は幾何分布に似ている。幾何分布がベルヌーイ試行が1回成功するまでに行ったベルヌーイ試行の回数を乱数として分布するのに対し、負の二項分布はベルヌーイ試行が$k$回成功するまでに行ったベルヌーイ試行の回数を乱数として分布する。

負の二項分布を具体的な例で考えよう。

  • コイントスを、10回、表が出るまで行った回数
  • 6面ダイスを、10回、1の目が出るまで振った回数
  • 確率1%で当たるくじ引きを、10回、アタリが出るまで引いた回数

幾何分布は負の二項分布で表現することもできる。

  • コイントスを、1回、表が出るまで行った回数
  • 6面ダイスを、1回、1の目が出るまで振った回数
  • 確率1%で当たるくじ引きを、1回、アタリが出るまで引いた回数

std::negative_binomial_distribution<IntType>はIntType型の乱数$i$, $i \geq 0$を以下の離散確率関数に従って分布する。

$$ P(i,|,k,p) = \binom{k+i-1}{i} \cdot p^k \cdot (1-p)^i \text{ .} $$

$p = 1$のときの$P(i,|,k,p)$は未定義だ。

変数の宣言:

std::negative_binomial_distribution<IntType> d( k, p ) ;

IntTypeは整数型でデフォルトはint、kはIntType型の値$0 < k$で成功させるベルヌーイ試行の回数、pはdouble型の確率$- < p \leq 1$だ。

使い方:

int main()
{
    // k == 1, p == 0.5
    std::negative_binomial_distribution d( 1, 0.5 ) ;
    d.k() ; // 1
    d.p() ; // 0.5

    std::mt19937 e ;
    d(e) ;
}

幾何分布と同じく、負の二項分布が生成する乱数iはk回のベルヌーイ試行を成功させるまでに失敗したベルヌーイ試行の数を返す。

例えば、コイントスで10回表が出るまでに失敗したコイントスの数を返す。コイントスがベルヌーイ試行で、表が成功だ。成功したベルヌーイ試行の数は返さない。そのため、結果の乱数は10以下、時には0であることもあり得る。0というのは10回コイントスをしたらすべて表になったので1回もベルヌーイ試行が失敗しなかった場合だ。

コイントスを10回表が出るまでに行ったコイントスの回数を返す関数count_10_coinflipsは以下のように書く。

// 10回表が出るまでに行ったコイントスの数
// 表が出たコイントスも含める
template < typename Engine >
auto count_10_coinflips( Engine & e )
{
    std::negative_binomial_distribution d( 10, 0.5 ) ;
    return d(e) + 10 ;
}

最後に+10しているのは、この関数は成功も含めたコイントスの回数を返すからだ。10回表が出るまでに失敗した、つまり裏になったコイントス回数がほしければ、そのままの値を使えばいい。

// 10回表が出るまでに行った失敗したコイントスの数
// 表が出たコイントスは含めない
template < typename Engine >
auto count_failed_coinflips_until_10_heads( Engine & e )
{
    std::negative_binomial_distribution d( 10, 0.5 ) ;
    return d(e) ;
}

参考までに、n回表が出るまでに行ったコイントスの回数を乱数で返す関数count_n_coinflipsは以下のとおり。

template < typename Engine >
auto count_n_coinflips( unsigned int n, Engine & e )
{
    std::negative_binomial_distribution d( n, 0.5 ) ;
    return d(e) + n ;
}

6面ダイスを10回、1の目が出るまで振った回数を乱数で返す関数count_10_rollsは以下のとおり。

template < typename Engine >
auto count_10_rolls( Engine & e )
{
    std::negative_binomial_distribution d( 10, 1.0/6.0 ) ;
    return d(e) + 10 ;
}

確率1%のくじを10回当てるまでくじを引いた回数を返す関数count_10_lootboxesは以下のとおり。

template < typename Engine >
auto count_10_lootboxes( Engine & e )
{
    std::negative_binomial_distribution d( 10, 0.01 ) ;
    return d(e) + 10 ;
}