Конструкция аппроксимации методом наименьших квадратов обычно требует, чтобы один имел в руках базиса для пространства, из которого должны быть аппроксимированы данные. Как иллюстрирует пример пространства «естественных» кубических сплайнов, явную конструкцию базиса не всегда простое.
В этом разделе четко указывается, что явный базис на самом деле не нужен; достаточно иметь некоторые доступные средства интерполяции в некоторой степени из пространства аппроксимаций. Для этого необходим тот факт, что функции Curve Fitting Toolbox™ spline поддерживают работу с векторно значимыми функциями.
В этом разделе рассматриваются эти аспекты приближения методом наименьших квадратов «естественными» кубическими сплайнами.
Вы хотите создать приближение методом наименьших квадратов к данным (x, y) из пространственного S «естественных» кубических сплайнов с заданными пропусками b(1) <... <b(l+1).
Если вы знаете базис, (f1, f2,..., fm), для линейного пространства S всех «естественных» кубических сплайнов с пропуском последовательностью b, тогда вы научились находить приближение методом наименьших квадратов в форме c(1)f1 + c(2)f2 +... + c(m)fm, с вектором 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 <reservedrangesplaceholder0>f2 <reservedrangesplaceholder0>... +c(m)fm. Ибо, при этом, можно получить, для j=1:m, а j-й столбец A как fnval(F(ej),x), с ej а j-й столбец eye(m), матрица тождеств порядка m.
Еще лучше, что функции сплайна Curve Fitting Toolbox могут обрабатывать векторные функции, поэтому вы должны иметь возможность создавать карту базиса F для обработки векторных коэффициентов c(i) также. Однако, по соглашению, в этом тулбоксе векторный коэффициент является вектором-столбцом, поэтому последовательность c обязательно является вектором-строкой векторов-столбцов, то есть матрицей. При этом F(eye(m)) - векторный сплайн, чей i-й компонент является базисным элементом f i, i=1:m. Следовательно, принимая вектор x сайтов данных, которые будут вектором-строкой, fnval(F(eye(m)),x) - матрица, чья (i,j)-entry - значение f i при x(j), т.е. транспонирование матрицы A вы ищете. С другой стороны, как только что было указано, ваша базовая карта F ожидает последовательность коэффициентов c быть вектором-строкой, т.е. транспонировать вектор A\y. Следовательно, принимая, соответственно, вектор y из значений данных, которые будут вектором-строкой, можно получить приближение наименьших квадратов от S к данным (x, y) как
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
Посмотрите на приближение методом «Natural» Cubic Splines With Three Interior Breaks, которая показывает, в густом синем цвете, полученное приближение, наряду с заданными данными.
Это выглядит как хорошее приближение, за исключением того, что это не похоже на «естественный» кубический сплайн. «Естественный» кубический сплайн, чтобы вспомнить, должен быть линейным слева от своего первого пропуска и справа от последнего пропуска, и это приближение не удовлетворяет ни одному из условий. Это связано со следующими фактами.
«Естественная» кубическая сплайн интерполяция к данным обеспечивается csape в pform, с интервалом, охватываемым сайтами данных, его основной интервал. С другой стороны, оценка pporm вне его основного интервала выполняется в MATLAB ppval или функцию сплайна Curve Fitting Toolbox fnvalпри помощи соответствующей полиномиальной концевой части ppform, то есть путем экстраполяции в полном порядке. В случае «естественного» кубического сплайна, вы хотите вместо этого экстраполяцию второго порядка. Это означает, что вы хотите, слева от первого разрыва, прямой пропуск, которая согласуется с кубическим сплайном по значению и наклон на первом пропуске. Такая экстраполяция обеспечивается fnxtr. Поскольку «естественный» кубический сплайн имеет нулевую вторую производную при первом пропуске, такая экстраполяция является даже третьим порядком, то есть удовлетворяет трем совпадающим условиям. Точно так же, за пределами последнего пропуска кубического сплайна, требуется прямая линия, которая согласуется со сплайном в значении и наклона на последнем пропуске, и это тоже обеспечивается fnxtr.
Наименьшие квадраты, Приближение «натуральными» кубическими сплайнами с тремя внутренними Пропусками
Следующий однострочный код обеспечивает правильное приближение методом наименьших квадратов к данным (x, y) «естественными» кубическими сплайнами с последовательностью пропусков b:
fnxtr(csape(b,y(:).'/ ...
fnval(fnxtr(csape(b,eye(length(b)),'var')),x(:).'),'var'))
Но это, надо признать, довольно длинная линия.
Следующий код использует эту правильную формулу и строит более тонкую, красную линию, получая приближение поверх более ранних графиков, как показано на приближении методом наименьших квадратов «Natural» Cubic Splines With Three Interior Break.
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, потому что можно полагаться на MATLAB spline. Ты знаешь это, с 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 к данным (x, y) :
spline(b,y(:)'/spline(b,eye(length(b)),x(:)'))