Интерполяция, экстраполяция, сплайны, регрессия

На практике часто возникает необходимость построения для заданной функции y=f(x) приближенных формул, когда нужно подобрать некоторую функцию φ(x), которая близка к y=f(x) и просто вычисляется. Близость φ(x) к y=f(x) достигается введением в аппроксимирующую функцию φ(x) свободных параметров и соответ­ствующим их выбором. Подбор удачного вида функциональной зависимости φ(x,a) – искусство.

Пусть функция y=f(x) известна только в узлах некоторой сетки x_i, т.е. задана табли­цей. Если потребовать, чтобы φ(x,a) совпадала с табличными значениями в n выбранных узлах сетки, то получим систему

,

из которой можно определить параметры a_k. Такой способ вы­бора параметров называется лагранжевой интерполя­цией. Термин интерполяция в узком смысле употребляют, если x заключено между крайними узлами интерполяции. Если значение x выходит из этих пределов, то говорят об экстраполяции.

Для построения интерполяции/экстраполяции в Mathcad имеются несколько встроенных функций, позволяющих соединить точки данных кривыми разной степени гладкости.

Самый простой вид интерполяции – линейная, которая представляет искомую зависимость A(x) в виде ломаной линии, т.е. интерполирующая функция A(x) состоит из отрезков прямых, соединяющих точки.

Для построения линейной интерполяции служит встроенная функция linterp(x,y,t), где x – вектор значений аргумента, y – вектор (того же размера) значений функции, t – значение аргумента, при котором вычисляется интерполирующая функция. Элементы вектора x должны быть определены в порядке возрастания, т.е. .

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

,

коэффициенты которого кусочно-постоянны и который в узлах интерполяции принимает заданные значения и непрерывен вместе со своими n-1 производными, т.е.

здесь N – число узлов интерполяции.

В Mathcad сплайн-интерполяция реализована при помощи встроенной функции interp(s,x,y,t), где x – вектор значений аргумента, элементы которого расположены в по­ряд­ке возрастания; y – вектор значений функции того же размера; t – значение аргу­мен­та, при котором вычисляется интерполирующая функция, а s – вектор коэффи­циентов сплай­на, созданный при помощи одной из сопутствующих встроенных функций тех же аргументов (x,y):

• lspline(x,y) – для линейного сплайна;
• pspline(x,y) – для квадратичного сплайна;
• cspline(x,y) – для кубического сплайна.

На практике наиболее часто используется кубическая сплайн-интерполяция.

Пример кубической сплайн-интерполяции.

Более сложный тип интерполяции – интерполяция B-сплайнами. В отличие от обычной сплайн-интерполяции, сшивка элементарных B-сплайнов производится не в узлах x_i, а в точках u_i, координаты которых нужно ввести пользователю. B-сплайны могут быть полиномами 1, 2 или 3 степени (линейные, квадратичные или кубические).

В Mathcad интерполяция B-сплайнами реализована точно так же, как и обычная сплайн-интерполяция, различие состоит в определении коэффициентов сплайна, которая осуществляется при помощи встроенной функции bspline(x,y,u,n), где x,y – векторы значений аргумента и функции в узлах; u – вектор значений аргумента, в которых производится сшивка B-сплайнов; n – порядок полиномов сплайновой интерполяции (1, 2 или 3).

Размер вектора u должен быть на 1, 2 или 3 меньше размера векторов x и y. Первый элемент вектора u должен быть меньше или равен первому элементу вектора x (, т.е. лежать за пределами левой границы интервала интерполирования). Последний элемент вектора u должен быть больше или равен последнему элементу вектора x (т.е. лежать за пределами правой границы интервала интерполирования).

Интерполяция позволяет легко аппроксимировать функцию y=f(x). Однако точ­ность такой аппроксимации гарантирована лишь в небольшом интервале порядка не­сколь­ких шагов сетки. Для другого интервала приходится заново вычислять коэффици­ен­ты интерполяционной формулы. Удобнее иметь единую приближенную формулу, пригодную для большого диапазона изменения аргумента. При интерполяции мы приравниваем значения функций f(x) и φ(x) в узлах x_i. Если значения y_i=f(x_i) определены неточно – например, из эксперимента, – то точное приравнивание неразумно. В этом случае целесо­об­разнее приближать функцию не по точкам, а в среднем.

Задачи математической регрессии имеют смысл приближения данных некоторой функцией φ(x,a₁,…,a_k), минимизирующей функционал

,

т.е. решается задача .

Таким образом, регрессия сводится к подбору неизвестных коэффициентов, опре­де­ляю­щих аналитическую зависимость φ(x,a₁,…,a_k). Как правило, регрессия очень эффективна, если заранее известен или хорошо угадывается закон распределения данных.

