はじめに

あるデータセットの分布(確率密度関数)を見たいときに、ヒストグラムを作るのは面倒くさい。最初に最小値と最大値を調べないといけないし、分散の小さい正規分布などでは、中央付近のビンは細かく、その他は粗くとらないときれいな分布が見えない。そういう場合は確率密度関数ではなく、累積分布関数を見るほうが楽だしきれいだ。累積分布関数はソート一発で求めることができる。

方法

データが配列dataに入っているとする。これをソートして、

(データの値, そのデータの位置の割合)

をプロットすると、それが累積分布関数になる。コードにするとこんな感じ。

# 配列dataにデータが入っている
data.sort!
n = data.size
n.times do |i|
  puts "#{data[i]} #{i.to_f/n.to_f}" # これが累積分布関数になる
end

原理

確率密度関数(probability density function、PDF) $f(x)$の原始関数 $F(x)$ が累積分布関数(cumulative distribution function, CDF)であり、以下のように定義される。

F(x) = \int_{-\infty}^x f(x) dx

累積分布関数の定義から明らかなように、ある確率変数$X$が、ある値$x$以下の値を取る確率が$F(x)$で与えられる。いま、$N$個のサンプルがあったとして、ある値$x$より小さいデータが$m$個あったとするならば、この確率変数の累積分布関数$F(x)$は$m/N$と推定できる。データが昇順にソートしてあれば、あるデータが$i$番目にあったならば、そのデータを含めてそれより小さいデータは$i$個あるのだから$i/N$が累積分布関数の推定値を与える。

サンプル

いくつかの分布について例を挙げる。

一様分布

rand()をそのままdataにpushすれば、一様分布になる。

  data = []
  10000.times do
    data.push rand()
  end

この時、分布関数は$f(x) = 1$であるので、累積分布関数は$F(x) = x$となる。プロットするとこんな感じ。

赤点がデータ、黒の実線が理論曲線。

指数分布

一様乱数を$r$として、$-\log(1-r)/\lambda$とすると、パラメータ$\lambda$の指数分布となる。

  data.push -Math::log(1.0-rand())/lambda

分布関数は$f(x) = \lambda e^{-\lambda x}$、累積分布関数は$F(x) = 1 - e^{-\lambda}$となる。

ガウス分布

一様乱数からガウス分布を作るには、ボックス＝ミュラー法を使う。一様乱数$r_1, r_2$から、平均0、標準偏差$\sigma$のガウス分布を作るには$\sigma \sqrt{-2 \log(r_1)} \cos(2 \pi r_2)$とすれば良い。

  r1 = rand()
  r2 = rand()
  data.push sigma*Math::sqrt(-2.0*Math::log(r1))*Math::cos(2.0*Math::PI*r2)

ソース

上記３つのデータを作る手抜きRubyスクリプトも掲載しておく。

dist.rb

def uniform(param)
  rand()
end

def exp_dist(lambda)
  -Math::log(1.0-rand())/lambda
end

def gaussian(sigma)
  r1 = rand()
  r2 = rand()
  sigma*Math::sqrt(-2.0*Math::log(r1))*Math::cos(2.0*Math::PI*r2)
end

def make_dist(dist,param,filename)
  puts filename
  data = []
  10000.times do 
    data.push dist.call(param)
  end
  data.sort!
  n = data.size
  open(filename,"w") do |f|
    n.times do |i|
      f.puts "#{data[i]} #{i.to_f/n.to_f}"
    end
  end
end

make_dist(method(:uniform),nil,"uniform.dat")
make_dist(method(:exp_dist),2.0,"exp.dat")
make_dist(method(:gaussian),2.0,"gaussian.dat")

上記をプロットするgnuplotファイル

plot.plt

set term png 
set out "uniform.png"
set xla "x"
set yla "F(x)"
set ytics 0.2
set xtics 0.2
unset key
p "uniform.dat" pt 7 lc rgb "red"\
, x lc rgb "black" lw 2

set out "exp.png"
set xtics 1.0
p "exp.dat" pt 7 lc rgb "red"\
, 1-exp(-x*2) lc rgb "black" lw 2

set out "gaussian.png"
set xtics 2.0
p "gaussian.dat" pt 7 lc rgb "red"\
, (1+erf(x/sqrt(2)/2))/2 lc rgb "black" lw 2

以下のようにすればグラフが描画される。

$ ruby dist.rb
$ gnuplot plot.plt

まとめ

データがどういう分布をしているか見るのには、とりあえず累積分布関数を見るのが良いと思う。ソートするだけだし、確率密度関数を見るよりきれいな場合が多い。確率密度関数が欲しいときには、累積分布関数を数値微分して出したほうが、ヒストグラムから作るよりもきれいで簡単だと思う。

どうでもいいけど

誤差関数の定義って

\mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt

なんだけどさ、なんで

\mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2/2} dt

としなかったんだろうね。おかげで上記のプロットファイルでも、ガウス分布の理論曲線描く時に、

p "gaussian.dat" pt 7 lc rgb "red"\
, (1+erf(x/sqrt(2)/2))/2 lc rgb "black" lw 2

って余計なsqrt(2)が必要になって面倒くさい。

累積分布関数をソートで求める