このエントリーをはてなブックマークに追加

 

C++11 乱数 std::random 入門
Copyright (C) 2014 by Nobuhide Tsuda

 

※ イラスト提供:かわいいフリー素材集「いらすとや」様

rand() の問題点

これまで広く使用されてきた rand() による乱数生成には以下の問題点がある。

  • 生成される範囲が [0, 32767] と狭い
  • rand() % N は一様では無い
  • 周期があまり長くない
  • 乱数生成アルゴリズムが固定(通常は線形合同法)
  • 正規分布など、一様でない乱数生成が面倒

rand() で生成される値の範囲は [0, 32767] と狭い。15ビットしかない。
(rand() << 15) + rand() とすれば、30ビット乱数にすることは可能だが、スマートではないし、 生成アルゴリズムの関係で、乱数に偏りが出る場合がある。

通常、乱数として欲しいのは [0, N) 範囲の値であることが多い。この場合は rand() % N とすることが一般的だ。 しかし、このようにして生成した乱数は一様でなく偏りがある。特に N が大きい時にその現象が顕著になる。
例えば、N が1万だとすると、rand() % N は [0, 2768) の範囲の値を返す確率は [2768, 10000) の範囲の値を返す確率より有意に大きくなる。 具体的には 4:3 の比率になる。100万回乱数を発生させるとすると、前者の期待値は 120回, 後者の期待値は 90回となる。
当然、この問題を回避するコードは記述可能なのだが、ライブラリがそれを用意するべきだ。

rand() の周期は 2^31 である。充分な長さだと思う読者もいるかもしれないが、適用分野によっては短い場合もある。

通常、rand() は線形合同法により乱数を生成する。このアルゴリズムは原始的であり、高速だが、あまりよい性質の乱数を生成してくれない。
場合によっては乱数系列を推測されてしまう危険性もある。
オンラインポーカーのシャフルにこの乱数を使い、もしユーザに次のカードを予測されてしまっては、大問題だ。

rand() は一様な乱数しか生成してくれない。正規分布の乱数を生成させたい場合は、そのような分布になるよう、プログラムを記述しなくていけない。

上記のような問題点を解消するために、C++11 で乱数ライブラリが追加された。
Visual Studio に於ける本格的なC++11対応は VS2013 からであるが、auto など一部の仕様は VS2010 から利用可能になっている。 嬉しいことに、この random ライブラリも利用可能なので、VS2010 以上を使用しているのであれば random ライブラリ使わなきゃ損だぜ。

演習問題:解答例

  1. rand() % 10000 により [0, 10000) の乱数を100万回発生させ、0, 1, 2, 9997, 9998, 9999 が出る回数をカウントし、表示しなさい。
  2. 同様に、rand() % 1000 の場合はどうなるか検証してみなさい。

C++ 乱数ライブラリ std::random の使い方

準備

乱数ライブラリを使用するには、以下のように <random> をインクルードする。

#include <random>

非決定的な乱数生成

非決定的な乱数生成器として「std::random_device」が用意されている。
「非決定的」とは、次の値が一意に決まらないという意味だ。内部的にはハードウェアの状態や時刻を用いて実現している。

使い方は「std::random_device rnd;」で乱数生成器(rnd は乱数生成器の名前、お好みの名前を付けてよい)を生成し、 「rnd()」をコールすると、乱数を生成してくれる。生成される乱数は 32bit で、[0, 0xffffffff] の範囲だ。

※ rnd は std::random_device クラスのオブジェクトなのに、それをあたかも関数のように「()」を付けてコール出来るの? と不思議に思った読者もいるかもしれない。 実はクラス定義で () 演算子(oprator()())をオーバーロードすることで、オブジェクト名() でその演算子をコールすることが出来るのだ。 これを「ファンクタ」と呼ぶ。これは初級者にはわかりづらいかもしれないので、よくわからない場合は、あまりこだわらず、そういうものだと思い、 先に進んで欲しい。

#include <random>
int main()
{
    std::random_device rnd;     // 非決定的な乱数生成器
    for (int i = 0; i < 10; ++i) {
        std::cout << std::hex << rnd() << "\n";
    }
    getchar();
    return 0;
}

上記は10個の乱数を生成し、表示するサンプルだ。

演習問題:

  1. 上記ソースコードをビルド・実行し、結果を確かめなさい。

