Suavizar un histograma

Este ejemplo muestra cómo usar los comandos de splines de Curve Fitting Toolbox™ para suavizar un histograma.

A continuación se muestra un histograma de valores aleatorios que podría representar datos recopilados en una medición.

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

Nos gustaría derivar de este histograma una aproximación más suave a la distribución subyacente. Lo hacemos construyendo una función de spline f cuyo valor promedio en cada intervalo de barras es igual a la altura de esa barra.

Si h es la altura de una de estas barras, y sus bordes izquierdo y derecho están en L y R, buscamos que el spline f satisfaga

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

o, con F como la integral indefinida de f, es decir, 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

De este modo, con n como el número de barras, t(i) como el borde izquierdo de la i-ésima barra, dt(i) como su anchura y h(i) como su altura, buscamos

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

o, estableciendo de forma arbitraria 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]);

Añada a esto las dos condiciones de contorno DF(t(1)) = 0 = DF(t(n+1)), y tenemos todos los datos que necesitamos para obtener F como interpolación por splines cúbicos completa.

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

Los dos valores de cero adicionales en el segundo argumento indican las condiciones de cero de la pendiente final.

Por último, la derivada f = DF del spline F es la versión suavizada del histograma.

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)]);