乱数

Python には何通りかの乱数が用意されています。ここではシミュレーション用の乱数を説明します。パスワード生成のためには安全な乱数をご覧ください。

標準ライブラリ random

random は標準ライブラリのモジュールですのでインストールする必要はありません。Mersenne Twister が使われています。

import random

random.random()  # 引数をとらない。区間 [0,1) の1個の乱数を返す
random.randrange(1, 7)  # 1から6までの整数の乱数
random.randint(1, 6)    # 1から6までの整数の乱数

次の NumPy による方法と違って、整数の乱数の範囲に制限はありません。

本格的なシミュレーションで大量の乱数を必要とする場合には次の NumPy による方法を使います。一方で、非常に大きい整数の乱数を生成するには、NumPy は不適です。

ライブラリ NumPy

numpy.random を使います。

古い方法

Mersenne Twister が使われています。

import numpy as np

np.random.random(10)             # 0以上1未満の一様乱数10個
np.random.randint(1, 7, size=10) # 1以上7未満の整数の乱数10個

np.random.randint() は64ビット整数の範囲でしか使えません。

新しい方法

2013年にリリースされた NumPy 1.17.0 からは Generator を使った新しいインターフェースが推奨されています（NumPy のバージョンは np.__version__ でわかります）。デフォルトのアルゴリズムは PCG64 です。これは次のようにして使います：

import numpy as np

rng = np.random.default_rng(348309421) # seedの設定（空も可）

rng.random()             # 0以上1未満の一様乱数
rng.random(10)           # 0以上1未満の一様乱数10個の配列
rng.normal(5, 3)         # 平均5，標準偏差3の正規分布の乱数
rng.normal(5, 3, 10)     # 平均5，標準偏差3の正規分布の乱数10個の配列
rng.standard_normal()    # 標準正規分布（平均0，標準偏差1）の乱数
rng.standard_normal(10)  # 標準正規分布の乱数10個の配列
rng.integers(1, 7)       # 1以上7未満の整数の乱数
rng.integers(1, 7, 10)   # 1以上7未満の整数の乱数10個の配列
rng.exponential(3)       # 平均3の指数分布の乱数
rng.exponential(3, 10)   # 平均3の指数分布の乱数10個の配列

rng.integers() は64ビット整数の範囲でしか使えません。

上の例のように「タネ（seed）」を与えるのがシミュレーションでの推奨です。同じタネからは同じ乱数列が生成されますので，あとで結果を再現できます。ゲームなどでは rng = np.random.default_rng() のように空にすると，再現性のない（毎回異なる）乱数列が得られます。

例えばモンテカルロ法で円周率を求めるには次のようにします：

import matplotlib.pyplot as plt

N = 10000
x = rng.random(N)
y = rng.random(N)
c = x*x + y*y < 1
plt.scatter(x, y, c=c)
plt.axis('scaled')
print('pi =', 4 * np.sum(c) / N)

ランダムウォークの例です：

N = 1000
x = np.cumsum(rng.random(N) - 0.5)
y = np.cumsum(rng.random(N) - 0.5)
plt.plot(x, y)
plt.axis('equal')

硬貨投げのシミュレーションは rng.integers(0, 2, size=100) のようにできます。表の枚数は np.sum(rng.integers(0, 2, size=100)) です。これを何回もやってヒストグラムを描いてみましょう。

a = [np.sum(rng.integers(0, 2, size=100)) for _ in range(10000)]

plt.hist(a, bins=range(min(a), max(a)+2), align="left", edgecolor="black")

plt.hist() のオプションがわずらわしいなら seaborn を使いましょう：

import seaborn as sns

sns.histplot(a, discrete=True)

速度比較

Mac mini (M1, 2020) 16GB Sonoma 14.5、Python 3.12.4、NumPy 2.0.0、PyTorch 2.3.1、MLX 0.15.1での速度比較です。一様乱数の配列を作って合計します。最初が標準ライブラリ、次がNumPy、PyTorch、最後がApple Silicon専用の MLX です。PyTorchとMLXは32ビット精度、他は64ビット精度です。

[2023-12-23] 最初 N=10000 でやっていましたが、N を増やして再計算しました。@mosko_mule さん、@BOBtakumi さんのご教示に感謝します。

import random

N = 10_000_000

%timeit sum([random.random() for _ in range(N)])

495 ms ± 485 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Numba 0.60.0 で高速化してみましょう。

from numba import jit

@jit
def foo():
    return sum([random.random() for _ in range(N)])

%timeit foo()

48.4 ms ± 122 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

NumPy です。

import numpy as np

%timeit np.sum(np.random.random(N))

40.1 ms ± 149 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

rng = np.random.default_rng(348309421)

%timeit np.sum(rng.random(N))

48.8 ms ± 12.6 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

PCG64よりMersenne Twisterのほうが若干速いようです。

PyTorch です。

import torch

%timeit torch.rand(N).sum()

31.4 ms ± 339 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

GPU（mps）を指定してみます。torch.rand(N, device="mps").sum() ではキャッシュされるようで、次のようにするのが正しいようです。

%timeit torch.rand(N, device="mps").sum().item()

2.96 ms ± 14.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

MLX です。

import mlx.core as mx

MLXは遅延評価のため、mx.sum(mx.random.uniform(shape=(10000,))) では速すぎる値になってしまいました。評価を強制してみます：

%timeit mx.sum(mx.random.uniform(shape=(N,))).item()

29.2 ms ± 41.5 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit mx.eval(mx.sum(mx.random.uniform(shape=(N,))))

29.2 ms ± 22.1 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

PyTorchでmpsを指定するよりかなり遅くなってしまいました。将来に期待。