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