El comando CSAPS

El comando csaps proporciona el spline de suavizado. Se trata de un spline cúbico que sigue más o menos la tendencia subyacente supuesta de datos ruidosos. El parámetro de suavizado que elija determinará cuánto se aproxima el spline de suavizado a los datos proporcionados. Aquí le presentamos la información básica, una versión abreviada de la documentación:

CSAPS Cubic smoothing spline.

VALUES = CSAPS(X, Y, P, XX)

Returns the values at XX of the cubic smoothing spline for the

given data (X,Y) and depending on the smoothing parameter P, chosen

from the interval [0 .. 1]. This smoothing spline f minimizes

P * sum_i W(i)(Y(i) - f(X(i)))^2 + (1-P) * integral (D^2 f)^2

Ejemplo: Datos ruidosos de un polinomio cúbico

Estas son algunas pruebas llevadas a cabo. Se comienza con los datos de un cúbico simple, q(x) := x^3, se contaminan los valores con ruido y se determina el valor del parámetro de suavizado en 0,5. Entonces se presentan los valores de suavizado resultantes junto con el cúbico subyacente y los datos contaminados.

xi = (0:.05:1);
q = @(x) x.^3;
yi = q(xi);
randomStream = RandStream.create( 'mcg16807', 'Seed', 23 );
ybad = yi+.3*(rand(randomStream, size(xi))-.5);
p = .5;
xxi = (0:100)/100;
ys = csaps(xi,ybad,p,xxi);
plot(xi,yi,':',xi,ybad,'x',xxi,ys,'r-')
title('Clean Data, Noisy Data, Smoothed Values')
legend( 'Exact', 'Noisy', 'Smoothed', 'Location', 'NorthWest' )

El suavizado aquí está demasiado exagerado. Eligiendo un valor del parámetro de suavizado p más cercano a 1, obtenemos un spline de suavizado más aproximado a los datos proporcionados. Probamos con p = .6, .7, .8, .9, 1 y representamos los splines de suavizado resultantes.

yy = zeros(5,length(xxi));
p = [.6 .7 .8 .9 1];
for j=1:5
   yy(j,:) = csaps(xi,ybad,p(j),xxi);
end
hold on
plot(xxi,yy);
hold off
title('Smoothing Splines for Various Values of the Smoothing Parameter')
legend({'Exact','Noisy','p = 0.5','p = 0.6','p = 0.7','p = 0.8', ...
        'p = 0.9', 'p = 1.0'}, 'Location', 'NorthWest' )

Vemos que el spline de suavizado varía mucho según el parámetro de suavizado elegido. Incluso con p = 0,9, el spline de suavizado sigue estando muy alejado de la tendencia subyacente, mientras que con p = 1, la interpolación se acerca a los datos (ruidosos).

De hecho, la formulación utilizada por csapi (pág. 235ff de A Practical Guide to Splines) es muy sensible al escalado de la variable independiente. Un simple análisis de las ecuaciones utilizadas revela que el rango de sensibilidad de p está en torno a 1/(1+epsilon), siendo epsilon := h^3/16 y h la diferencia media entre los sitios vecinos. Sobre todo, podemos esperar una mayor aproximación a los datos cuando p = 1/(1+epsilon/100) y un suavizado más satisfactorio cuando p = 1/(1+epsilon*100).

La representación siguiente muestra el spline de suavizado para valores de p cuando se acerca a este número mágico, 1/(1+epsilon). En este caso, sería más conveniente observar 1-p, ya que el número mágico, 1/(1+epsilon), se acerca mucho a 1.

epsilon = ((xi(end)-xi(1))/(numel(xi)-1))^3/16;
1 - 1/(1+epsilon)

ans = 7.8124e-06

plot(xi,yi,':',xi,ybad,'x')
hold on
labels = cell(1,5);
for j=1:5
   p = 1/(1+epsilon*10^(j-3));
   yy(j,:) = csaps(xi,ybad,p,xxi);
   labels{j} = ['1-p= ',num2str(1-p)];
end
plot(xxi,yy)
title('Smoothing Splines for Smoothing Parameter Near Its ''Magic'' Value')
legend( [{'Exact', 'Noisy'}, labels], 'Location', 'NorthWest' )
hold off

