連続一様分布に従う浮動小数点数のサンプルが従う離散分布

概要

(連続) 一様分布に従う浮動小数点数のサンプルは (離散) 一様分布には従わないという話.

前提知識

浮動小数点数の大まかな仕様と基礎的な確率論の知識を前提とする.

解説

コンピュータでは連続濃度の実数を扱いきれないため浮動小数点数による離散集合で近似され, 連続分布に従うサンプルも実際には浮動小数点数の集合上の離散分布に従うサンプルで実現される.
ここでは連続一様分布を例に, それに従う浮動小数点数のサンプルが具体的にどのような離散分布に従っているのか見てみる.

浮動小数点数の記述

$2$ 以上の整数 $b$ と正整数 $r$ を固定した時, 浮動小数点数 (のうち正規化数) の一般形は整数 $s, e, k$ (それぞれ符号部, 指数部, 仮数部と呼ばれる) に対して次で与えられる:

$$
f(s, e, k) := (-1)^s\left(1+(b-1)\frac{k}{r}\right)b^e.
$$

一般には $s = 0, 1$ 及び $0 \leq k \leq r-1$ という制限が与えられ, その制限の下では $s$ で符号を決めて $e$ によって実数全体を区間 $[b^e, b^{e+1})$ または $(-b^{e+1}, -b^e]$ に区切り, 各区間を $r$ 等分して $k$ でインデクス付ける表現方法と説明できる.
以降は特に言及がなければこの制限を仮定して議論を進めるが, 後に制限から外れた記述も行うので覚えておこう.

このように定まる浮動小数点数は (前述の制限の下では) 数値のスケールに比例して隣り合った数値の間隔が大きくなるという特性を得られるが $0$ を扱えないという欠点があり, さらにコンピュータで効率的に扱うため $e$ の範囲を制限すると, $0$ 周辺で解像度が大きくなってしまうという現象が発生する.
実際, $m \leq e \leq M$ 及び $k \geq 0$ に制限すると正の最小値 $f(0, m, 0)$ とその前後の数は,

$$
\begin{align*}
f(1, m, 0) & = -b^m, \\
f(0, m, 0) & = b^m, \\
f(0, m, 1) & = \left(1+\frac{b-1}{r}\right)b^m, \\
\end{align*}
$$

であり, 正の最小値の両隣の数との差はそれぞれ,

$$
\begin{align*}
f(0, m, 0)-f(1, m, 0) & = 2b^m, \\
f(0, m, 1)-f(0, m, 0) & = \frac{b-1}{r}b^m, \\
\end{align*}
$$

となって比率は $2r:b-1$ となる.
標準的には (IEEE 754等) $b$ は高々10の小さい値である一方で $r$ は1000を超える非常に大きい値を取るため, 0周辺に全く値を取れない区間が存在してしまうことが分かる.

そのため $e$ の取り得る値として $m-1$ を許容し,

$$
f(s, m-1, k) := (-1)^s\frac{k}{r}b^m,
$$

によって区間 $[-b^m, b^m]$ を $2r$ 等分する形式と組み合わせて浮動小数点数を構成する方法が標準的に取られる¹.
この時, $e \geq m$ で表される数を正規化数, $e = m-1$ で表される数を非正規化数と呼ぶ.
非正規化数は $0$ を挟んだ区間を等分するという点で固定小数点数という形式と等価である.

浮動小数点数の集合

浮動小数点数を与える関数 $f(s, e, k)$ は, $s = 0, 1$ 及び $0 \leq k \leq r-1$ の制限の下では $e = m-1, k = 0$ の時 $s = 0, 1$ によらず $0$ を取る以外は単射である (演習).
つまり,

$$
f(s, e, k) = f(s', e', k') \neq 0 \ \Longrightarrow \ (s, e, k) = (s', e', k'),
$$

であり, $f_{b, r}(0, m-1, 0) = +0, f_{b, r}(1, m-1, 0) = -0$ 等と書いて $\pm0$ を区別することにすれば $(s, e, k)$ の組み合わせによって重複なく数を表現できる形式となる.

これにより浮動小数点数の集合は,

$$
\mathcal{F} := \left\{f(s, e, k) \in \mathbb{R} \ \middle\vert \ \begin{array}{l} s = 0, 1, \\ m-1 \leq e \leq M, \\ 0 \leq k \leq r-1 \end{array}\right\},
$$

と表され, 前述のように $\pm0$ を区別すれば,

