Кубические сплайны сглаживания

Этот пример показывает, как использовать csaps и команды spaps от Curve Fitting Toolbox™, чтобы создать кубические сплайны сглаживания.

Команда 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, мы получаем interpolant к (шумным) данным.

На самом деле формулировка, используемая csapi (p.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 может обеспечить линейное или quintic, сглаживающий сплайн, в дополнение к кубическому сплайну сглаживания.

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

Для просмотра документации необходимо авторизоваться на сайте