2013/11/23
PostgreSQL > ボックス=ミュラー法を SQL で（正規乱数）
PostgreSQL > ボックス=ミュラー法を SQL で（正規乱数）
昨日はモンテカルロ法の一例で円周率を近似したが、今日は正規分布にしたがう乱数を得るボックス＝ミュラー法を SQL で行う。同手法の詳細は下記。実行環境は 2013/11/02 を参照。
Wikipedia ボックス＝ミュラー法
上記のとおり、二つの乱数群を用意して sin か cos どちらか好きな式で計算すれば、標準正規分布に従う乱数群が得られる。平均と標準偏差を変えたければ、乱数に標準偏差を掛けて平均を足す。下が、cos を使って平均 2、標準偏差 4 の正規乱数 10個を得る様子。
使いやすいようストアド関数化した。引数で平均、標準偏差、乱数の数を渡す。
↓ 上記ストアド関数を使い、1000個の正規乱数を得て平均と標準偏差（不偏分散の平方根の方）を算出するのを10回繰り返した様子。繰り返しは再帰クエリ WITH RECURSIVE で行った。冒頭のように数を10個得るだけでは本当に正規分布に従っているか分からないが、下の結果を見ると平均と標準偏差は概ね合っている。
さらに R で可視化して確認するため、下記の PL/R のストアド関数を作った。引数で、平均、標準偏差、乱数の個数に加え、繰り返し数と出力ファイルパスを渡す。繰り返し数が NULL なら一回実行してヒストグラムを描き、本来の正規分布曲線を重ね描きする。繰り返し数が指定されたら、その数だけ実行し、各回の乱数群の平均と標準偏差をプロットする。
↓ 標準正規分布の乱数を 1000個得て分布を見たようす。ヒストグラムの棒が PostgreSQL による乱数。赤い曲線が、R の dnorm 関数と curve 関数で本来の標準正規分布に従う確率。確かに正規乱数になっている。
↓ 上と同様に、平均 2、標準偏差 4 で正規乱数 1000個を出して確認。
↓ 平均 2、標準偏差 4 の正規乱数 1000個得るのを 100回繰り返し、各回の数値群に対し平均と標準偏差を算出してプロットした。概ね本来の値を中心に分布している。
↓ 同じ正規乱数を、今度は 10000個得るのを 100回繰り返した。上の1000個の時より、本来の平均・標準偏差に近く寄っている。
以上、PostgreSQL でボックス=ミュラー法による正規乱数を得ることができた。昨日のモンテカルロ法もそうだが、ランダムな手法を PostgreSQL で少しずつ試していきたい。
Wikipedia ボックス＝ミュラー法
上記のとおり、二つの乱数群を用意して sin か cos どちらか好きな式で計算すれば、標準正規分布に従う乱数群が得られる。平均と標準偏差を変えたければ、乱数に標準偏差を掛けて平均を足す。下が、cos を使って平均 2、標準偏差 4 の正規乱数 10個を得る様子。
使いやすいようストアド関数化した。引数で平均、標準偏差、乱数の数を渡す。
CREATE OR REPLACE FUNCTION "201311"."23_box_muller"(float, float, int)
RETURNS setof float LANGUAGE sql IMMUTABLE AS $$
/* 任意の平均, 標準偏差の正規乱数を返す。ボックス＝ミュラー法。
引数 : 平均, 標準偏差, 返す乱数の数 */
SELECT sqrt(-2 * ln(random())) * cos(2 * pi() * random()) * $2 + $1
FROM generate_series(1, $3) t (foo) ;
$$ ;
-- Example
SELECT "201311"."23_box_muller"(0, 1, 1000) ;
SELECT "201311"."23_box_muller"(2, 4, 1000) ;
↓ 上記ストアド関数を使い、1000個の正規乱数を得て平均と標準偏差（不偏分散の平方根の方）を算出するのを10回繰り返した様子。繰り返しは再帰クエリ WITH RECURSIVE で行った。冒頭のように数を10個得るだけでは本当に正規分布に従っているか分からないが、下の結果を見ると平均と標準偏差は概ね合っている。
さらに R で可視化して確認するため、下記の PL/R のストアド関数を作った。引数で、平均、標準偏差、乱数の個数に加え、繰り返し数と出力ファイルパスを渡す。繰り返し数が NULL なら一回実行してヒストグラムを描き、本来の正規分布曲線を重ね描きする。繰り返し数が指定されたら、その数だけ実行し、各回の乱数群の平均と標準偏差をプロットする。
CREATE OR REPLACE FUNCTION "201311"."23_box_muller_plr"(
float, float, int, int, text)
RETURNS void LANGUAGE plr IMMUTABLE AS $$
# 正規乱数を得る関数 23_box_muller の可視化。
# 第一〜三引数はそのまま同関数に渡す（平均, 標準偏差, 個数）。
# 第四引数が繰り返し回数で、NULL ならヒストグラムを描画。
# NULL 以外なら指定回数繰り返して平均と標準偏差を描画。
graphics.off()
png(file=arg5, width=444, height=444)
if(is.null(arg4)){
sql = paste(sep='', 'SELECT "201311"."23_box_muller"('
, arg1, ',', arg2, ',', arg3, ')')
dat = pg.spi.exec(sql)[,1]
hist(dat, prob=T, plot=T, xlab="Value"
, main="Distribution by Box-Muller method (on PostgreSQL)")
curve(dnorm(x, mean=arg1, sd=arg2), col=2, add=T, xpd=NA)
} else {
sql = paste(sep='', '
WITH RECURSIVE a AS (
SELECT 0 cnt, ', arg4, ' lim
, NULL :: float avg, NULL :: float std
UNION ALL (
SELECT cnt + 1, lim, b.avg, b.std
FROM a, (
SELECT avg(num), stddev_samp(num) std
FROM "201311"."23_box_muller"('
, arg1, ',', arg2, ',', arg3, ') t (num)
) b
WHERE cnt < lim ORDER BY cnt DESC LIMIT 1
)
)
SELECT avg, std FROM a WHERE cnt > 0
')
dat = pg.spi.exec(sql)
ylm = range(c(dat$avg, dat$std))
par(mai=par("mai") * c(1, 1, 0, 1))
plot(dat$avg, pch=1, ylim=ylm, xlab="Count", ylab="avg or stddev_samp")
points(dat$std, pch=2)
}
dev.off()
$$
-- Example
SELECT "201311"."23_box_muller_plr"(0, 1, 1000, NULL, 'foo.png') ;
SELECT "201311"."23_box_muller_plr"(2, 4, 1000, 100, 'bar.png') ;
↓ 標準正規分布の乱数を 1000個得て分布を見たようす。ヒストグラムの棒が PostgreSQL による乱数。赤い曲線が、R の dnorm 関数と curve 関数で本来の標準正規分布に従う確率。確かに正規乱数になっている。
↓ 上と同様に、平均 2、標準偏差 4 で正規乱数 1000個を出して確認。
↓ 平均 2、標準偏差 4 の正規乱数 1000個得るのを 100回繰り返し、各回の数値群に対し平均と標準偏差を算出してプロットした。概ね本来の値を中心に分布している。
↓ 同じ正規乱数を、今度は 10000個得るのを 100回繰り返した。上の1000個の時より、本来の平均・標準偏差に近く寄っている。
以上、PostgreSQL でボックス=ミュラー法による正規乱数を得ることができた。昨日のモンテカルロ法もそうだが、ランダムな手法を PostgreSQL で少しずつ試していきたい。