Команда CSAPS

Команда csaps обеспечивает сглаживающий сплайн. Это кубический сплайн, который более или менее следует предполагаемому основному тренду в шумных данных. Параметр сглаживания, выбираемый пользователем, определяет, насколько точно сглаживающий сплайн следует заданным данным. Ниже приведена основная информация, сокращенная версия документации:

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

Пример: Шумные данные из кубического многочлена

Вот несколько пробных прогонов. Мы начинаем с данных из простого кубического, q(x) := x^3, загрязните значения некоторым шумом и выберите значение параметра сглаживания, равное .5. Затем постройте график полученных сглаженных значений вместе с нижележащими кубическими и загрязненными данными.

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' )

Сглаживание здесь сильно переусердствует. Путем выбора параметра сглаживания p ближе к 1 получим сглаживающий сплайн ближе к заданным данным. Мы стараемся p = .6, .7, .8, .9, 1и постройте график полученных сглаживающих сплайнов.

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' )

Мы видим, что сглаживающий сплайн может быть очень чувствительным к выбору параметра сглаживания. Даже для p = 0,9, сглаживающий сплайн все еще далек от основного тренда, в то время как для p = 1, получаем интерполятор к (шумным) данным.

Фактически, состав, используемый csapi (стр. 235ff Практического руководства по сплайнам) очень чувствителен к масштабированию независимой переменной. Простой анализ используемых уравнений показывает, что чувствительный диапазон для p находится вокруг 1/(1+epsilon), с epsilon := h^3/16, и h средняя разница между соседними сайтами. В частности, вы ожидаете близкого следования данных, когда p = 1/(1+epsilon/100) и некоторое удовлетворительное сглаживание, когда p = 1/(1+epsilon*100).

На графике ниже показан сглаживающий сплайн для значений p рядом с этим волшебным числом 1/(1+epsilon). Для этого случая более информативно смотреть на 1-p начиная с магического числа, 1/(1+epsilon), очень близко к 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

В этом примере сглаживающий сплайн очень чувствителен к изменению параметра сглаживания вблизи магического числа. Тот, кто дальше от 1, кажется лучшим выбором из них, но вы можете предпочесть тот, который выше этого.

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

Вы также можете поставлять csaps с весами ошибок, чтобы уделять больше внимания некоторым точкам данных, чем другим. Кроме того, если вы не предоставляете сайты анализа xx, то csaps возвращает ppform сглаживающего сплайна.

Наконец, csaps может также обрабатывать векторные данные и даже многомерные данные с привязкой к сетке.

Команда SPAPS

Кубический сглаживающий сплайн, предоставляемый командой spaps отличается от построенного в csaps только таким образом, как он выбран. Вот сокращенная версия документации для 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)).

Допуск по сравнению с параметром сглаживания

Может быть проще предоставить подходящий допуск для spaps чем параметр сглаживания p требуется по csaps. В нашем предыдущем примере мы добавили равномерно распределенный случайный шум из интервала 0.3*[-0.5 .. 0.5]. Следовательно, мы можем оценить разумное значение для tol как значение измерения ошибки e при таком шуме.

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

На этом графике показан результирующий сглаживающий сплайн, созданный spaps. Обратите внимание, что веса ошибок указаны как однородные, что является их значением по умолчанию в 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' )

Заголовок рисунка показывает значение p вы бы использовали в csaps для получения именно этого сглаживающего сплайна для этих данных.

Здесь, кроме того, находится сглаживающий сплайн, предоставляемый csaps если параметр сглаживания не задан. В этом случае csaps выбирает параметр с помощью определенной специальной процедуры, которая пытается найти область, где сглаживающий сплайн наиболее чувствителен к параметру сглаживания (аналогично предыдущему обсуждению).

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 и SPAPS

csaps и spaps команды различаются тем, как задается определенный сглаживающий сплайн с помощью параметра сглаживания и допуска. Другая разница в том, что spaps может обеспечить линейный или квинтичный сглаживающий сплайн в дополнение к кубическому сглаживающему сплайну.

Квинтичный сглаживающий сплайн лучше кубического сглаживающего сплайна в ситуации, когда требуется, чтобы вторая производная перемещалась как можно меньше.