決定的な乱数(擬似乱数)生成

「std::random_device」は若干速度が遅いので、通常は決定的な乱数生成器を使用する。
「決定的」とは、現在の値により次の値が一意に決まるという意味。擬似乱数とも呼ばれ、数学的に完全な乱数ではない。
が、rand() の srand(int) 同様、乱数のシードを指定できるので、決定的でも実用上問題は無い

決定的な乱数生成アルゴリズムとしては、線形合同法メルセンヌ・ツイスタなどがあり、選択することができる。
線形合同法は大昔から使われていたアルゴリズムで、単純・高速ではあるが、乱数の質はあまり良くない。
メルセンヌ・ツイスタは1990年代後半に日本人により開発されたもので、周期が 2^219937-1 もあるなど、高品質の乱数を高速に生成するように設計されている。 かなりお薦めだ。

本稿では、この メルセンヌ・ツイスタ を用いた乱数生成について説明する。

まずは、シード生成のために「std::random_device rnd;」のようにして、非決定的乱数生成器を生成する。 これまで通り、time(0) を使ってもよい。が、お薦めはしない。

これを引数として指定し、「std::mt19937 mt(rnd());」の様にメルセンヌ・ツイスタ生成器を生成する(mt は生成器の名前。自由に命名してよい)。
後は「mt()」により乱数を生成するだけだ。

生成される乱数の範囲は 32ビット、[0, 0xffffffff] の範囲だ。
64ビット範囲の乱数を生成したい場合は std::mt19937_64 を使用する。

    std::random_device rnd;     // 非決定的な乱数生成器
    std::mt19937 mt(rnd());            // メルセンヌ・ツイスタの32ビット版、引数は初期シード
    //std::mt19937 mt((int)time(0));        //  メルセンヌ・ツイスタの32ビット版、引数は初期シード
    for (int i = 0; i < 10; ++i) {
        std::cout << std::hex << mt() << "\n";
    }

シードの指定はコンストラクタだけでなく、seed(val) メソッドにより指定することも可能。
下記はその例。先にデフォルトコンストラクタでオブジェクトを構築し、後でシードを指定している。
自分で定義するクラスのメンバ変数として std::mt19937 を保持する場合は、この方法を使いコンストラクタなどでシードを指定するとよい。

    std::mt19937 mt;            // メルセンヌ・ツイスタの32ビット版
    std::random_device rnd;     // 非決定的な乱数生成器
    mt.seed(rnd());            // シード指定
    for (int i = 0; i < 10; ++i) {
        std::cout << std::hex << mt() << "\n";
    }

演習問題:

  1. 上記ソースコードをビルド・実行し、結果を確かめなさい。
  2. 上記ソースコードで、初期シードを指定しないようにして、ビルド・実行し、結果を確かめなさい。
  3. mt19937 を mt19937_64 に変更・ビルド・実行し、結果を確かめなさい。

範囲指定乱数生成

これまで、std::random_device、std::mt19937 を生で使ってきたが、これらは範囲が固定だった。
rand() の時と同様に「% N」演算を行うことで [0, N) の範囲の(ある程度)一様乱数を得ることは可能だが、 もっと正確な一様乱数を得る仕組みが用意された。それが「uniform_int_distribution」だ。

「std::uniform_int_distribution<> rand100(0, 99);」のように宣言すると、[0, 99] の範囲の一様分布器を生成できる。
これを先の乱数生成器と組み合わせ、「rand100(mt)」のように記述することで一様分布乱数を生成できる。
※ 引数が「mt()」ではなく「mt」であることに注意されたい。生成された乱数値を渡しているではなく、乱数生成器自体を渡している。 乱数生成器が持つ情報を参照可能だから、正確に一様分布する乱数を発生可能なのである。

    std::random_device rnd;     // 非決定的な乱数生成器を生成
    std::mt19937 mt(rnd());     //  メルセンヌ・ツイスタの32ビット版、引数は初期シード値
    std::uniform_int_distribution<> rand100(0, 99);        // [0, 99] 範囲の一様乱数
    for (int i = 0; i < 20; ++i) {
        std::cout << rand100(mt) << "\n";
    }

上記は、初期シードに非決定的乱数生成器で生成した値を使ったメルセンヌ・ツイスタ乱数生成器と、 [0, 100) の範囲の一様分布器を組み合わせて、20個の乱数を生成する例。

