円周上・円周内・球面上・球面内の一様乱数

円周上・円周内の一様乱数および球面上・球面内の一様乱数について記述します。簡略化のために円と球の半径は$1$としています。

円周上の一様乱数

半径$1$の円周上の点は次のように表すことができます。

\left\{
\begin{array}{ll}
x = \cos\theta \\
y = \sin\theta
\end{array}
\right. \\
0 \leqq \theta \leqq 2\pi

円周上に一様に分布させるには$[0, 1]$の一様乱数$u$を用いて$\theta$を次のように求めます。

\theta = 2\pi u

円周内の一様乱数

半径$1$の円周内の点は次のように表すことができます。

\left\{
\begin{array}{ll}
x = r\cos\theta \\
y = r\sin\theta
\end{array}
\right. \\
0 \leqq \theta \leqq 2\pi, 0 \leqq r \leqq 1

円周内に一様に分布させるには$[0, 1]$の一様乱数$u, v$を用いて$r, \theta$を次のように求めます。

\theta = 2\pi u \\
r = \sqrt{v} 

球面上の一様乱数

半径$1$の球面上の点は次のように表すことができます。

\left\{
\begin{array}{lll}
x = \sin\theta \cos\phi \\
y = \sin\theta \sin\phi \\
z = \cos\theta
\end{array}
\right. \\
0 \leqq \theta \leqq \pi, 0 \leqq \phi \leqq 2\pi

球面上に一様に分布させるには$[0, 1]$の一様乱数$u, v$を用いて$\cos\theta, \phi$を次のように求めます。

\cos\theta = -2u + 1 \\
\phi = 2\pi v \\

$\sin^2\theta+\cos^2\theta=1$かつ$0 \leqq \theta \leqq \pi$より$\sin\theta=\sqrt{1 - \cos^2\theta}$となるので、球面上に一様に分布する点は$u, v$を用いて次のように表すことができます。

\left\{
\begin{array}{lll}
x = \sqrt{1 - z^2} \cos2\pi v \\
y = \sqrt{1 - z^2} \sin2\pi v \\
z = -2u + 1
\end{array}
\right. \\

球面内の一様乱数

半径$1$の球面内の点は次のように表すことができます。

\left\{
\begin{array}{lll}
x = r\sin\theta \cos\phi \\
y = r\sin\theta \sin\phi \\
z = r\cos\theta
\end{array}
\right. \\
0 \leqq \theta \leqq \pi, 0 \leqq \phi \leqq 2\pi, 0 \leqq r \leqq 1

球内に一様に分布させるには$[0, 1]$の一様乱数$u, v, w$を用いて$\cos\theta, \phi, r$を次のように求めます。

\cos\theta = -2u + 1 \\
\phi = 2\pi v \\
r = \sqrt[3]{w}

よって球内に一様に分布する点は$u, v, w$を用いて次のように表すことができます。

\left\{
\begin{array}{lll}
x = \sqrt[3]{w}\sqrt{1 - z^2} \cos2\pi v \\
y = \sqrt[3]{w}\sqrt{1 - z^2} \sin2\pi v \\
z = \sqrt[3]{w}(-2u + 1)
\end{array}
\right. \\

実装

JavaScriptで実装してみます。サンプルイメージはp5.jsで生成しています。

円周上の一様乱数

function randomOnCircle() {
  const theta = 2.0 * Math.PI * Math.random();
  return [Math.cos(theta), Math.sin(theta)];
}

random on circle with p5.js

円周内の一様乱数

function randomInCircle() {
  const theta = 2.0 * Math.PI * Math.random();
  const radius = Math.sqrt(Math.random());
  return [radius * Math.cos(theta), radius * Math.sin(theta)];
}

random in circle with p5.js

球面上の一様乱数

function randomOnSphere() {
  const cosTheta = -2.0 * Math.random() + 1.0;
  const sinTheta = Math.sqrt(1.0 - cosTheta * cosTheta);
  const phi = 2.0 * Math.PI * Math.random();
  return [sinTheta * Math.cos(phi), sinTheta * Math.sin(phi), cosTheta];
}

random on sphere with p5.js

球面内の一様乱数

function randomInSphere() {
  const cosTheta = -2.0 * Math.random() + 1.0;
  const sinTheta = Math.sqrt(1.0 - cosTheta * cosTheta);
  const phi = 2.0 * Math.PI * Math.random();
  const radius = Math.pow(Math.random(), 1.0 / 3.0);
  return [radius * sinTheta * Math.cos(phi), radius * sinTheta * Math.sin(phi), radius * cosTheta];
}

random in sphere with p5.js

導出

任意の確率密度関数に従う乱数を生成するためには逆関数法を用います。
逆関数法では対象の確率密度関数から累積密度関数の逆関数を求めます。$[0, 1]$の一様乱数を求めた逆関数で変換すると対象の確率密度関数に従った乱数になります。

[逆関数法を用いた乱数生成の証明と例 | 高校数学の美しい物語]
(https://mathtrain.jp/invsampling)

(筆者は数学にそこまで詳しくないので、以下の内容は本当に正しいのかあまり自信がないです...)

円周上の一様乱数

変数$\theta$について考えます。変数$\theta$は$[0, 2\pi]$の一様分布になることが想定できるので、確率密度関数$p(\theta)$は次のようになると考えられます($a$は定数)。

p(\theta) = a

確率密度関数の対象区間の積分は$1$になることから、$a$を求めると次のようになります。

\int_0^{2\pi} p(\theta) d\theta =
\int_0^{2\pi} a d\theta =
\left[a\theta \right]^{2\pi}_0 =
2\pi a = 1 \\
a = \frac{1}{2\pi} 

これにより変数$\theta$の確率密度関数$p(\theta)$は次のようになることがわかりました。

p(\theta) = \frac{1}{2\pi}

確率密度関数$p(\theta)$から累積密度関数$P(\theta)$を求めます。

P(\theta) =
\int_0^{\theta} p(\theta') d\theta' =
\int_0^{\theta} \frac{1}{2\pi} d\theta' =
\left[\frac{1}{2\pi}  \theta' \right]^{\theta}_0 =
\frac{1}{2\pi} \theta

$u=P(\theta)$とすると、

u = \frac{1}{2\pi} \theta \\
\theta = 2\pi u

となるので、累積密度関数$P(\theta)$の逆関数$P^{-1}(u)$は次のようになります。

P^{-1}(u) = 2\pi u

これにより円周上に一様に分布する点を求めるには$[0, 1]$の一様乱数$u$を用いて$\theta = 2\pi u$とすればよいことがわかりました。

円周内の一様乱数

円周内の一様乱数では２つの変数$\theta, r$があります。$\theta$は円周上の一様乱数で求めたものと同じなので、ここでは$r$について考えます。

半径$r$の円の円周は$2\pi r$となるように、円周は円の半径に比例します。このことから$r$に関する確率密度関数$p(r)$は次のように$r$に比例することが想定されます($a$は定数)。

p(r) = ar

$a$を求めると、

\int_0^{1} p(r) dr = \int_0^{1} ar dr = \left[\frac{1}{2}ar^2 \right]^{1}_0 = \frac{1}{2}a = 1 \\
a = 2

となるので確率密度関数$p(r)$は次のようになります。

p(r) = 2r

確率密度関数$p(r)$の累積密度関数$P(r)$は次のようになります。

P(r) =
\int_0^{r} p(r') dr' =
\int_0^{r} 2r' dr' =
\left[r'^2 \right]^{r}_0 =
r^2

$u = P(r)$とすると、

u = r^2 \\
r = \sqrt{u}

より、累積密度関数$P(r)$の逆関数$P^{-1}(u)$となります。

P^{-1}(u) = \sqrt{u}

これにより円周内に一様に分布する点を求めるのに、$[0, 1]$の一様な乱数$u$に用いて$r =\sqrt{u}$とすればよいことがわかりました。

球面上の一様乱数

変数$\phi$は円周上の一様乱数と同じ方法で導出できるので省略して、ここでは変数$\theta$について考えます。

任意の$\theta$に対して、XY平面に射影される円の半径は$\sin\theta$となります。複数の円に対して円周上に同じ密度で点を打つには半径に比例させればよいので、$\theta$の確率密度関数$p(\theta)$は次のようになると考えられます。

p(\theta) = a\sin\theta

$a$を求めると、

\int_0^{\pi} p(\theta) d\theta =
\int_0^{\pi} a\sin\theta d\theta =
\left[-a\cos\theta \right]^{\pi}_0 =
2a = 1 \\
a =\frac{1}{2} 

より確率密度関数$p(\theta)$は次のようになります。

p(\theta) = \frac{1}{2} \sin\theta

確率密度関数$p(\theta)$に対する累積密度関数$P(\theta)$は次のようになります。

P(\theta) =
\int_0^{\theta} p(\theta') d\theta' =
\int_0^{\theta} \frac{1}{2} \sin\theta' d\theta' =
\left[-\frac{1}{2} \cos\theta' \right]^{\theta}_0 =
-\frac{1}{2}\cos\theta + \frac{1}{2}

$u = P(\theta)$とすると、

\cos \theta = - 2u + 1

となり、$[0, 1]$の一様乱数$u$を用いると$\cos \theta = - 2u + 1$となることがわかりました。

球面内の一様乱数

ここでは変数$r$について考えます。

半径$r$の球の表面積は$4\pi r^2$となるように$r^2$に比例します。このことから$r$に関する確率密度関数$p(r)$は次のようになることが想定されます。

p(r) = ar^2

$a$を求めると、

\int_0^{1} p(r) dr =
\int_0^{1} ar^2 dr =
\left[\frac{1}{3}ar^3 \right]^{1}_0 =
\frac{1}{3}a =
1 \\
a = 3

となることから確率密度関数$p(r)$は次のようになることがわかりました。

p(r) = 3r^2

確率密度関数$p(r)$に対する累積密度関数$P(r)$次のようになります。

P(r) =
\int_0^{r} p(r') dr' =
\int_0^{r}  3r^2 dr' =
\left[r'^3 \right]^{r}_0 =
r^3

$u = P(r)$とすると

u = r^3 \\
r = \sqrt[3]{u}

となります。これより、$[0, 1]$の一様乱数$u$を用いると$r = \sqrt[3]{u}$となることがわかりました。