julia での度数分布の求め方

julia で、数値ベクトルの度数分布を求めてみましょう。文字列なら、StatsBase.countmap() で足りることが多いでしょう。

まずは、パッケージの読み込みと、テストサンプルデータを作ります。これからしばらく2次元データを使うので、2次元で作ります。 

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];
ydata = data[:,2];

はじめに、10000個ずつ3通りの正規分布を縦に並べて、最後の t 分布はノイズで入れています。そして、この中から、ランダムで選ぶ方法を採りました。ランダムで採り直す方が現実に近いと考えています。

1列目を xdata として、使います。

次に、xdata のヒストグラムを作ります。 

histogram(xdata, nbins=100)

このヒストグラムのそれぞれの階級のデータの数 (棒グラフの高さ) を求めたいということです。

それぞれのデータの個数は、fit(Histogram, xdata, ...) を使います。 

julia> x_hist = fit(Histogram, xdata, nbins=100)
┌ Warning: Default for keyword argument "closed" has changed from :right to :left.
│ To avoid this warning, specify closed=:right or closed=:left as appropriate.
│   caller = ip:0x0
└ @ Core :-1
Histogram{Int64,1,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}
edges:
  -10.0:0.2:6.2
weights: [6, 0, 4, 0, 2, 7, 8, 4, 7, 7  …  866, 626, 328, 159, 87, 75, 27, 13, 3, 3]
closed: left
isdensity: false

"x_hist.edges" などとして、それぞれの値を取り出します。

x_hist.edges は、階級の分け方です。-10 から、0.2 ずつ 6.2 までです。100階級を指定しましたが、81階級になっています。指定の仕方が、データに合わないのでしょう。julia は比較的、切れの良い階級を作ろうとする為と考えます。階級自身を指定する方法もありますが、これは、StatsBase.fit(Histogram, ... ) のマニュアルを見て下さいね。

この範囲指定を取り出す時は、さらに、 

julia> x_hist.edges[1]
-10.0:0.2:6.2

[1] を付けるようです。collect() ではダメなようです。

次に、それぞれの階級内のデータの数です。weights に格納されています。 

julia> x_hist.weights
81-element Array{Int64,1}:
   6
   0
   4
   0
   2
   7
   8
   4
   7
   7
   ⋮
 626
 328
 159
  87
  75
  27
  13
   3
   3

こちらは既に、Array ですので、そのまま使えるでしょう。

B! LINE