julia のヒストグラムに、カーネル密度推定からの予測曲線を重ねる方法

julia でヒストグラムを書いて、カーネル密度推定からの予測度数の曲線を重ねてみます。

カーネル密度推定は、KernelDensity.jl を用いて行います。推定値は、確率 (0から1の値) で出てきます。ヒストグラムに重ねる時には、ヒストグラムの最も観測数が多かった階級の観測数を求める必要があります。この直前のブログの fit(Histogram, ...) というのを使います。

PyPlot と seaborn を呼び出して、 sns.distplot(data) で終りでしょう。と、おっしゃる方がいらっしゃると思います

そのとおりです。上記の図をヒストグラムと曲線に分けて書いてみようというわけです。

では、パッケージとランダムサンプルの作製をします。サンプルは、2次元ですが、1列目だけを今回は使います。グラフは最近注目の Plots.jl を使うことにします。 

using Plots
using KernelDensity
using Distributions
using StatsBase

data = vcat(
  rand(Normal(),10000,2),
  rand(Normal(3,0.8),10000,2),
  rand(Normal(-2,2),10000,2).+[0 2],
  rand(TDist(10),1000,2)
  )[rand(1:31000,100000) ,:];

xdata = data[:,1];

まずは、fit(Histogram, ...) を使って最頻階級の度数 xn_max と、観測値の最小値 x_min と、最大値 x_max を求めておきます。カーネル密度推定は、観測範囲外ではできないので、x_mimは大きめ, x_min は小さめにしておきます 

x_hist = fit(Histogram, xdata, nbins=100)
xn_max = maximum(x_hist.weights)
x_min = floor(minimum(xdata))+1
x_max = floor(maximum(xdata))

カーネル密度推定を、KernelDensity.jl を用いて行います。実際に曲線を描く時は、観測値に対して推定値を返すのではなく、グラフのx軸上の値に対して、推定値を求めていきます。カーネル密度推定を行う x_range の範囲を決めて、0.1 きざみで、推定値を求め、上で求めた xn_max をかけてヒストグラムと高さを合せます。 

x_range = [x_min:0.1:x_max;];
x_kde = kde(xdata)
x_pdf = pdf(x_kde, x_range);
x_ndf = x_pdf./(maximum(x_pdf)) .* xn_max

予測度数の曲線が (x_range, x_ndf) で書けるようになりました。 

histogram(xdata, nbins=100, xlims = (x_min-1, x_max+1))
plot!(x_range, x_ndf, lw=3)

以上です。少しだけ、カーネル密度推定が理解できたような気がします。

B! LINE