28.5 Полиномиальная интерполяция

Octave поддерживает различные виды интерполяции, большинство из которых описаны в Интерполяции. Простым вариантом функций, описанных в вышеупомянутой главе, является подгонка одного полинома или кускового полинома (сплайна) к заданным точкам данных. Чтобы избежать сильного колебания полинома, часто требуется подгонка полинома низкой степени к данным. Обычно это означает, что необходимо подстроить полином в смысле наименьших квадратов, что и делает функция polyfit.

: p = polyfit (x, y, n)

: [p, s] = polyfit (x, y, n)

: [p, s, mu] = polyfit (x, y, n)

Возвращает коэффициенты полинома p(x) степени n, минимизирующие ошибку наименьших квадратов подгонки к точкам [x, y].

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

Коэффициенты полинома возвращаются в строчном векторе.

Необязательный вывод s — это структура, содержащая следующие поля:

‘R’

Треугольный множитель R из QR-разложения.

‘X’

Матрица Вандермонда, используемая для вычисления коэффициентов полинома.

‘C’

Немасштабированная матрица ковариации, формально равная обратной матрице x’*x, но вычисляемая таким образом, чтобы минимизировать распространение ошибки округления.

‘df’

Степени свободы.

‘normr’

Норма остатков.

‘yf’

Значения полинома для каждого значения x.

Второй выходной параметр может быть использован функцией polyval для вычисления статистических пределов погрешностей предсказанных значений. В частности, стандартное отклонение коэффициентов p определяется выражением

sqrt (diag (s.C)/s.df)*s.normr.

Когда присутствует третий выходной параметр mu, коэффициенты p относятся к полиному в

xhat = (x - mu(1)) / mu(2)
где mu(1) = среднее (x), а mu(2) = стандартное отклонение (x).

Это линейное преобразование x улучшает числовую устойчивость подгонки.

См. также: polyval, polyaffine, roots, vander, zscore.

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

: pp = splinefit (x, y, breaks)

: pp = splinefit (x, y, p)

: pp = splinefit (…, "periodic", periodic)

: pp = splinefit (…, "robust", robust)

: pp = splinefit (…, "beta", beta)

: pp = splinefit (…, "order", order)

: pp = splinefit (…, "constraints", constraints)

Подгоняет кусковой кубический сплайн с разрывами (узлами) breaks к шумным данным x и y.

x — вектор, а y — вектор или N-мерный массив. Если y является N-мерным массивом, то x(j) сопоставляется с y(:,…,:,j).

p — положительное целое число, определяющее число интервалов вдоль x, а p+1 — число разрывов. Количество точек в каждом интервале отличается не более чем на 1.

Необязательное свойство periodic — логическое значение, указывающее, применяется ли периодическое граничное условие к сплайну. Длина периода составляет max (breaks) - min (breaks). Значение по умолчанию равно false.

Необязательное свойство robust — логическое значение, указывающее, применяется ли робастная подгонка для уменьшения влияния выбросов. Выполняется три итерации взвешенных наименьших квадратов. Веса вычисляются из предыдущих остатков. Чувствительность к идентификации выбросов регулируется свойством beta. Значение beta ограничено диапазоном 0 < beta < 1. Значение по умолчанию — beta = 1/2. Значения, близкие к 0, придают всем данным одинаковый вес. Увеличение значений beta уменьшает влияние выбросов. Значения, близкие к единице, могут вызвать нестабильность или недостаточность ранга.

Подгоняемый сплайн возвращается как кусковой полином pp и может быть оценен с помощью ppval.

Сплайны построены из полиномов степени order. По умолчанию — кубический, order=3. Сплайн с P кусками имеет P+order степеней свободы. При периодических граничных условиях степени свободы сокращаются до P.

Необязательное свойство constraints — структура, определяющая линейные ограничения на подгонку. Структура имеет три поля, "xc", "yc", и "cc".

"xc"

Вектор расположения ограничений x.

"yc"

Значения ограничений в местах xc. По умолчанию — массив нулей.

"cc"

Коэффициенты (матрица). По умолчанию — массив единиц. Количество строк ограничено степенью кусковых полиномов order.

Ограничения являются линейными комбинациями производных порядка от 0 до order-1 в соответствии с

cc(1,j) * y(xc(j)) + cc(2,j) * y'(xc(j)) + ... = yc(:,...,:,j).

См. также: interp1, unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps.

Число breaks (или узлов), используемых для построения кускового полинома, является важным фактором для подавления шума в входных данных x и y. Это демонстрируется примером ниже.

