ヒストグラムの平滑化
この例では、Curve Fitting Toolbox™ のスプライン コマンド群を使用して、ヒストグラムを平滑化する方法を示します。
次に、ある測定で収集されたデータを表す乱数値のヒストグラムを示します。
y = randn(1,5001); hist(y);
このヒストグラムから、基の分布のより滑らかな近似を引き出します。このため、各バー区間の平均値がそのバーの高さと等しくなるスプライン関数 f
を作成します。
これらのバーの 1 つの高さを h
、その左端と右端を L
と R
としたとき、次を満たすスプライン f
を求めます。
integral {f(x) : L < x < R}/(R - L) = h
,
あるいは、F
を f
の不定積分、つまり 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
したがって、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)]);