Самый простой и наиболее часто используемый вид регрессии – линейная. В этом случае приближение данных осуществляется линейной функцией φ(x)=b+ax, которая на координатной плоскости изображается прямой линией. Линейную регрессию часто назы­вают методом наименьших квадратов.

Для расчета линейной регрессии b+ax в Mathcad имеются встроенные функции, реализующие два дублирующих друг друга способа:

• Функция line(x,y) – возвращает вектор (b,a) коэффициентов линейной регрессии.
• Функция intercept(x,y) – возвращает коэффициент b линейной регрессии, а функция slope(x,y) – возвращает коэффициент a линейной регрессии.

Здесь x, y – векторы данных одинакового размера.

Кроме того, в Mathcad имеется альтернативный алгоритм, реализующий медиан-медианную линейную регрессию для расчетов коэффициентов a и b через встроенную векторную функцию medfit(x,y).

В Mathcad реализована регрессия одним полиномом, отрезками нескольких полиномов и двумерная регрессия массива данных, причем она осуществляется в комбинации со встроенной функцией полиномиальной интерполяции interp(s,x,y,t) (рассмотренной нами ранее). Для этого используются следующие встроенные функции

• Функция regress(x,y,k) – возвращает вектор коэффициентов для построения полиномиальной регрессии данных.
• Функция loess(x,y,span) – возвращает вектор коэффициентов для построения регрессии данных фрагментами полиномов.

Здесь x – вектор значений аргумента, элементы которого расположены в порядке возрастания; y – вектор значений функции того же размера; k – степень полинома регрессии (целое положительное число); span – положительный параметр сглаженности данных, определяющий размер фрагментов полиномов (хорошие результаты дает значение span≈0.75).

Для построения регрессии полиномом степени k необходимо наличие, по крайней мере, (k+1)-ой точек данных.

Пример. Полиномиальная регрессия квадратичной параболой.

x:=(0 1 2 3 4 5 6)^T

y:=(4.1 2.1 3.2 3.6 4.3 5.2 6.0)^T

k:=2

s:=regress(x,y,k)

A(t):=interp(s,x,y,t)

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

Кроме рассмотренных, в Mathcad встроено еще несколько видов регрессии специального вида. Это трехпарамет­ри­чес­кой регрессии (с искомыми параметрами a, b, c), и регрессия в виде линейной комбинации C₁f₁+C₂f₂+…, где f_i(x) – любые функции пользователя, а C_i – подлежащие определению коэффициенты. К тому же, имеется путь проведения регрессии общего вида, когда комбинацию функций и искомых коэффициен­тов задает сам пользователь.

Для построения специальной регрессии используются следующие встроенные функции:

• expfit(x,y,g) – регрессия экспонентой f(x)=a·e^bx+c;
• lgsfit(x,y,g) – регрессия логистической функцией f(x)=a/(1+b·e^-cx);
• sinfit(x,y,g) – регрессия синусоидой f(x)=a·sin(x+b)+c;
• pwfit(x,y,g) – регрессия степенной функцией f(x)=a·x^b+c;
• logfit(x,y,g) – регрессия логарифмической функцией f(x)=a·ln(x+b)+c;
• lnfit(x,y) – регрессия двухпараметрической логарифмической функцией f(x)=a·ln(x)+b;
• linfit(x,y,F) – регрессия в виде линейной комбинации функций пользователя;
• genfit(x,y,g,G) – регрессия в виде функции пользователя общего вида.

Здесь x – вектор значений аргумента данных, y – вектор значений функции данных, g – вектор, задающий начальные значения параметров регрессии, F – пользовательская векторная функция скалярного аргумента, G – векторная (N+1)-мерная функция, составленная из функции пользователя и ее N частных производных по каждому из параметров.

Впрочем, начиная с Mathcad 13, для функции регрессии общего вида genfit доста­точ­но определить в качестве последнего аргумента только саму функцию регрессии F, не тратя время на ввод ее производных.

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

Пример. Расчет экспоненциальной регрессии

x:=(0 1 2 3 4 5 6)^T

y:=(4.1 2.1 3.2 3.6 4.3 5.2 6.0)^T

g:=(1 1 1)^T

C:=expfit(x,y,g)

F(t):=C₀·exp(C₁·t)+C₂

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

Виктор Семенович Болдасов

Задания практики по Mathcad

и методические указания по их выполнению

В авторской редакции

Тематический план 2007 г., позиция

Подписано в печать Формат 60 ´ 84 ¤ 16

Бумага Печать офсетная.

Усл.печ.л. Уч.-изд.л.

Тираж экз. Заказ

Издательство МГУП.

Отпечатано в ИПК МГУП.

127550, Москва, ул. Прянишникова, 2а.