x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
## Uniform breaks
breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces
pp1 = splinefit (x, y, breaks);
## Breaks interpolated from data
pp2 = splinefit (x, y, 10);  % 11 breaks, 10 pieces
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"})

Результат которого показан на рисунке 28.1.

Рисунок 28.1: Сравнение подгонки кускового полинома с 41 разрывом с подгонкой с 11 разрывами. Подгонка с большим числом разрывов демонстрирует быстрые колебания, отсутствующие в исходной функции.

Подгонка кускового полинома, предоставляемая функцией splinefit, имеет непрерывные производные до order-1. Например, кубическая подгонка имеет непрерывные первую и вторую производные. Это демонстрируется кодом

## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
## Piecewise constant
pp1 = splinefit (x, y, 8, "order", 0);
## Piecewise linear
pp2 = splinefit (x, y, 8, "order", 1);
## Piecewise quadratic
pp3 = splinefit (x, y, 8, "order", 2);
## Piecewise cubic
pp4 = splinefit (x, y, 8, "order", 3);
## Piecewise quartic
pp5 = splinefit (x, y, 8, "order", 4);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
y4 = ppval (pp4, xx);
y5 = ppval (pp5, xx);
plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
axis tight
ylim auto
legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})

Результат которого показан на рисунке 28.2.

Рисунок 28.2: Сравнение кусковых постоянных, линейных, квадратичных, кубических и квартатных полиномов с 8 разрывами со случайными данными. Решения более высокой степени более точно представляют исходную функцию, но сопряжены с повышенной вычислительной сложностью.

Если подгонка должна быть к периодической функции, splinefit способна применить граничные условия, необходимые для получения периодической подгонки. Это демонстрируется кодом ниже.

## Data (100 points)
x = 2 * pi * [0, (rand (1, 98)), 1];
y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
## No constraints
pp1 = splinefit (x, y, 10, "order", 5);
## Periodic boundaries
pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "no constraints", "periodic"})

Результат которого показан на рисунке 28.3.

Рисунок 28.3: Сравнение подгонки кускового полинома к шумной периодической функции с периодическими и без периодических граничных условий.

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

## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (2 * x) + 0.1 * randn (size (x));
## Breaks
breaks = linspace (0, 2 * pi, 10);
## Clamped endpoints, y = y' = 0
xc = [0, 0, 2*pi, 2*pi];
cc = [(eye (2)), (eye (2))];
con = struct ("xc", xc, "cc", cc);
pp1 = splinefit (x, y, breaks, "constraints", con);
## Hinged periodic endpoints, y = 0
con = struct ("xc", 0);
pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "clamped", "hinged periodic"})

Результат которого показан на рисунке 28.4.

Рисунок 28.4: Сравнение двух периодических кусковых кубических подгонок к шумному периодическому сигналу. Одна подгонка имеет зафиксированные концы, а вторая — шарнирно закрепленные.

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

