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 обычно является целым числом ≥ 0, определяющим степень аппроксимирующего полинома. Если n является логическим вектором, он используется в качестве маски для выборочного включения или исключения соответствующих коэффициентов полинома.

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

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

‘yf’

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

‘X’

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

‘R’

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

‘C’

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

‘df’

Число степеней свободы.

‘normr’

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

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

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

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

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

Пример 1: логическое n и целое n

f = @(x) x.^2 + 5;   # data-generating function
x = 0:5;
y = f (x);
## Fit data to polynomial A*x^3 + B*x^1
p = polyfit (x, y, logical ([1, 0, 1, 0]))
⇒ p = [ 0.0680, 0, 4.2444, 0 ]
## Fit data to polynomial using all terms up to x^3
p = polyfit (x, y, 3)
⇒ p = [ -4.9608e-17, 1.0000e+00, -3.3813e-15, 5.0000e+00 ]

См. также: 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.

Подгонка кусочно-полиномиальной функции, предоставленная 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.

Когда подлежащая функция для предоставления подгонки является периодической, 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.

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

## 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.

Функция 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.

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

: [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.