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

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

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

生成

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

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

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

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

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

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

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