範囲指定付き疑似ガウス乱数

2020/07/07 in 逆引き PROCESSING
python processing プログラミングジェネラティブアート数学乱数
4 min read

Table of Content

はじめに

　Processing には randomGaussian という平均値 0、分散 1 の正規乱数を返す関数があります。この関数の戻り地は、その仕様の通り、非常に大きな値や絶対値が非常に大きなマイナスの値を返す可能性があります（もちろん、そのようなことが起きる確率は非常に小さい訳ですが）。

　多くの場合、この randomGaussian 関数で問題は無いかもしれません。しかし、もし範囲指定付きの正規乱数のようなものが存在すれば、それらの乱数を得た後に、ある範囲に乱数が収まっているか否かを調べるコードは不要となります。

このページでは、得られる乱数の範囲を指定できる疑似正規乱数を作る方法を説明します。

原理

　まず、正規乱数について少し説明をします。中心極限定理という、ものすごく難しそうな定理があります。これは、大雑把に言ってしまうと、一様乱数を生成する乱数発生装置から発生される乱数を、いくつか足し合わせて平均をとると、それらは正規乱数に近づくよ、というものです。

　統計的な話については、私の説明なんかよりも、以下のページが遥かに優れていると思います。なので、理論的な話に興味のある人は、こちらのページをご覧下さい：

「一様乱数の平均値を正規乱数として代用する」という話をゆるふわ統計的に検証する

　とにかく、一様乱数を足し合わせて平均をとれば良いのですが、これはサイコロを考えてみると分かりやすいと思います。

ちゃんとしたサイコロの出る目は、どの目も確率 1/6 です。これは 1 から 6 の整数を返す乱数とみると一様乱数となります。

次に、2 個のサイコロ振り、出た目の合計を考えるとどうでしょうか？ 1 つのサイコロの場合、1 から 6 の整数値が出る確率はどれも 1/6 で一様乱数でしたが、 2 個のサイコロの出た目の合計となると、圧倒的に 7 が出る確率が多くなります。

もちろんこれは、合計が 7 となる組み合わせの数が多くなるからに他ありません。一方、2 や 12 が出る組み合わせは、 2 つのサイコロの目が共に 1 であったり 6 であったりする場合に限ります。 2 つの目の和として 7 が出る場合と比べると、圧倒的にレアというか珍しい組み合わせのときのみです。

　一様乱数のグラフはどの目も確率 1/6 でしたので、単なる水平線になりますが、 2 個のサイコロの目の和の確率のグラフは、 7 のところが盛り上がったグラフとなります。

サイコロ 1 個の時と比べると、より正規分布に形が近づいていると言えます。

　とまあ、こんな感じで、一様乱数を足し合わせて、その平均を取ると正規乱数に近づいて行くことが分かるかと思います。

実際のコード

　使いやすいように、最小値 m、最大値 n の範囲の疑似正規乱数を返す関数 pseudoGaussian は以下のようになります：

def pseudoGaussian(m,n):
    t=0;N=10
    for i in [0]*N:
	    t+=random(n-m)/N
    return t+m

　プログラムの構造としては、範囲 0 から n-m の疑似正規乱数を生成し、最後に、得られた乱数に m を足しています。

このソースでは N 個の乱数を足し合わせて疑似正規乱数を生成しています。 N の値を大きくすればするほど、正規乱数に近づいていきますが、処理負荷は高くなります（なので、疑似正規乱数が得られるまで、より多くの時間がかかることとなります）。

　for ループの範囲指定のところは range(N) でも良いのですが、カウンタ変数 i を 0 から N-1 まで変化させる必要がなければ、このような書き方でも動きます（文字数が少ない書き方となりますので、これはつぶやき Processing むけの書き方です）。

実行例

　この範囲指定付きの疑似正規乱数を使ったテスト動画がありますので、以下に示します：

#つぶやきProcessing 範囲指定付き疑似 randomGaussian
def setup():size(500,500);blendMode(ADD);clear();fill(-1,10)

def N(s,e):
a=0;m=10
for i in [0]*m:a+=random(e-s)/m
return a+s

def draw():
for i in [0]*100:rect(N(50,450),N(50,450),2,2)

N(最小値,最大値) として使って下さい。 pic.twitter.com/LzQtleEJAA
— Koji Saito (@KojiSaito) May 20, 2020

　ソースコードは省略するため、関数名など変更されていますが、原理は全く同じものです。範囲指定の中央部分（=250=(50+450)/2）に多くの乱数が分布していることが確認できます。