ここではPythonを用いて任意の正規分布に従う乱数を生成する方法を説明します。Pythonを用いた乱数の生成は標準モジュールを用いた方法に加えて、numpyモジュール、scipyモジュールを用いた方法などがあります。ここではそれぞれの方法の使い方や違いなどについて説明していきます。
開発環境
- numpy 1.19.2
- python 3.7.9
- scipy 1.6.0
Contents
正規分布に従う乱数を1つ生成する方法
randomモジュールを用いる方法
標準ライブラリのrandomモジュールのgauss関数に正規分布の平均と標準偏差を指定すると、それに従うランダムな1つの値を取得することができます。例えば平均0、標準偏差1の標準正規分布に従うランダムな値を1つ取得してみましょう。
import random
rnd_val = random.gauss(0,1)
これで指定した正規分布に従うランダムな値を1つ取得することができます。
正規分布に従う乱数を大量に生成する方法
続いて正規分布に従うランダムな値を大量に取得する場合を見ていきましょう。実際にシミュレーションなどで使う場合はランダムな値を大量に用意して、1つ1つ計算していくような場合が多いです。そのような場合は先ほどのgauss関数をforループで繰り返すよりもNumPyやSciPyに用意されている関数を使った方が高速に処理を行うことができます。
numpyモジュールを用いる方法
NumPyでランダムな値をndarrayとして取得する関数はnumpy.randomモジュールのGeneratorクラスのnormalメソッドを用います。まずdefault_rng関数でGeneratorクラスのインスタンスを取得し、normalメソッドに平均や標準偏差、サンプル数を指定して実行します。なお、平均と標準偏差はデフォルトでは標準正規分布となっています。
それでは標準正規分布に従うランダムな値を10000個生成してみましょう。
import numpy as np
generator = np.random.default_rng()
rnd = generator.normal(size=10000)
print('サンプル数\t' + str(len(rnd)))
print('平均\t\t' + str(rnd.mean()))
print('標準偏差\t' + str(rnd.std()))
サンプル数 10000 平均 -0.006732355755997109 標準偏差 0.9979441628505185
これで平均0、標準偏差1の正規分布に従うランダムな値10000個が生成されました。
また、例えば平均10、標準偏差3の正規分布に従う乱数10000個は次のように生成できます。
generator = np.random.default_rng()
rnd = generator.normal(loc=10, scale=3, size=10000)
なお、正規分布に従うランダムな要素の10×10の二次元配列を取得したい場合は、size引数にsize=(10, 10)のように配列の形状を指定することもできます。
NumPy 1.17以前
なお、NumPy 1.17以前はnumpy.random.normal関数を用いていたので、現在でもこの書き方で説明される場合もあります。
import numpy as np
rnd = np.random.normal(size=10000)
print('サンプル数\t' + str(len(rnd)))
print('平均\t\t' + str(rnd.mean()))
print('標準偏差\t' + str(rnd.std()))
サンプル数 10000 平均 0.005879449052299143 標準偏差 0.9957912941811267
これでも同様の結果が得られるのですが、現在では以下のようにGeneratorクラスを用いる方法が推奨されているのでご注意ください。
scipyモジュールを用いる方法
SciPyを用いても同様に正規分布に従う乱数を取得することができます。正規分布を表すnormオブジェクトを用いて、そのrvsメソッドで乱数を生成することができます。こちらもデフォルトでは標準正規分布となっているので、分布の形を変える場合はloc引数に平均を、scale引数に標準偏差を指定しましょう。
それでは標準正規分布に従うランダムな値を10000個生成してみましょう。
from scipy.stats import norm
rnd = norm.rvs(size=10000)
print('サンプル数\t' + str(len(rnd)))
print('平均\t\t' + str(rnd.mean()))
print('標準偏差\t' + str(rnd.std()))
サンプル数 10000 平均 -0.007398366360305928 標準偏差 1.0040620356362158
これで平均0、標準偏差1の正規分布に従うランダムな値10000個が生成されました。
また、例えば平均10、標準偏差3の正規分布に従う乱数10000個は次のように生成できます。
rnd = norm.rvs(loc=10, scale=3, size=10000)
それぞれの方法の比較
10000個程度の乱数の取得では上記のどの方法を用いても実行速度に大差はありませんが、乱数のサイズが大きくなると差が出てきます。それでは具体的にどのくらいの差があるのかを比較してみましょう。
NumPyを用いる方法
まずは、NumPyを用いて標準正規分布に従う乱数1,000,000,000個を生成してみます。
import time
import numpy as np
t_start = time.time()
generator = np.random.default_rng()
rnd = generator.normal(size=1000000000)
t_stop = time.time()
print(t_stop - t_start)
10.88458514213562
約10.9秒で1,000,000,000個の乱数を生成できました。
NumPyを用いる方法 (NumPy 1.17 以前)
続いて、NumPy 1.17以前の方法でも同様に計測してみましょう。
import time
import numpy as np
t_start = time.time()
rnd = np.random.normal(size=1000000000)
t_stop = time.time()
print(t_stop - t_start)
19.3181574344635
こちらは約19.3秒かかっており、現在推奨されているGeneratorクラスを用いる方法の方が高速であると言えます。
SciPyを用いる方法
最後にSciPyを使って同様に1,000,000,000個の乱数を生成してみましょう。
import time
from scipy.stats import norm
t_start = time.time()
rnd = norm.rvs(size=1000000000)
t_stop = time.time()
print(t_stop - t_start)
22.874540328979492
この方法では約22.9秒かかってしまい、最も時間のかかる方法となってしまいました。
以上から一番高速に乱数を生成できる方法は現在推奨されているNumPyのGeneratorクラスを用いる方法だと言えます。
これは従来の方法では乱数生成にメルセンヌ・ツイスタ(MT19937)という生成器を用いていたのですが、新しい方法ではPCG64という新しくてより高速な方法に更新されたことによります。なお、GeneratorクラスではPCG64以外の乱数生成器も選択することが可能となっているので、従来からのメルセンヌ・ツイスタを選択することも可能となっています。
SciPyを用いる方法の乱数生成器はデフォルトでメルセンヌ・ツイスタ(MT19937)なので、NumPy1.17以前と同様に時間がかかってしまいます(実際には他の処理も加わるせいかさらに遅くなっています)。SciPyを用いてPCG64で乱数生成するためには次のようにrandom_state引数に指定してください。
import time
from scipy.stats import norm
import numpy as np
t_start = time.time()
generator = np.random.default_rng()
rnd = norm.rvs(size=1000000000, random_state=generator)
t_stop = time.time()
print(t_stop - t_start)
12.016556978225708
ただし、他の処理も加わるせいか純粋にNumPyで乱数生成するよりも若干時間がかかってしまっています。速度を重視するならNumPyで乱数生成した方がよさそうです。