$$
\begin{align*}
\mathcal{F} & = \left(\bigsqcup_{e=m-1}^M\mathcal{F}_e^+\right)\sqcup\left(\bigsqcup _{e=m-1}^M\mathcal{F}_e^-\right), \\
\mathcal{F}_e^+ & := \{f(0, e, k) \ \vert \ 0 \leq k \leq r-1\}, \\
\mathcal{F}_e^- & := \{f(1, e, k) \ \vert \ 0 \leq k \leq r-1\}, \\
\end{align*}
$$

と分解される.

浮動小数点数による連続分布の実現

現代のコンピュータにおける標準的な方法では離散的なデータしか扱えないため, 連続濃度集合である実数集合は浮動小数点数による離散濃度集合で近似される.
よって, 実数上の連続分布に従うサンプルも浮動小数点数上の離散分布に従うサンプルで近似される.

具体的には, 実数の部分集合 $\Omega \subset \mathbb{R}$ と $\Omega$ 上の連続分布 $\mu$ があった時, $X \subset \Omega$ を近似する浮動小数点数の集合 $\mathcal{F}_X := X\cap\mathcal{F}$ に対して,

$$
\mu(X) \sim \mu_\mathcal{F}(\mathcal{F}_X) := \sum _{x \in \mathcal{F}_X}p_x,
$$

という近似が成り立つ $\mathcal{F}_\Omega$ 上の離散分布 $\mu _\mathcal{F}$ を考える².
近似 $\sim$ は $|\mu(X)-\mu _\mathcal{F}(\mathcal{F}_X)|$ が常に小さい程度の意味付けを与えておくと, これにより $\mu _\mathcal{F}$ に従うサンプルによって $\mu$ に従うサンプルの近似が得られる.

$x \in \mathcal{F}_\Omega$ を取る確率質量 $p_x$ は, 適当な $\Omega$ から $\mathcal{F} _\Omega$ への丸め関数 $\rho: \Omega \rightarrow \mathcal{F} _\Omega$ を用いて $p_x := \mu(\rho^{-1}(x))$ で与えればいい.
丸め関数については特に定義はないが,

