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 преобразуется в двумерную матрицу размера [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.