整数ではなく、double の一様分布器は「uniform_real_distribution」で生成する。

一様乱数ではなく、「normal_distribution」を使って、指定された平均値、分散値の正規分布乱数を生成することもできる。
コードは以下の様になる。

    std::random_device rnd;     // 非決定的な乱数生成器でシード生成機を生成
    std::mt19937 mt(rnd()); //  メルセンヌツイスターの32ビット版、引数は初期シード
    //std::uniform_int_distribution<> rand100(0, 99);     // [0, 99] 範囲の一様乱数
    std::normal_distribution<> norm(50.0, 10.0);       // 平均50, 分散値10の正規分布
    for (int i = 0; i < 20; ++i) {
        //std::cout << rand100(mt) << "\n";
        std::cout << norm(mt) << "\n";
    }ma

演習問題:解答例

  1. 上記ソースコードをビルド・実行し、結果を確かめなさい。
  2. 上記ソースコードを、[1.0, 2.0] の範囲の一様なdouble型乱数を20個表示するよう修正・ビルド・実行してみなさい。
  3. 上記ソースコードを、平均 0.0、分散 0.2 の正規分布乱数を20個表示するように修正してみなさい。
  4. [0, 9999] の範囲の一様な整数乱数を100万回生成し、0, 1,2 と 9997, 9998, 9999 の値が出る回数がほぼ等しいことを確認しなさい。
  5. rand() % 100 と、本稿の両方の方法で [0, 99] の範囲の一様な乱数を1000回生成し、平均・標準偏差を計算し、それらを比較してみなさい。

まとめ

ここに紹介した以外にも、random ライブラリには便利なクラス・機能がたくさん用意されている。 詳細は cpprefjp - random(C++11) や、 <random> などを読んで自分で勉強してね。

これまでの、srand(), rand() に慣れている人には、ちょっと面倒に感じるかもしれないが、
慣れればたいしたことないし、パズルの問題生成・カードシャフル・敵機の動きなど、ちゃんとした乱数にした方がよいことは確かなので、 random ライブラリ使いこなせるようになっておくれやす。


演習問題解答例

 

  rand() の問題点

  1. rand() % 10000 により [0, 10000) の乱数を100万回発生させ、0, 1, 2, 9997, 9998, 9999 が出る回数をカウントし、表示しなさい。
  2. #include <iostream>
    #include <time.h>
    #include <vector>
    
    int main()
    {
        srand((int)time(0));
        
        const int N = 10000;      // [0, N) の乱数を発生
        std::vector<int> vc(N, 0);     // 回数を保存する vector
        for (int i = 0; i < 1000*1000; ++i) {       //  百万回ループ
            int r = rand() % N;     // 乱数生成
            vc[r] += 1;           //  生成された乱数の個数をインクリメント
        }
        std::cout << "vc[0] = " << vc[0] << "\n";
        std::cout << "vc[1] = " << vc[1] << "\n";
        std::cout << "vc[2] = " << vc[2] << "\n";
        std::cout << "vc[" << N-3 << "] = " << vc[N-3] << "\n";
        std::cout << "vc[" << N-2 << "] = " << vc[N-2] << "\n";
        std::cout << "vc[" << N-1 << "] = " << vc[N-1] << "\n";
        getchar();
        return 0;
    }
    
  3. 同様に、rand() % 1000 の場合はどうなるか検証してみなさい。
  4. #include <iostream>
    #include <time.h>
    #include <vector>
    
    int main()
    {
        srand((int)time(0));
        
        const int N = 1000;      // [0, N) の乱数を発生
        std::vector<int> vc(N, 0);     // 回数を保存する vector
        for (int i = 0; i < 1000*1000; ++i) {       //  百万回ループ
            int r = rand() % N;     // 乱数生成
            vc[r] += 1;           //  生成された乱数の個数をインクリメント
        }
        std::cout << "vc[0] = " << vc[0] << "\n";
        std::cout << "vc[1] = " << vc[1] << "\n";
        std::cout << "vc[2] = " << vc[2] << "\n";
        std::cout << "vc[" << N-3 << "] = " << vc[N-3] << "\n";
        std::cout << "vc[" << N-2 << "] = " << vc[N-2] << "\n";
        std::cout << "vc[" << N-1 << "] = " << vc[N-1] << "\n";
        getchar();
        return 0;
    }
    