$$
\left\{\begin{align*}
x \in \mathcal{F} & \ \Longrightarrow \ \rho(x) = x, \\
x \leq y & \ \Longrightarrow \ \rho(x) \leq \rho(y), \\
\end{align*}\right.
$$

という「丸め」に自然な条件を満たしていれば常に $p_x$ が計算可能になる³.
四捨五入や切り捨て, 切り上げといった一般的な丸め関数は全てこれらの条件を満たしている.

連続一様分布を近似する浮動小数点数集合上の離散分布

一般に連続一様分布は, 半開区間 $\Omega = [x_{\min}, x_{\max})$ 上で $\mathcal{U}([x, y)) = \frac{y-x}{x_{\max}-x_{\min}}$ から生成される確率分布として定義される.
開区間, 閉区間を用いてもほぼ同じ分布が得られるが後の議論を簡単にするためここでは半開区間を使って書き, さらに簡単のため浮動小数点数のパラメータ $e$ の範囲を定める $m, M$ について $m \leq -1 \leq M$ と条件付けて, $x_{\min} = 0, x_{\max} = 1$ に固定すると,

$$
\mathcal{F}_{[0, 1)} = [0, 1)\cap\mathcal{F} = \bigsqcup _{e=m-1}^{-1}\mathcal{F}_e^+,
$$

の上の離散確率分布 $\mathcal{U}_\mathcal{F}$ で $\mathcal{U}$ を近似することになる.
これにより以降は $s = 0$ に制限して話を進めるため, $f^+(e, k) := f(0, e, k)$ という記号を導入しておこう.

この時, 丸め関数として,

$$
\rho(x) = \mathrm{floor}_\mathcal{F}(x) := \max{\{y \in \mathcal{F} \ \vert \ y \leq x \}},
$$

で定まる実数を浮動小数点数に切り下げる関数を取るとこれは前節の条件を満たし, $f^+(e, r) = f^+(e+1, 0)$ に注意すれば,

$$
\rho^{-1}(f^+(e, k)) = [f^+(e, k), f^+(e, k+1)),
$$

より,

$$
p_{f^+(e, k)} = \mathcal{U}(\rho^{-1}(f^+(e, k))) = \begin{cases} (b-1)\frac{1}{r}b^e, & m \leq e \leq -1, \\ \frac{1}{r}b^m, & e = m-1, \end{cases}
$$

が得られる⁴.
これは $e$ によって $f^+(e, k)$ を取り得る確率 $p_{f^+(e, k)}$ が変わるため, 離散一様分布ではない.

逆に $\mathcal{U}_\mathcal{F}$ を離散一様分布とした場合, ある $\alpha$ ($0 \lt \alpha \lt 1$) があって任意の $x \in \mathcal{F} _{[0, 1)}$ に対して $p_x = \alpha$ となり, $[b^e, b^{e+1}) \subset [0, 1)$ ($m \leq e \leq -1$) に対して,

$$
\begin{align*}
& \hphantom{=} \ \ \mathcal{U}([b^e, b^{e+1})) = b^{e+1}-b^e \\
& \sim \mathcal{U}_{\mathcal{F}}(\mathcal{F} _{[b^e, b^{e+1})}) = \sum _{x \in \mathcal{F}_e^+}p _{f^+(e, k)} = r\alpha, \\
\end{align*}
$$

と近似することになる.
これは, $b^{e+1}-b^e$ が $e$ に依存して単調に増加する一方でその近似値である $r\alpha$ は一定であり, 良い近似ではないことが分かる.

ただし, 先に書き下したように丸め関数を $\mathrm{floor}_\mathcal{F}(x)$ とした時の $\mathcal{U} _\mathcal{F}$ の確率質量は, $e$ には依存するが $k$ には依らない値となる.
これは半開区間 $[b^e, b^{e+1})$ (及び $[0, b^m)$) 上の連続一様分布であれば $\mathcal{F} _{[b^e, b^{e+1})} = \mathcal{F}_e^+$ 上の離散一様分布でうまく近似できることを示している⁵.

整数型一様乱数に基づく連続一様乱数 (の近似) の生成

前節の丸め関数を $\rho(x) = \mathrm{floor}_\mathcal{F}(x)$ とした時の離散分布 $\mathcal{U} _\mathcal{F}$ は一様ではないが, 各点を取る確率 $p_x$ が有理数となるため整数上の離散一様分布に従う乱数を変換して生成することが出来る.

$x \in \mathcal{F}_{[0, 1)}$ に対して $p_x$ が取り得る値は $\epsilon = \frac{1}{r}b^m$ の整数倍となっており,

$$
p_{f^+(e, k)} = \begin{cases} (b-1)\frac{1}{r}b^e = (b-1)b^{e-m}\epsilon, & m \leq e \leq -1, \\ \frac{1}{r}b^m = \epsilon, & e = m-1, \end{cases}
$$

と書ける.
この $\epsilon$ の係数 $(b-1)b^{e-m}$ 及び $1$ の $f^+(e, k) \in \mathcal{F}_{[0, 1)}$ における総和は,

$$
N := \sum _{x \in \mathcal{F} _{[0, 1)}}\frac{p_x}{\epsilon} = \frac{1}{\epsilon} = rb^{-m},
$$

で得られて, $\{0, 1, \dots, N-1\}$ 上の離散一様分布に従うサンプルを $p_{f^+(e, k)}$ に従う比率で $\mathcal{F}_{[0, 1)}$ に変換すれば $\mathcal{U} _\mathcal{F}$ に従うサンプルが生成できるということになる.

この変換 $\varphi: \{0, 1, \dots, N-1\} \rightarrow \mathcal{F}_{[0, 1)}$ としては例えば,

$$
\varphi(z) = \mathrm{floor}_\mathcal{F}(z\epsilon),
$$

が考えられ,

$$
\begin{align*}
\varphi^{-1}(f^+(e, k)) & = \{z \ \vert \ f^+(e, k) \leq z\epsilon \lt f^+(e, k+1)\} \\
& = \{z \ \vert \ Nf^+(e, k) \leq z \lt Nf^+(e, k+1)\}, \\
\end{align*}
$$

の要素数 $N(f^+(e, k+1)-f^+(e, k))$ は $(b-1)b^{e-m}$ 及び $1$ に一致する.
この時, $z \in \{0, 1, \dots, N-1\}$ に対応する浮動小数点数のパラメータ $e_z, k_z$ は,

$$
f^+(e_z, k_z) \leq z\epsilon \lt f^+(e_z, k_z+1),
$$

を解いて,

$$
(e_z, k_z) = \begin{cases} \left(m+\mathrm{int}\left(\log_b{\frac{z}{r}}\right), \left(zb^{-\mathrm{int}\left(\log_b{\frac{z}{r}}\right)}-r\right)//(b-1)\right), & z \geq r, \\ (m-1, z), & z \leq r-1, \end{cases}
$$

で得られる.
ただし $x//y$ は $x$ を $y$ で割り切った整数値で $\mathrm{int}(x)$ は $x$ を切り捨てた整数値, つまり,

$$
\begin{gather*}
x//y := \mathrm{floor}\left(\frac{x}{y}\right) = \max{\{y \in \mathbb{Z} \ \vert \ y \leq x \}}, \\
\mathrm{int}(x) := x//1,
\end{gather*}
$$

である.

実験

1. 整数乱数から浮動小数点数パラメータを経て浮動小数点数乱数を生成

単純に整数乱数を前節の $\varphi$ で変換する方法ではよく知られた一様乱数生成法が得られるだけ⁶なのでパラメータ $e, k$ を計算して byte 列から浮動小数点数に変換する方法を試してみる.
Python の float 型は標準的な 64bit (8byte) の浮動小数点数で, 大域的なパラメータ $b, r, m$ 及び $N$ は,

$$
\begin{cases}
b = 2, \\
r = 2^{52} = 4503599627370496, \\
m = -2^{10}+2 = -1022, \\
N = 4503599627370496\times2^{1022} = 2^{1074} \\
\hphantom{N} = 202402253307310618352495346718917307049556649764142118356901358027430339567995346891960383701437124495187077864316811911389808737385793476867013399940738509921517424276566361364466907742093216341239767678472745068562007483424692698618103355649159556340810056512358769552333414615230502532186327508646006263307707741093494784, \\
\end{cases}
$$

である.
$N$ の値がおかしなことになっているが, Python の整数は多倍長なので問題なく扱える.
浮動小数点数と byte 列の相互変換はこの辺りを参照.

import struct, random

# 符号部, 指数部, 仮数部から浮動小数点数の構築
# bit shift と加算で結合して bytes 化して float 化
def construct_float(sgn, exp, frc):
    b = ((sgn<<63)+(exp<<52)+frc).to_bytes(8, 'big')
    return struct.unpack('>d', b)[0]

# 浮動小数点数の大域的パラメータ設定
b = 2
r = 2**52
m = -2**10+2
N = r*b**(-m)

# 整数型一様乱数 random.randrange による連続一様乱数生成
# 整数乱数生成して e_z, k_z 計算して float 構築
def random_uniform():
    z = random.randrange(N)
    if z<r:
        e_z = m-1
        k_z = z
    else:
        zr = z//r
        ilog = zr.bit_length()-1
        e_z = m+ilog
        k_z = (z>>ilog)-r
    return construct_float(0, e_z-(m-1), k_z)

$10^4$ 個ほどサンプル生成してヒストグラムを描いてみると一様に分布しているらしいことは見て取れる.

import matplotlib.pyplot as plt

n = 10**4
x = []
for i in range(n):
    x.append(random_uniform())

plt.hist(x, bins=50, density=True)
plt.show()

なお, 遅い.

%timeit -n 100 -r 100 random_uniform()
%timeit -n 100 -r 100 random.random()

4.77 µs ± 935 ns per loop (mean ± std. dev. of 100 runs, 100 loops each)
68.6 ns ± 21.3 ns per loop (mean ± std. dev. of 100 runs, 100 loops each)

2. 浮動小数点数乱数のパラメータの分布を確認

今度は一般的な方法 (random.random()) で生成した乱数を浮動小数点数パラメータに変換してその分布を見てみる.
まず random.random() で生成した浮動小数点数乱数を byte 列に変換して整数を経て bit 列に変換する.
ここは回りくどいけどショートカットする方法が分からなかった.
符号部が常に $0$ だという前提で得られた bit 列から指数部と仮数部を切り出して再び整数に戻す.

プロットすると確率質量が指数部に依存し, 仮数部に関しては一定になっていることが見て取れる.

n = 10**4
e = []
f = []
for i in range(n):
    z = random.random()
    b = bin(int.from_bytes(struct.pack('>d', z), 'big'))
    e.append(int(b[2:-52], 2))
    f.append(int(b[-52:], 2))

plt.hist(e, bins=50, density=True)
plt.show()

plt.hist(f, bins=50, density=True)
plt.show()

参考資料

754-2008 - IEEE Standard for Floating-Point Arithmetic | IEEE Standard | IEEE Xplore

IEEE 754 - Wikipedia

1. 最大値以上の区間についても $e = M+1$ を許容して無限大や nan 等を表現する方法も標準的に採用されるがここでは割愛する. ↩

2. 測度論的には $X \subset \Omega$ として可測集合に限定する必要があるがここでは無視しておく. ↩

3. つまりこの丸め関数は $\Omega \subset \mathbb{R}$ 上の標準的な可測空間において可測関数となる. ↩

4. この他の IEEE 754 等に基づく丸め関数から得られる離散分布の書き下しは読者に任せる. ↩

5. 取得可能な値は $[0, 1)$ を近似した場合の $-\frac{1}{m}$ 通りとなり, サンプルの解像度は大きく落ちる. ↩

6. 浮動小数点数の標準的な丸めが切り捨てでなく偶数への最近接丸めなので追加の処理は必要になる. ↩