## Data
x = linspace (0, 2*pi, 200);
y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
## Add outliers
x = [x, linspace(0,2*pi,60)];
y = [y, -ones(1,60)];
## Fit splines with hinged conditions
con = struct ("xc", [0, 2*pi]);
## Robust fitting, beta = 0.25
pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
## Robust fitting, beta = 0.75
pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
## No robust fitting
pp3 = splinefit (x, y, 8, "constraints", con);
## Plot
xx = linspace (0, 2*pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
plot (x, y, ".", xx, [y1; y2; y3])
legend ({"data with outliers","robust, beta = 0.25", ...
         "robust, beta = 0.75", "no robust fitting"})
axis tight
ylim auto

Результат которого показан на рисунке 28.5.

Рисунок 28.5: Сравнение двух уровней робастной подгонки (beta = 0,25 и 0,75) к шумным данным с выбросами. Также включена обычная подгонка без робастной подгонки (beta = 0).

Очень специфической формой полиномиальной интерпретации является аппроксимация Паде. Для систем управления непрерывную задержку можно очень просто смоделировать с помощью аппроксиманта.

: [num, den] = padecoef (T)

: [num, den] = padecoef (T, N)

Вычислить паде-приближение N-го порядка для задержки в непрерывном времени T в форме передаточной функции.

Паде-приближение для exp (-sT) определяется следующим уравнением

Pn(s)
exp (-sT) ~ -------
             Qn(s)

Где оба Pn(s) и Qn(s) являются рациональными функциями N-го порядка, определёнными следующими выражениями

N    (2N - k)!N!        k
Pn(s) = SUM --------------- (-sT)
        k=0 (2N)!k!(N - k)!

Qn(s) = Pn(-s)

Вводные данные T и N должны быть неотрицательными числовыми скалярами. Если N не указано, оно по умолчанию равно 1.

Выходные векторы-строки num и den содержат коэффициенты числителя и знаменателя в порядке убывания степеней s. Оба являются многочленами N-го порядка.

Например:

t = 0.1;
n = 4;
[num, den] = padecoef (t, n)
⇒ num =

      1.0000e-04  -2.0000e-02   1.8000e+00  -8.4000e+01   1.6800e+03

⇒ den =

      1.0000e-04   2.0000e-02   1.8000e+00   8.4000e+01   1.6800e+03

Функция ppval, оценивает куски многочленов, созданные с помощью mkpp или другими средствами, и unmkpp возвращает подробную информацию о кусковом многочлене.

Следующий пример показывает, как объединить две линейные функции и квадратную в одну функцию. Каждая из этих функций выражена на смежных интервалах.

x = [-2, -1, 1, 2];
p = [ 0,  1, 0;
      1, -2, 1;
      0, -1, 1 ];
pp = mkpp (x, p);
xi = linspace (-2, 2, 50);
yi = ppval (pp, xi);
plot (xi, yi);

: pp = mkpp (breaks, coefs)

: pp = mkpp (breaks, coefs, d)

Создать структуру кускового многочлена (pp) из точек выборки breaks и коэффициентов coefs.

breaks должен быть вектором строго возрастающих значений. Количество интервалов задаётся ni = length (breaks) - 1.

Когда m это порядок многочлена, coefs должен иметь размер: ni-by-(m + 1).

i-я строка coefs, coefs(i,:), содержит коэффициенты для многочлена на i-м интервале, упорядоченные от старшей (m) до младшей (0) степени.

coefs также может быть многомерным массивом, определяющим векторнозначный или многомернозначный многочлен. В этом случае порядок многочлена m определяется длиной последнего измерения coefs. Размер первых измерений задаётся скаляром или вектором d. Если d не указан, он устанавливается в 1. В этом случае p(r, i, :) содержит коэффициенты для r-го многочлена, определённого на интервале i. В любом случае coefs преобразуется в матрицу 2-х мерного размера [ni*prod(d) m].

Примечание для программистов: ppval оценивает многочлены в xi - breaks(i), т.е. он вычитает нижнюю границу текущего интервала из xi. Это необходимо учитывать при создании объектов кусковых многочленов с помощью mkpp.

См. также: unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps.

: [x, p, n, k, d] = unmkpp (pp)

Извлечь компоненты структуры кускового многочлена pp.

Эта функция является обратной к mkpp: она извлекает входные данные для mkpp необходимые для создания структуры кускового многочлена PP. Код ниже делает это отношение явным:

[breaks, coefs, numinter, order, dim] = unmkpp (pp);
pp2  = mkpp (breaks, coefs, dim);

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

Компоненты:

x

Точки выборки или точки разрыва.

p

Коэффициенты многочлена для точек в интервале выборки. p(i, :) содержит коэффициенты для многочлена на интервале i, упорядоченные от старшей до младшей степени. Если d > 1, то p является матрицей размера [n*prod(d) m], где i + (1:d) строки - коэффициенты всех d многочленов в интервале i.

n

Количество кусков многочлена или интервалов, n = length (x) - 1.

k

Порядок многочлена плюс 1.

d

Количество многочленов, определённых для каждого интервала.

См. также: mkpp, ppval, spline, pchip.

: yi = ppval (pp, xi)

Оценить структуру кускового многочлена pp в точках xi.

Если pp описывает скалярную функцию многочлена, результат является массивом того же размера, что и xi. В противном случае размер результата равен [pp.dim, length(xi)] если xi является вектором, или [pp.dim, size(xi)] если это многомерный массив.

См. также: mkpp, unmkpp, spline, pchip.

: ppd = ppder (pp)

: ppd = ppder (pp, m)

Вычислить кусковую m-ую производную структуры кускового многочлена pp.

Если m опущено, вычисляется первая производная.

См. также: mkpp, ppval, ppint.

: ppi = ppint (pp)

: ppi = ppint (pp, c)

Вычислить интеграл структуры кускового многочлена pp.

c, если указано, является константой интегрирования.

См. также: mkpp, ppval, ppder.

: jumps = ppjumps (pp)

Оценить скачки на границах кускового многочлена.

Если есть n интервалов и размерность pp равна d, результирующий массив имеет размерность [d, n-1].

См. также: mkpp.