Приближение наименьших квадратов естественными кубическими сплайнами

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

Этот раздел ясно дает понять, что явный базис не на самом деле необходим; достаточно иметь доступный некоторые средние значения интерполяции некоторым способом от пробела аппроксимирующих функций. Для этого то, что поддержка функций сплайна Curve Fitting Toolbox™ работает с функциями с векторным знаком, важно.

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

Проблема

Вы хотите создать приближение наименьших квадратов к определенным данным (xY) от пробела S “естественных” кубических сплайнов с данными пропусками b(1) <... <b(l+1).

Общее разрешение

Если вы знаете базис, (f1, f2..., из), для линейного пробела S всех “естественных” кубических сплайнов с последовательностью пропуска b, затем вы учились находить приближение наименьших квадратов в форме c(1)f1 + c(2)f2 +... + c(m)из, с векторным c решение методом наименьших квадратов к линейной системе A*c = y, чьей матрицей коэффициентов дают

A(i,j) = fj(x(i)),   i=1:length(x),  j=1:m .

Другими словами, c = A\y.

Потребность в базисной карте

Общее решение, кажется, требует, чтобы вы знали базис. Однако для того, чтобы создать содействующую последовательность c, только необходимо знать матричный A. Для этого достаточно иметь под рукой базисную карту, а именно, функциональный F скажите, так, чтобы F(c) возвращает сплайн, данный конкретной взвешенной суммой c(1)f1+c(2)f2+... +c(m)из. Поскольку, с этим можно получить для j=1:m, j- столбец th A как fnval(F(ej),x), с ej j- столбец th eye(m), единичная матрица порядка m.

Еще лучше функции сплайна Curve Fitting Toolbox могут обработать функции с векторным знаком, таким образом, необходимо смочь создать базисную карту F обрабатывать коэффициенты с векторным знаком c(i) также. Однако по соглашению, в этом тулбоксе, коэффициент с векторным знаком является вектор-столбцом, следовательно последовательность c является обязательно вектором-строкой из вектор-столбцов, т.е. матрицей. С этим, F(eye(m)) сплайн с векторным знаком чей i- компонент th является базисным элементом fi, i=1:m. Следовательно, принимая векторный x из сайтов данных, чтобы быть вектором-строкой, fnval(F(eye(m)),x) матрица чей (i,j)- запись является значением fi в x(j), т.е. транспонирование матричного A вы ищете. С другой стороны, как только указано, ваша базисная карта F ожидает содействующую последовательность c быть вектором-строкой, т.е. транспонированием векторного A\y. Следовательно, принятие, соответственно, векторный y из значений данных, чтобы быть вектором-строкой, можно получить приближение наименьших квадратов от S до данных (xY) как

F(y/fnval(F(eye(m)),x))

Безусловно, если вы хотели быть подготовленными к x и y чтобы быть произвольными векторами (той же длины), вы использовали бы вместо этого

F(y(:).'/fnval(F(eye(m)),x(:).'))

Базисная карта для “естественных” кубических сплайнов

Что точно требуется базисной карты F для линейного пробела S “естественных” кубических сплайнов с последовательностью пропуска b(1) < ... < b(l+1)? Принятием размерности этого линейного пробела является m, карта F должен настроить линейное взаимно-однозначное соответствие между m- векторы и элементы S. Но это точно что csape(b, . ,'var') делает.

Чтобы быть явными, рассмотрите следующий функциональный F:

function s = F(c)
s = csape(b,c,'var');

Для данного векторного c (той же длины как b), это предоставляет уникальному “естественному” кубическому сплайну последовательность пропуска b, который принимает значение c(i) в b(i), i=1:l+1. Уникальность является ключевой. Это гарантирует что соответствие между векторным c и получившийся сплайн F(c) является непосредственным. В частности, m равняется length(b). Больше, чем это, потому что значение f (t) функционального f в точке t зависит линейно от f, эта уникальность, гарантируют тот F(c) зависит линейно от c (потому что c равняется fnval(F(c),b) и инверсия обратимой линейной карты является снова линейной картой).

Короткое решение

Соединяя все это, вы прибываете в следующий код

csape(b,y(:).'/fnval(csape(b,eye(length(b)),'var'),x(:).'),...
'var')

для приближения наименьших квадратов “естественными” кубическими сплайнами с последовательностью пропуска b.

