Main Content

ヒストグラムの平滑化

この例では、Curve Fitting Toolbox™ のスプライン コマンド群を使用して、ヒストグラムを平滑化する方法を示します。

次に、ある測定で収集されたデータを表す乱数値のヒストグラムを示します。

y = randn(1,5001);
hist(y);

Figure contains an axes object. The axes object contains an object of type patch. This object represents y.

このヒストグラムから、基の分布のより滑らかな近似を引き出します。このため、各バー区間の平均値がそのバーの高さと等しくなるスプライン関数 f を作成します。

これらのバーの 1 つの高さを h、その左端と右端を LR としたとき、次を満たすスプライン f を求めます。

integral {f(x) : L < x < R}/(R - L) = h,

あるいは、Ff の不定積分、つまり DF = f とすると、

F(R) - F(L) = h*(R - L).

[heights,centers] = hist(y);
hold on
ax = gca;
ax.XTickLabel = [];
n = length(centers);
w = centers(2)-centers(1);
t = linspace(centers(1)-w/2,centers(end)+w/2,n+1);
p = fix(n/2);
fill(t([p p p+1 p+1]),[0 heights([p p]),0],'w')
plot(centers([p p]),[0 heights(p)],'r:')
h = text(centers(p)-.2,heights(p)/2,'   h');
dep = -70;
tL = text(t(p),dep,'L');
tR = text(t(p+1),dep,'R');
hold off

Figure contains an axes object. The axes object contains 6 objects of type patch, line, text. This object represents y.

したがって、n をバーの数、t(i)i 番目のバーの左端、dt(i) をその幅、h(i) をその高さとすると、求めるのは

F(t(i+1)) - F(t(i)) = h(i) * dt(i), for i = 1:n,

あるいは、任意に F(t(1)) = 0 と設定すると次のようになります。

F(t(i)) = sum {h(j)*dt(j) : j=1:i-1}, for i = 1:n+1.

dt = diff(t);
Fvals = cumsum([0,heights.*dt]);

2 つの端点条件 DF(t(1)) = 0 = DF(t(n+1)) にこれを追加すると、完全 3 次スプライン内挿として F を求める必要のあるすべてのデータが揃います。

F = spline(t, [0, Fvals, 0]);

2 番目の引数にもう 2 つゼロ値がありますが、これはゼロ端点の傾き条件を示します。

最後に、スプライン F の微分 f = DF が、ヒストグラムの平滑化されたバージョンとなります。

DF = fnder(F);  % computes its first derivative
h.String = 'h(i)';
tL.String = 't(i)';
tR.String = 't(i+1)';
hold on
fnplt(DF, 'r', 2)
hold off
ylims = ylim; 
ylim([0,ylims(2)]);

Figure contains an axes object. The axes object contains 7 objects of type patch, line, text. This object represents y.