円内に一様分布する乱数を生成する時は、俺、絶対忘れないよヤコビアンのこと

頭出し

ある半径Rの円内に一様分布する乱数を生成する時には注意しないといけないことがありますよというお話。所謂「一度はやってしまうミス」系でもある。この手の話は円に限ったわけではなく、円の高次元版である球、あるいは超球(次元>3)、あるいは任意の座標変換をかませてそこにヤコビアンが出て来るときでも同じ。

使うライブラリは以下の２つ。なければinstall.packages関数でインストールしておく。

library(ggplot2)
library(dplyr)

また、以下のように定数を２つ定義しておく。意味はコメントにある通りだ。

#サンプル数
N <- 10^4
#半径のサイズ
R <- 4

本題

さて、問題のある半径Rの円内に一様分布する乱数を生成するにはどうしたらいいのかというと、非常に単純に考えた場合、以下のように(俺は)思考する。

X方向の成分として[-R, R]の間の一様乱数を生成する
Y方向の成分として[-R, R]の間の一様乱数を生成する
X^2 + Y^2 < R^2を満たすようなデータだけを残す

こうすることで、半径R内に一様に分布した乱数が得られそうだぞと、で、実際にやってみる。

#一様乱数から生成
df <- data.frame(x=runif(N, min=-R, max=R), y=runif(N, min=-R, max=R)) %>% filter(sqrt(x^2+y^2)<R)
ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")

散布図を見る限り、確かに一様に分布していそうだぞと、これでいいぞと。

ただ、このやり方はちょっとだけ無駄があって、「X^2 + Y^2 < R^2」を満たすデータのみを取得するというフィルタリングの処理を噛ませているために、発生させた乱数の全てを使っているわけではなく、一部（円からはみ出てしまった部分）を捨ててしまっている。そこでもうちょっと効率的なやり方を考えたい。高校数学でお勉強した"極座標"を思い出すとx, y座標と動径方向r、角度Θは

$x = r\cos(\theta)$
$y = r\sin(\theta)$

という関係で結びつけられていた。従って、円内に一様に分布する乱数を得るために

動径方向rの成分として[0, R]の間の一様乱数を生成する
角度Θ方向の成分として[0, 2π]の間の一様乱数を生成する
XとYを上記の数式に従って変換し、生成する

というようにしてやれば半径R内に一様に分布した乱数が得られそうだぞと、で、実際にやってみる。

#動径・角度方向に分けてサンプリング
theta <- runif(N, min=0, max=2*pi)
r <- runif(N, min=0, max=R)
df <- data.frame(x=r*cos(theta), y=r*sin(theta))
ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")

あ、あれ？この結果の図を見ると明らかに中心部分の方にサンプルが集中していて、外側が手薄になっている、つまり、一様に分布していないぞと、これはよくないぞということが散布図からわかる。
この原因は今度は大学の数学で習う、重積分の変数変換公式における面積(要)素の変換、

$dxdy = r dr d\theta$

のrを完全に見ていない、つまりそのヤコビアンを無視しているからだ*1。従って、動径方向に対しては半径rに比例するように乱数を引いてやらないといけなくて、[0, R]の一様な乱数uを逆関数数を用いて、

$r = \sqrt{2u}$

と変換して動径方向の乱数rを生成すればよいことがわかる。これは動径方向の累積分布関数P(R)が

$P(R) = \int_0^{R} r dr = \frac{1}{2}R^2$

とかける事から、このP(R)の逆関数をとることで証明される。

これを実際にやってみると

theta <- runif(N, min=0, max=2*pi)
r <- sqrt(2*runif(N, min=0, max=R))
df <- data.frame(x=r*cos(theta), y=r*sin(theta))
ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")

散布図より、一様に分布してるのがわかる。

追記(2015/1/18)

ほんぎゃああああああああああああああ！！！えらいミスしてて、上の図を良く見ると半径３くらいの円にしかなってなくて、俺が欲しかった、半径４の円になっていない。これを修正するためには、上のコードの

r <- sqrt(2*runif(N, min=0, max=R))

となっている箇所を、

r <- sqrt(2*runif(N, min=0, max=0.5*R^2))

とすればOK。これは、uの最大値とRの最大値を調正するためにmax引数をいじる必要があったということだ。これでPLOTすれば心穏やかに正しい結果となる。

*1:そもそも"一様"というのは考えている空間における抽象的な意味での面積素上に一様に分布するということなんで、こいつを無視してはいかん