En este ejemplo, el spline de suavizado es muy sensible a la variación del parámetro de suavizado cerca del número mágico. Parece que la mejor opción es la que se aleja más del 1, pero puede que quiera elegir una más alejada.

p = 1/(1+epsilon*10^3);
yy = csaps(xi,ybad,p,xxi);
hold on
plot( xxi, yy, 'y', 'LineWidth', 2 )
title( sprintf( 'The Smoothing Spline For 1-p = %s is Added, in Yellow', num2str(1-p) ) )
hold off

También puede proporcionar ponderaciones de errores a csaps para prestar más atención a ciertos puntos de datos. Además, si no proporciona los sitios de evaluación xx, entonces csaps devolverá el spline de suavizado en formato ppform.

Por último, csaps también puede gestionar datos con valor vectorial e, incluso, datos de malla y multivariados.

El comando SPAPS

El spline cúbico de suavizado proporcionado por el comando spaps difiere de aquel construido en csaps únicamente en la forma en que se selecciona. Aquí le presentamos una versión abreviada de la documentación de spaps:

SPAPS Smoothing spline.

[SP,VALUES] = SPAPS(X,Y,TOL) returns the B-form and, if asked,

the values at X, of a cubic smoothing spline f for the given

data (X(i),Y(:,i)), i=1,2, ..., n.

The smoothing spline f minimizes the roughness measure

F(D^2 f) := integral ( D^2 f(t) )^2 dt on X(1) < t < X(n)

over all functions f for which the error measure

E(f) := sum_j { W(j)*( Y(:,j) - f(X(j)) )^2 : j=1,...,n }

is no bigger than the given TOL. Here, D^M f denotes the M-th

derivative of f. The weights W are chosen so that E(f) is the

Composite Trapezoid Rule approximation for F(y-f).

f is constructed as the unique minimizer of

rho*E(f) + F(D^2 f),

with the smoothing parameter RHO so chosen so that E(f) equals

TOL. Hence, FN2FM(SP,'pp') should be (up to roundoff) the same

as the output from CPAPS(X,Y,RHO/(1+RHO)).

Tolerancia frente a parámetro de suavizado

Puede resultar más sencillo proporcionar una tolerancia adecuada a spaps en lugar del parámetro de suavizado p requerido por csaps. En el ejemplo anterior, hemos añadido ruido aleatorio distribuido uniformemente del intervalo 0.3*[-0.5 .. 0.5]. Por tanto, podemos estimar un valor razonable para tol como valor de la medida de error e en dicho ruido.

tol = sum((.3*(rand(randomStream, size(yi))-.5)).^2);

Esta representación muestra el spline de suavizado resultante que se ha construido mediante spaps. Tenga en cuenta que se ha especificado una ponderación de errores uniforme, que es el valor predeterminado en csaps.

[sp,ys,rho] = spaps(xi,ybad,tol,ones(size(xi)));
plot(xi,yi,':',xi,ybad,'x',xi,ys,'r-')
title( sprintf( 'Clean Data, Noisy Data, Smoothed Values (1-p = %s )', num2str(1/(1+rho)) ) );
legend( {'Exact','Noisy','Smoothed'}, 'location', 'NorthWest' )

El título de la figura muestra el valor de p que se usaría en csaps para obtener exactamente este spline de suavizado para estos datos.

Además, también se muestra el spline de suavizado obtenido de csaps cuando no se proporciona un parámetro de suavizado. En este caso, csaps elige el parámetro de acuerdo con un procedimiento ad hoc determinado cuyo objetivo es localizar la región en la que el spline de suavizado es más sensible al parámetro de suavizado (similar a como hemos comentado antes).

hold on
plot(xxi,fnval(csaps(xi,ybad),xxi),'-')
title('Clean Data, Noisy Data, Smoothed Values')
legend({'Exact' 'Noisy' 'spaps, specified tolerance' ...
        'csaps, default smoothing parameter'}, 'Location', 'NorthWest' )
hold off

CSAPS frente a SPAPS

Los comandos csaps y spaps difieren en la forma en la que se especifica un spline de suavizado determinado, mediante un parámetro de suavizado frente a una tolerancia. Otra diferencia es que spaps puede proporcionar un spline de suavizado lineal o de quinto grado, además del spline cúbico de suavizado.

El spline de suavizado de quinto grado es mejor que el spline cúbico de suavizado en situaciones en las que se desee que la segunda derivada se mueva lo menos posible.