一様から何でも
乱数生成器は0から1まで一様に散らばった数を返します。モデルは別の形を必要とします。30%の確率で勝つコイン、約定間の待ち時間、正規分布の日次リターンなどです。やり方は同じで、一様な数に変換をかけます。最も一般的なのが逆変換法で、欲しい分布の累積分布関数を F とすると、一様な u に対する F⁻¹(u) はちょうどその分布に従います。配列全体を一度に変換できます。
import math
import numpy as np
# Exponential waiting times with rate 2 (a mean of 0.5): F(x) = 1 - exp(-rate * x), so x = -ln(1 - u) / rate
u = np.array([0.1, 0.5, 0.9])
print((-np.log(1 - u) / 2).round(4).tolist()) # [0.0527, 0.3466, 1.1513] small u, short wait
rng = np.random.default_rng(42)
waits = -np.log(1 - rng.random(100000)) / 2
print(round(waits.mean(), 2)) # 0.5 the mean is 1 / rate離散的な結果と正規分布
離散的な結果の表なら、rng.choice(values, size=n, p=probs) が逆変換を代わりに行い、各値をその確率で引きます。たとえば -2%、0%、+3% の日をそれぞれ20%、50%、30%の確率で引く場合です。正規分布では分布関数の逆関数に単純な式がないため、古典的な手法がボックス＝ミュラー法です。独立な一様乱数 u1、u2 から、z = sqrt(-2 ln u1) * cos(2 pi u2) が標準正規乱数になります。実際には rng.standard_normal(n) のほうが速いですが、一度書いてみれば魔法ではないとわかります。
rng = np.random.default_rng(3)
days = rng.choice(np.array([-0.02, 0.0, 0.03]), size=100000, p=[0.2, 0.5, 0.3])
print(round(days.mean(), 3), round(0.2 * -0.02 + 0.3 * 0.03, 3)) # 0.005 0.005 the sample mean matches the table
rng = np.random.default_rng(42)
u1, u2 = 1 - rng.random(100000), rng.random(100000) # 1 - random() is in (0, 1]: never log(0)
z = np.sqrt(-2 * np.log(u1)) * np.cos(2 * math.pi * u2)
print(abs(z.mean()) < 0.02, round(z.std(), 2)) # True 1.0 mean 0 and standard deviation 1標準正規乱数から価格へ
標準正規乱数 z はショックで、それが価格をどう動かすかはモデルが決めます。株の標準的なモデルである幾何ブラウン運動では、t 年後の価格は S = S0 * exp((mu - 0.5 sigma²) t + sigma sqrt(t) z) です。これは対数正規で、常に正であり、実際の価格のように右に長い裾を持ちます。-0.5 sigma² の項は平均を補正し、期待価格がちょうど mu で成長するようにします：E[S] = S0 * exp(mu t)。中央値は平均より低く、典型的な結果が平均的な結果より悪い理由です。
s0, mu, sigma, t = 100, 0.08, 0.2, 1
z = np.array([-1.0, 0.0, 1.0]) # a one-sigma loss, no shock, a one-sigma gain
prices = s0 * np.exp((mu - 0.5 * sigma ** 2) * t + sigma * math.sqrt(t) * z)
print(prices.round(2).tolist()) # [86.94, 106.18, 129.69] all positive, the upside is larger
draws = np.random.default_rng(0).standard_normal(100000)
simulated = s0 * np.exp((mu - 0.5 * sigma ** 2) * t + sigma * math.sqrt(t) * draws)
print(round(simulated.mean(), 1), round(s0 * math.exp(mu * t), 1)) # 108.3 108.3 the mean grows at mu
print(round(float(np.median(simulated)), 1)) # 106.1 the median is lower# NumPy ships ready-made samplers too; the transformations above show what they do inside
z = np.random.default_rng(0).standard_normal(100000)
print(abs(z.mean()) < 0.02, round(z.std(), 2)) # True 1.0 standard_normal: mean 0, standard deviation 1S0 * exp(mu t) です。符号の誤りや平方根の書き忘れは、もっともらしく見えて間違った数を生み、こうした確認だけがそれを明らかにします。- 逆変換法：一様乱数に分布関数の逆関数を適用する。指数分布は
-ln(1 - u) / rate rng.choice(values, size=n, p=probs)は離散表から抽出し、ボックス＝ミュラーsqrt(-2 ln u1) * cos(2 pi u2)は正規乱数を作るt年後の価格はS0 * exp((mu - 0.5 sigma²) t + sigma sqrt(t) z)（対数正規）。大標本は常に正確な平均と照らして確認する
次はモンテカルロ推定です。式のない平均や確率を、乱数の抽出で計算します。
ラボ：採点される内容
sample_discreteで結果の表から抽出する25点exponential_from_uniformで分布関数の逆関数を使う25点box_mullerで釣鐘型の乱数を作る25点terminal_pricesでショックを将来価格に変換する25点