倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

単位球体内の一様分布 : 3次元

単位球体内の一様分布を生成する方法を見ていくシリーズ(目次)。 今回は3次元。 やることは2次元の場合と同じ。 ちなみに3次元の単位球体は普通の「球」です。

極座標

3次元の極座標は以下のように与えられるのでした:

  { \displaystyle\begin{align*}
    &\begin{cases}
        x = r \sin\theta \cos\varphi \\
        y = r \sin\theta \sin\varphi \\
        z = r\cos\theta
    \end{cases} &
    \begin{pmatrix}
        0 \le r < \infty \\
        0 \le \theta \le \pi \\
        0 \le \varphi < 2\pi
    \end{pmatrix}
\end{align*}}

単位球内での積分は以下のようになります:

  { \displaystyle\begin{align*}
    \int_0^1 r^2dr \int_0^\pi \sin\theta d\theta \int_0^{2\pi} d\varphi
\end{align*}}

一様分布する変数

一様分布する変数への変数変換は以下の通り:

  { \displaystyle\begin{align*}
  (r,\, \theta,\, \varphi )
    \rightarrow (R,\, \Theta,\, \varphi)
    = (\tfrac{1}{3} r^3,\, \cos\theta,\, \varphi ) \in [0,\, \tfrac{1}{3}] \times [-1,\, 1] \times [0,\,2\pi )
\end{align*}}

このとき、単位球体内の積分は以下のようになります:

  { \displaystyle\begin{align*}
    \int_0^{1/3} dR \int_{-1}^1 d\Theta \int_0^{2\pi} d\varphi
\end{align*}}

直交座標の計算

極座標の変数を一様分布する変数で表すと

  { \displaystyle\begin{align*}
    (r,\, \theta,\, \varphi) = (\sqrt[3]{3R},\, \cos^{-1}\Theta,\, \varphi)
\end{align*}}

となるので、直交座標は以下のようになります:

  { \displaystyle\begin{align*}
    \begin{cases}
        x = \sqrt[3]{3R} \sqrt{1-\Theta^2} \cos\varphi \\
        y = \sqrt[3]{3R} \sqrt{1-\Theta^2} \sin\varphi \\
        z = \sqrt[3]{3R} \Theta
    \end{cases}
\end{align*}}

これを用いて3次元球体内の一様乱数を生成するには、効率化のために { R } に3を吸収させて(定義域にも)、以下のようにします:

  { \displaystyle\begin{align*}
    \begin{cases} x = \sqrt[3]{\tilde{R}} \sqrt{1-\tilde{\Theta}^2} \cos\tilde{\varphi} \\ y = \sqrt[3]{\tilde{R}} \sqrt{1-\tilde{\Theta}^2} \sin\tilde{\varphi} \\ z = \sqrt[3]{\tilde{R}} \tilde{\Theta} \end{cases} \begin{pmatrix} 0 \le R \le 1 \\ -1 \le \Theta \le 1 \\ 0 \le \varphi < 2\pi \end{pmatrix}
\end{align*}}

ただし、変数 { X \in [a,\,b] } に対して { \tilde{X} }区間 { [a,\,b] } 間に一様分布する変数とします。

Java コード

実際に Java で3次元球体内の一様乱数を生成してみましょう。 Java コードは以下のようになります:

import static java.lang.Math.*;

public class UniformDistributionB3{

    public static void main(String... args){

        for(int i = 0; i < 30000; i++){
            double r = cbrt(random());
            double cosTheta = random() * 2 - 1;
            double sinTheta = sqrt(1- cosTheta*cosTheta);
            double phi = random() * 2 * PI;
            
            double x = r * sinTheta * cos(phi);
            double y = r * sinTheta * sin(phi);
            double z = r * cosTheta;
            
            System.out.println(x+" "+y+" "+z);
        }
    }
}

結果をプロットすると下図のようになります: