Python には何通りかの乱数が用意されています。ここではシミュレーション用の乱数を説明します。パスワード生成のためには安全な乱数をご覧ください。
random は標準ライブラリのモジュールですのでインストールする必要はありません。Mersenne Twister が使われています。
import random random.random() # 引数をとらない。区間 [0,1) の1個の乱数を返す random.randrange(1, 7) # 1から6までの整数の乱数 random.randint(1, 6) # 1から6までの整数の乱数
本格的なシミュレーションで大量の乱数を必要とする場合には次の 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個
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個の配列
上の例のように「タネ(seed)」を与えるのがシミュレーションでの推奨です。同じタネからは同じ乱数列が生成されますので,あとで結果を再現できます。ゲームなどでは rng = np.random.default_rng()
のように空にすると,再現性のない(毎回異なる)乱数列が得られます。
乱数のタネとしては、その日の日付でも何でもいいのですが、たまたま自分の都合の良い結果が出るようなタネを採用したのではないかと疑われないように、自分のORCIDをタネに使うのも一つの手です。ORCID(Open Researcher and Contributor ID)は研究者用のIDで、https://orcid.org に登録すれば無料で発行してもらえます。私のORCIDは 0000-0003-4830-9421 ですので、頭の0と区切りのマイナスを取り除いた 348309421 を上の例で使っています。
例えばモンテカルロ法で円周率を求めるには次のようにします:
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.2.1、Python 3.11.7での速度比較です。一様乱数の配列を作って合計します。最初が標準ライブラリ、次が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)])
432 ms ± 1.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
import numpy as np rng = np.random.default_rng(348309421) %timeit np.sum(rng.random(N))
48.7 ms ± 12 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
import torch %timeit torch.rand(N).sum()
28.5 ms ± 23.7 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
GPU(mps)を指定してみます。
%timeit torch.rand(N, device="mps").sum()
The slowest run took 4.03 times longer than the fastest. This could mean that an intermediate result is being cached. 214 µs ± 129 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
キャッシュされたようなので速すぎます。次のほうが正しいようです。
%timeit torch.rand(N, device="mps").sum().item()
4.03 ms ± 107 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
import mlx.core as mx %timeit mx.sum(mx.random.uniform(shape=(10000,)))
5.95 µs ± 42.5 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
MLXは遅延評価のため、おかしな値になってしまいました。評価を強制してみます:
%timeit mx.sum(mx.random.uniform(shape=(N,))).item()
29.2 ms ± 18.5 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit mx.eval(mx.sum(mx.random.uniform(shape=(N,))))
29.1 ms ± 14.4 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
PyTorchでmpsを指定するよりかなり遅くなってしまいました。まだMLX 0.0.6なので将来に期待。