はじめに

こんにちは、watnowアドベントカレンダー20日目担当のゆいぴです。私ごとですが最近コロナに罹患して咳やばすぎて泣きながらQiita書いてます。めっちゃしんどい。

まあでもそれはどうでもいいので本題に入ります。

シミュレータとの出会い

この前から研究室の方針的に、シミュレータで使う乱数の発生方法について勉強することになったのですが、わたくしの貧相な頭では理解出来なさすぎて研究室の先輩を小一時間捕まえる羽目になりました。先輩ありがとう、、、、
その時に自分なりに理解して咀嚼したことを覚書として書こう！というのが今回のモチベーションです。初学者のざっくり理解とざっくり解説なのでお手柔らかに。。

確率分布

一様分布、指数分布、Zipf分布の３つの分布について今回は紹介しようと思います。

一様分布

世の中には確率に従うものがたくさんあります。例えば、一般的なサイコロを振ってそれぞれの目が出る確率は1/6ですね。サイコロのように、どの目も絶対に同じ確率の場合（一様分布）は、考え方が簡単です。

確率質量関数f(x)は、xの値をとる確率、つまりどの値でも1/6です。
確率分布関数（累積分布関数）は、x以下の値をとる確率なので、例えばサイコロで2以下の目が出る確率は2/6で1/3になりますね。確率分布関数は、xまでのすべての確率質量関数の和になります。

シミュレータとは、例えばこのサイコロのすべての目が1/6の確率で出ることを実際に実験して確かめることです。
確率質量関数の求め方は、10000回サイコロを投げてそれぞれの目が何回出たのか回数を数えて、全体の試行回数で割ればいいですね。

簡単に下でC言語のシミュレータを載せておきます。（今回一様乱数の発生のためにメルセンヌ・ツイスタの関数を利用してgenrandを使っています。これについては割愛します。）

simuration.c

#include <stdio.h>
#include <stdlib.h>
#include "MT.h"

int num;
int count[6] = {0};
double probability[6] = {0};

int main(void){
    for(int i = 0; i < 10000; i++){
        num = (int)(genrand_real2()*6.0) //0から5の乱数発生
        count[num]++; //出た目の数を数える
    }
    for(int i = 0; i < 6; i++){
        probability[i] = (double)count[i] / 10000; //試行回数で出た回数を割る
    }
}

指数分布

指数関数とは、ランダムに発生する事象が生じる発生時間間隔などに用いられる確率分布です。しかし、genrandでは一様な乱数しか発生できないので指数分布に従う乱数を作ることは難しそうですね。こういう時に使える考え方が逆関数法です。
確率分布関数F(x)は、0から1の範囲をとります（確率の総和だから）。なのでy = F(x)の逆関数F^-1(y)を考えて、yに0から1をとる一様乱数をねじ込めば出てくるxはこの確率分布関数に従うものになります。

指数関数

F(x) = 1 - e^-λx

これを逆関数にすると

x =  -  \frac{log(1 - F(x))}{λ}

になります。

ここでのλは平均発生間隔なので、任意に設定することができます。F(x)の部分に一様乱数genrand_real2()を書くことで、指数分布に従う乱数を発生させることができます。
指数分布に従う乱数(指数乱数)を、逆関数法により1,000,000回発生させた時の頻度分布を作成するプログラムを例として載せておきます。平均発生間隔は10、値の間隔1、区間数100にしました。

exponential.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "MT.h"

