球面上に一様分布する乱数の生成

世の中には様々な事情が存在し、超球面上に一様分布する乱数が必要になることがあります。
３次元の場合はすでに記事がありますが（[こちら](http://qiita.com/clomie/items/e5dd35dcfcba082b2a7f)）、本記事では一般の場合を考えます。
サンプルプログラムはpythonです。

ではやっていきましょう。


# 生成
いくつか手法を書き連ねます。

+ 単位超立方体上の一様分布からサンプルし、ノルムが1以下なら正規化して採用する

```python:python
import numpy as np

def f1(dim):
  while True:
    x = np.random.random(dim)*2.-1.
    r = np.linalg.norm(x)
    if r != 0. and r < 1.:
      return x/r
```

ノルムのチェックをしないと立方体の"角"にあたる部分の影響を受け、一様ではなくなってしまうので注意が必要です。
このアルゴリズムには次元の呪いという明確な欠点があります。単位超立方体中に単位超球が占める体積の割合は次元の増加に伴ってどんどん小さくなりますから、ノルムチェックで弾かれる確率はどんどん高くなり、結果ループが止まらなくなるというものです。実用的とはいえないでしょう。

<!--
+ 座標変換による方式

初めに紹介した[こちらの記事](http://qiita.com/clomie/items/e5dd35dcfcba082b2a7f)で採用されている方法は多次元に拡張できるでしょうか。
n次元単位球面上の点の極座標表示はたとえば

```math
\begin{align*}
\begin{cases}
 x_1 &= \cos\theta_1 \\
 x_2 &= \sin\theta_1 \cos\theta_2 \\
 &\vdots \\
 x_{n-1} &= \sin\theta_1\sin\theta_2 \dots \sin\theta_{n-2}\cos\phi \\
 x_n &= \sin\theta_1\sin\theta_2 \dots \sin\theta_{n-2}\sin\phi
\end{cases} & \quad
\begin{pmatrix}
0\le\theta_i\le\pi\\
0\le\phi < 2\pi\\
\end{pmatrix}
\end{align*}
```

のように書け、これのヤコビアンは$\prod_{i=1}^{n-2}\sin^{n-2-i}\theta_i$となります。


>>>>

``` python:python
import numpy as np

def f2(dim):
  x = np.ones(dim)
  for i in range(dim-2):
    z = np.random.random()*2.-1.
    x[i] *= z
    x[i+1:] *= np.sqrt(1.-z*z)
  phi = np.random.random()*2.*np.pi
  x[dim-2] *= np.cos(phi)
  x[dim-1] *= np.sin(phi)
  return x
```

初めに紹介した[こちらの記事](http://qiita.com/clomie/items/e5dd35dcfcba082b2a7f)で採用されている方法を多次元に拡張したものです。
n次元単位球面上の点の極座標表示を考えると、

```math
\begin{align*}
\begin{cases}
 x_1 &= \cos\theta_1 \\
 x_2 &= \sin\theta_1 \cos\theta_2 \\
 &\vdots \\
 x_{n-1} &= \sin\theta_1\sin\theta_2 \dots \sin\theta_{n-2}\cos\phi \\
 x_n &= \sin\theta_1\sin\theta_2 \dots \sin\theta_{n-2}\sin\phi
\end{cases} & \quad
\begin{pmatrix}
0\le\theta_i\le\pi\\
0\le\phi < 2\pi\\
\end{pmatrix}
\end{align*}
```

とかけ、これのヤコビアンは$\prod_{i=1}^{n-2}\sin\theta_i$
(→$\prod_{i=1}^{n-2}\sin^{n-i-1}\theta_i$なので以下に説明されるこの変数変換は一般の$n$次元ではうまくいきません。密度関数が$\sin^{n-i-1}\theta_i$に比例する分布から乱数を生成する必要が出てくるからです)
となります。
これは単に$\theta_i$を一様乱数からとると$\sin\theta_i$の分だけ歪んでしまうことを意味します。実際に一様乱数からサンプルしてみると、$\sin$関数において絶対値が小さい、すなわち$\theta_i$の変化に対し$x_i$の変化が小さくなる場所（３次元であれば極付近）が濃くなっていることが確認できます。
これを回避するため変数変換 $z_i=\cos\theta_i$を行います。

```math
\begin{align*}
\begin{cases}
 x_1 &= z_1 \\
 x_2 &= z_2\sqrt{1-z_1^2} \\
 &\vdots \\
 x_{n-1} &= \sqrt{1-z_1^2}\sqrt{1-z_2^2} \dots \sqrt{1-z_{n-2}^2}\cos\phi \\
 x_n &= \sqrt{1-z_1^2}\sqrt{1-z_2^2} \dots \sqrt{1-z_{n-2}^2}\sin\phi \\
\end{cases} &　\quad
\begin{pmatrix}
-1 \le z_i \le 1\\
0 \le \phi < 2\pi\\
\end{pmatrix}
\end{align*}
```

$dz_i = -\sin\theta_id\theta_i$であることからヤコビアンは$(-1)^{n-2}$となり、すなわち$z_i$と$\phi$を一様乱数からとれば、偏りのない分布、すなわち球面上の一様乱数が得られることがわかります。

この方法は球面上以外での分布にも応用することができますが、アルゴリズムとしては少々複雑になり、上の実装では$O(dim)$となります。
今回の目的では次の方法がよりよい選択でしょう。
-->

+ 多次元ガウス分布からのサンプルを正規化する

```python:python
import numpy as np

def f3(dim):
  while True:
    x = np.random.randn(dim)
    r = np.linalg.norm(x)
    if r != 0.:
      return x/r
```

零ベクトルを平均、単位行列を分散共分散行列に持つ多次元正規分布の密度関数は、$Z$を正則化項として$f(x)= \frac{1}{Z}\exp(-\frac{1}{2}x^\top x)$です。$x$のノルムを$r$とおけば$x^\top x=r^2$ですから、各点の生成確率はそのノルム$r$のみに依存すること、すなわち方向はランダムであることがわかります。よってサンプルを正則化し$r$の影響を除去することで球面上一様分布乱数が得られます。


# おまけ：超球内一様分布
超球の内部に一様分布する乱数も考えてみましょう。
表面上での分布は上記の方法でとれますから、これに加えて動径方向をどうするか考えれば十分です。
具体的には、ノルムが$r$の点が出現する確率が、動径が$r$の超球の表面積に比例するようにすればよいことがわかります。すなわち球面上一様乱数のノルムを、一様乱数$r$を用いて$\sqrt[n]{r}$に調整することで超球内一様乱数が得られます。

```python:python
import numpy as np
def g(dim):
  x = f3(dim)
  r = np.power(np.random.random(), 1./dim)
  return x*r
```

![figure_1-1.png](https://qiita-image-store.s3.amazonaws.com/0/148948/4c3a94fd-8407-483f-7cef-2f9805f3995c.png)