範囲指定乱数生成
  1. 上記ソースコードをビルド・実行し、結果を確かめなさい。
  2. 解答例省略

  3. 上記ソースコードを、[1.0, 2.0] の範囲の一様なdouble型乱数を20個表示するよう修正・ビルド・実行してみなさい。
  4.     std::random_device rnd;     // 非決定的な乱数生成器を生成
        std::mt19937 mt(rnd());     //  メルセンヌ・ツイスタの32ビット版、引数は初期シード値
        std::uniform_real_distribution<> rand12(1.0, 2.0);    // [1.0, 2.0] 範囲の一様乱数
        for (int i = 0; i < 20; ++i) {
            std::cout << rand12(mt) << "\n";
        }  
    
  5. 上記ソースコードを、平均 0.0、分散 0.2 の正規分布乱数を20個表示するように修正してみなさい。
  6.     std::random_device rnd;     // 非決定的な乱数生成器を生成
        std::mt19937 mt(rnd());     //  メルセンヌ・ツイスタの32ビット版、引数は初期シード値
        std::normal_distribution<> norm(0.0, 0.2);       // 平均0.0, 分散値0.2の正規分布
        for (int i = 0; i < 20; ++i) {
            std::cout << norm(mt) << "\n";
        }  
    
  7. [0, 9999] の範囲の一様な整数乱数を100万回生成し、0, 1,2 と 9997, 9998, 9999 の値が出る回数がほぼ等しいことを確認しなさい。
  8. #include <iostream>
    #include <vector>
    int main()
    {
        const int N = 10000;     //  [0, N) の乱数を発生
        std::random_device rnd;     // 非決定的な乱数生成器を生成
        std::mt19937 mt(rnd());     //  メルセンヌ・ツイスタの32ビット版、引数は初期シード値
        std::uniform_int_distribution<> randN(0, N-1);        // [0, N) 範囲の一様乱数
        std::vector vc(N, 0);     // 回数を保存する vector
        for (int i = 0; i < 1000*1000; ++i) {     // 百万回ループ
            int r = randN(mt);      // 乱数生成
            vc[r] += 1;         //  生成された乱数の個数をインクリメント
        }
        std::cout << "vc[0] = " << vc[0] << "\n";
        std::cout << "vc[1] = " << vc[1] << "\n";
        std::cout << "vc[2] = " << vc[2] << "\n";
        std::cout << "vc[" << N-3 << "] = " << vc[N-3] << "\n";
        std::cout << "vc[" << N-2 << "] = " << vc[N-2] << "\n";
        std::cout << "vc[" << N-1 << "] = " << vc[N-1] << "\n";
        getchar();
        return 0;
    }
    
  9. rand() % 100 と、本稿の両方の方法で [0, 99] の範囲の一様な乱数を1000回生成し、平均・標準偏差を計算し、それらを比較してみなさい。
  10. #include <iostream>
    #include <time.h>
    #include <vector>
    #include <random>
    int main()
    {
        const int N = 1000;
        {
            srand((int)time(0));
            int sum = 0;      // Σx
            int sum2 = 0;    // Σx^2
            for (int i = 0; i < N; ++i) {
                int r = rand() % 100;
                sum += r;
                sum2 += r * r;
            }
            double ave = (double)sum / N;
            std::cout << "ave = " << ave << "\n";
            std::cout << "σ = " << sqrt((double)sum2/N - ave * ave) << "\n";
        }
        {
            int sum = 0;      // Σx
            int sum2 = 0;    // Σx^2
            std::random_device rnd;     // 非決定的な乱数生成器を生成
            std::mt19937 mt(rnd());     //  メルセンヌ・ツイスタの32ビット版、引数は初期シード値
            std::uniform_int_distribution<> rand100(0, 99);        // [0, 99] 範囲の一様乱数
            for (int i = 0; i < N; ++i) {
                int r = rand100(mt);
                sum += r;
                sum2 += r * r;
            }
            double ave = (double)sum / N;
            std::cout << "ave = " << ave << "\n";
            std::cout << "σ = " << sqrt((double)sum2/N - ave * ave) << "\n";
        }
        getchar();
        return;
    }