int main(void)
{
    int rambda = 10; // 平均発生間隔
    double x = 0.0;                     // 指数乱数
    int n = 100;                        // 調べる範囲
    int sample = 1000000;               // サンプル数(1,000,000回)
    int count[100] = {0};               // 回数を数える配列
    double probability[100] = {0.0};    // countを元に頻度分布を入れる配列

  
    for (int j = 0; j < sample; j++)
    {
        x = -log(1 - genrand_real2()) * rambda;
            
    if (x > n)
    {
        count[n]++;
    }
    else
    {
        for (int k = 0; k < n; k++)
        {
            if ((int)x == k)
            {
                count[k]++;
            }
        }
    }

    for (int j = 0; j < n; j++)
    {
        probability[j] = (double)count[j] / sample;
    }
}

Zipf分布

指数分布に比べて全然聞きなじみない名前で最初に見た時なんだこれ？と思ったんですけど、コンテンツの要求発生数などが従う分布らしいです。もう少し噛み砕いていうと、YouTubeなどでコンテンツを人気順に並べたときに、そのコンテンツがそれぞれどのくらいの確率で要求されるか(つまり、どれくらいの確率で見たい！とクリックされるか)が従う分布です。当たり前ですが、人気のコンテンツは要求される確率が高くなりますし、人気のないコンテンツはあんまり要求されない、そんな分布だと考えてください。
この分布を紹介したのは、実はこの分布は逆関数法が使えないからです。逆関数法は離散的な値のみをとる場合や、連続的でも確率分布関数が得られないなどの条件下では使うことができません。Zipf分布は離散的な値をとるため、逆関数法が使えないのです。あーあ。

じゃあどうしよう？という話ですが、やるべき手順は決まっています。

1.　各離散値が生じる確率を計算し，分布関数の値を計算し，配列に用意 : 一番人気のコンテンツの要求される確率がxx、二番人気のコンテンツの要求される確率がxx、...とそれぞれのコンテンツの要求される確率を求め、「分布関数」を求めます。
2.　0から1の一様乱数を発生させて、生じた値>分布関数の値となる最小のところの値が発生したことにする: 例えば、分布関数がP(0) = 0.20,P(1) = 0.35,P(2) = 0.45,...となっているとき、乱数で0.30が発生したら0番目のコンテンツ、0.38が出たら1番目のコンテンツが要求されたとみなす、といった考え方です。

これで逆関数法が使えない分布に従うシミュレータも作ることができます。個人的には高校の時にやったガウス記号のグラフになんか雰囲気似てるな〜と思います。0.20~0.35の間は何が出ても0！0.35~0.45の間は何が出ても1！って感じで数字は対応一致しないんでほんとに雰囲気だけですけどね。

Zipf分布において、k番目のコンテンツに対する要求発生確率は

p_k = \frac {1/k^\theta}{\left( \sum_{j=1}^M 1/j^\theta \right)}

で表すことができます。θは、偏りを表すパラメータで大きいほど偏りが大きいことを示します。
この確率に基づいて、分布関数を求め、それに基づいて一様乱数の値を固定化することでZipf分布に従う乱数も作成することができます！

Zipf関数に従う乱数の頻度分布も作成しました。

Zipf.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "MT.h"

int main(void)
{

    int M = 100;       // コンテンツ数
    int try = 1000000; // 試行回数
    double x = 0.0;

    int count[100] = {0};
    double n_theta[100] = {0.0};     // 分布関数の分子を考える
    double j_theta = {0.0};          // 分布関数の分母を考える
    double zipf_p[100] = {0.0};      // 分布関数に基づく離散値
    double theta = 0.8;              // 偏りを表すパラメータ
    double probability[100] = {0.0}; // 頻度分布

    // 累積分布関数zipf_p[k]を求める
    n_theta[0] = 1.0 / pow(1.0, theta);
    for (int k = 1; k < M; k++)
    {
        n_theta[k] = n_theta[k - 1] + 1.0 / pow((double)k + 1, theta);
    }
    for (int j = 1; j <= M; j++)
    {
        j_theta += 1 / pow((double)j, theta);
    }

    for (int k = 0; k < M; k++)
    {
        zipf_p[k] = n_theta[k] / j_theta;
    }

    // 頻度分布 probabilityを求める
    for (int j = 0; j < try; j++)
    {
        x = genrand_real2();
        for (int k = 0; k < M; k++)
        {
            if (x < zipf_p[k])
            {
                count[k]++;
                break;
            }
        }
    }
    for (int j = 0; j < M; j++)
    {
        probability[j] = (double)count[j] / try;
    }
}

おわりに

ほぼほぼ備忘録的な感じになってしまって誰のためになるん！って感じですが、シミュレータはこんな感じで使いたい乱数をうまく使えるようにどうにかこうにかして、使える形に整えてあげて、それからなんとかしよう！みたいなノリらしいです。むずいよ〜〜。
でもなんとなくシミュレータへの理解は進んだので卒業研究頑張ります（？）、何かを求めてこの記事にやってきたあなたも一緒に頑張りましょう。それでは。