🚀

Box-Muller変換の証明とシュミレーション

2024/08/09に公開

Box-Muller変換

はじめに以下の確率変数を定義する.

\begin{align*} & U_{1}, U_{2}, \text { i.i.d. } \sim U(0,1) \end{align*}

以下の確率変数変換を考える．

\begin{align*} X &= R\cos(\Theta) = (-2\log U_1) \cos(2\pi U_2)\\ Y &= R\sin(\Theta) = (-2\log U_1) \sin(2\pi U_2) \end{align*}

このとき，以下が成り立つ．

\begin{align*} X, Y, \text { i.i.d. } \sim N(0,1) \end{align*}

証明

はじめに以下の確率変数を定義する.

\begin{align*} & U_{1}, U_{2}, \text { i.i.d. } \sim U(0,1) \end{align*}

この確率変数について，以下の変数変換を考える.

\begin{align*} \left\{ \begin{array}{l} R = (-2\log U_1)^{-\frac{1}{2}} \\ \theta = 2\pi U_2 \end{array} \right . \end{align*}

逆変換は以下である．

\begin{align*} \left\{ \begin{array}{l} U_1 = \exp({-\frac{R^2}{2}}) \\ U_2 = \frac{\theta}{2\pi} \end{array} \right . \end{align*}

ここで，ヤコビアンは以下のように計算される.

\begin{align*} J((r, \theta) \rightarrow (u_1, u_2)) &=\operatorname{det}\left(\begin{array}{cc} - r e^{-\frac{r^{2}}{2}} & 0 \\ 0 & (2 \pi)^{-1} \end{array}\right) \\ & = -(2 \pi)^{-1} r e^{-\frac{r^{2}}{2}} \cdot 1 \\ & = -(2 \pi)^{-1} r e^{-\frac{r^{2}}{2}} \\ \end{align*}

$U(0,1)$ の密度が範囲内で $1$ であることに注意して変形を行う．

\begin{align*} f_{R \theta}(r, \theta) &=f_{U_{1} U_{2}}\left(u_{1}, u_{2}\right)\left|J\left(\left(r, \theta\right) \rightarrow(u_1, u_2)\right)\right| \\ & =\frac{r e^{-\frac{r^{2}}{2}}}{2 \pi} \\ \end{align*}

これで， $R$ , $\Theta$ の同時分布が求まった．
次に，曲座標変換を用いる.

\begin{align*} & \left\{ \begin{array}{l} X=R \cos \theta \\ Y=R \sin \theta \end{array} \right . \end{align*}

ここで， $R$ , $\Theta$ の逆変換を求めるのは容易ではない．そこで，以下の公式を使う.

J((x, y) \rightarrow(r, \theta)) = \frac{1}{J((r, \theta) \rightarrow(x, y))}

上式を用いると，ヤコビアンを計算できる．

\begin{align*} J((x, y) \rightarrow(r, \theta)) &= J((r,\theta) \rightarrow (x,y))^{-1} \\ &= \operatorname{det}\left(\begin{array}{lr} \cos \theta & -r \sin \theta \\ \sin \theta & r \cos \theta \end{array}\right)^{-1} \\ &=r^{-1} \end{align*}

よって求めるX, Yの分布は以下のように定まる．

\begin{align*} f_{X Y}(x, y) &=f_{R \theta}(r, \theta)|J((x, y) \rightarrow(r, \theta))| \\ & =\frac{r e^{-\frac{r^{2}}{2}}}{2 \pi} \cdot r^{-1} \\ & =\frac{1}{2 \pi} e^{-\frac{x^{2}+y^{2}}{2}} \\ & =\frac{1}{\sqrt{2 \pi}} e^{-\frac{x^{2}}{2}} \cdot \frac{1}{\sqrt{2 \pi}} e^{-\frac{y^{2}}{2}} \\ &= f_X(x) f_Y(y) \end{align*}

したがって、 $X$ と $Y$ はそれぞれ独立に標準正規分布に従う．

シュミレーション

実際にこの変換がうまくいくのか検証を行う．

プログラムを作成する

描画する部分以外は標準関数とmath mouduleのみで動く．実装に際して，generate_samples(n: int)でn // 2を行なっているが，これは，Box-Muller変換が一度に二つの確率変数を生成するためである．

from typing import Tuple

import random
import math
import matplotlib.pyplot as plt
import seaborn as sns


def generate_uniform_random_pair() -> Tuple[float, float]:
    """ 独立に一様分布に従う2つの乱数を生成する。 """
    u1 = random.random()
    u2 = random.random()
    return u1, u2

def box_muller_transform(u1: float, u2: float) -> Tuple[float, float]:
    """ Box-Muller 変換を使用して、独立に一様分布に従う2つの乱数を正規分布に確率変数変換する。 """
    z0 = math.sqrt(-2.0 * math.log(u1)) * math.cos(2.0 * math.pi * u2)
    z1 = math.sqrt(-2.0 * math.log(u1)) * math.sin(2.0 * math.pi * u2)
    return z0, z1

def generate_samples(n: int):
    """ n 個の正規分布に変換された乱数を生成する。 """
    samples = []
    for _ in range(n // 2):
        u1, u2 = generate_uniform_random_pair()
        z0, z1 = box_muller_transform(u1, u2)
        samples.append(z0)
        samples.append(z1)
    return samples[:n]

# Seaborn and Matplotlibの準備
sns.set()
sns.set_style("whitegrid", {'grid.linestyle': '--'})
sns.set_context("talk", 0.8, {"lines.linewidth": 0})
colors = sns.color_palette("cool", 5)

# サンプルサイズを指定
sample_sizes = [10, 50, 100, 1000, 10000]

# それぞれのサイズでシュミレーションを行い, 結果を描画する
for i, n in enumerate(sample_sizes):
    samples = generate_samples(n)
    plt.figure(figsize=(10, 6))  # Fixed size for all figures
    sns.histplot(samples, kde=True, color=colors[i], binwidth=0.5)
    plt.title(f'Box-Muller Transformed Samples for n={n}')
    plt.xlabel('Value')
    plt.ylabel('Frequency')
    plt.show()