Потребность в соответствующей экстраполяции

Давайте попробуем его на некоторых данных, данные о переписи, скажем, которые обеспечиваются в MATLAB® командой

load census

и который предоставляет годы, 1790:10:1990, как cdate и значения как pop. Используйте последовательность пропуска 1810:40:1970.

b = 1810:40:1970;
s = csape(b, ...
pop(:)'/fnval(csape(b,eye(length(b)),'var'),cdate(:)'),'var');
fnplt(s, [1750,2050],2.2);
hold on
plot(cdate,pop,'or');
hold off

Взгляните на Приближение Наименьших квадратов “Естественными” Кубическими Сплайнами С Тремя Внутренними Пропусками, которое показывает, в толстом синем, получившемся приближении, наряду с определенными данными.

Это похоже на хорошее приближение, - за исключением того, что оно не похоже на “естественный” кубический сплайн. “Естественный” кубический сплайн, чтобы вспомнить, должен быть линейным слева от его первого пропуска и справа от его последнего пропуска, и это приближение не удовлетворяет никакому условию. Это происходит из-за следующих фактов.

“Естественный” кубический сплайн interpolant к определенным данным обеспечивается csape в ppform, с интервалом, заполненным по условию, располагает его основной интервал. С другой стороны, оценка ppform вне ее основного интервала сделана в MATLAB ppval или Curve Fitting Toolbox шлицует функциональный fnval, при помощи соответствующей полиномиальной части конца ppform, т.е. экстраполяцией полного порядка. В случае “естественного” кубического сплайна вы хотите вместо этого экстраполяцию второго порядка. Это означает, что вы хотите, слева от первого пропуска, прямая линия, которая соглашается с кубическим сплайном в значении и наклоне в первом пропуске. Такая экстраполяция обеспечивается fnxtr. Поскольку “естественный” кубический сплайн имеет нулевую вторую производную в своем первом пропуске, такая экстраполяция является даже третьим порядком, т.е. это удовлетворяет трем соответствующим условиям. Таким же образом, вне последнего пропуска кубического сплайна, вы хотите прямую линию, которая соглашается со сплайном в значении и наклоне в последнем пропуске, и это также предоставляется fnxtr.

Приближение наименьших квадратов “естественными” кубическими сплайнами с тремя внутренними пропусками

Правильное короткое решение

Следующий короткий код предоставляет правильное приближение наименьших квадратов данным (xY) “естественными” кубическими сплайнами с последовательностью пропуска b:

fnxtr(csape(b,y(:).'/ ...
    fnval(fnxtr(csape(b,eye(length(b)),'var')),x(:).'),'var'))

Но это - по общему признанию, довольно длинная линия.

Следующий код использует эту правильную формулу и графики, в разбавителе, красной линии, получившемся приближении сверху более ранних графиков, как показано в Приближении Наименьших квадратов “Естественными” Кубическими Сплайнами С Тремя Внутренними Пропусками.

 ss = fnxtr(csape(b,pop(:)'/ ...
  fnval(fnxtr(csape(b,eye(length(b)),'var')),cdate(:)'),'var'));
hold on, fnplt(ss,[1750,2050],1.2,'r'),grid, hold off
legend('incorrect approximation','population', ...
'correct approximation')

Приближение наименьших квадратов кубическими сплайнами

Короткое решение работает отлично, если вы хотите аппроксимировать пробелом S всех кубических сплайнов с данной последовательностью пропуска b. Вы не должны даже использовать функции сплайна Curve Fitting Toolbox для этого, потому что можно использовать spline MATLAB. Вы знаете это с c последовательность, содержащая еще две записи, чем, делает b, spline(b,c) предоставляет уникальному кубическому сплайну последовательность пропуска b это принимает значение c(i+1) в b(i), весь i, и берет наклонный c(1) в b(1), и наклонный c(end) в b(end). Следовательно, spline(b,.) базисная карта для S.

Больше, чем это, вы знаете тот spline(b,c,xi) вводит значение в xi из этого сплайна интерполяции. Наконец, вы знаете тот spline может обработать данные с векторным знаком. Поэтому следующий короткий код создает приближение наименьших квадратов кубическими сплайнами с последовательностью пропуска b к данным (xY) :

spline(b,y(:)'/spline(b,eye(length(b)),x(:)'))