Paired piecewise linear regression model with fixed nodes

Full Text

ВВЕДЕНИЕ

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

Исходными данными парной регрессионной модели являются наборы значений объясняющей и зависимой переменных, которые геометрически изображаются точками координатной плоскости, образуя корреляционное поле.

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

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

На заключительном этапе моделирования проверяется адекватность построенной модели, т. е. соответствует ли построенная модель априорным данным.

Постмодельный этап использует основные предпосылки регрессионного анализа для проверки значимости параметров аппроксимирующей функции — достаточно ли объясняющей переменной для описания зависимой переменной.

Парная линейная регрессионная модель является самой распространенной ввиду простоты получения ее оценок как точечных, так и интервальных. Тем не менее реальные физические и социально-экономические процессы, как правило, имеют более сложную, нелинейную природу. Существуют априорные приемы выявления класса нелинейных функций (степенных, показательных, логарифмических и т. д.), к которому принадлежит аппроксимирующая функция. Парные нелинейные регрессионные модели, перечисленные выше, относятся к числу линеаризующихся моделей. Они надлежащей заменой переменных приводятся к линейной вспомогательной модели. Нахождение оценок параметров нелинейной модели, а также проверка значимости этих параметров осуществляется на основании вспомогательной линейной модели.

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

Более эффективный путь — использование кусочно-непрерывной аппроксимации с известными кусками. Если в качестве кусков брать линейные функции, то получатся кусочно-линейные модели. Они различаются методами построения аппроксимирующих кусочно-линейных функций.

В работе (Боровской, Костелей, 2017, с. 73) применяется алгоритм кусочно-линейной аппроксимации при прогнозировании финансовых временных рядов. Здесь строится набор трендовых линий, разделяющих область интервала и не соединенных между собой. Затем для получения непрерывной трендовой модели используется механизм сглаживания. В результате получается трендовая модель, состоящая из последовательно соединенных трендов.

В работах (Носков, 2013, с. 135; Иванова, Лебедева, Носков, 2016, с. 107; Носков, Лоншаков, 2008, с. 63; Носков, Хоняков, 2019, с. 47) рассмотрена кусочно-линейная регрессионная модель, называемая также производственной функцией Леонтьева, или функцией с нулевой эластичностью замены ресурсов. Идентификация параметров этой модели осуществлена при помощи метода наименьших модулей, реализация которого сводится к задаче линейно-булевого программирования.

В работе (Носков, 2023, с. 218) осуществляется формализация вложенной кусочно-линейной регрессии двух типов, являющихся комбинациями традиционной кусочно-линейной регрессии и функции риска. Нахождение оценок параметров вложенной кусочно-линейной регрессии сводится к решению задачи линейно-булевого программирования.

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

1. КЛАСС КУСОЧНО-ЛИНЕЙНЫХ ФУНКЦИЙ С ФИКСИРОВАННЫМИ УЗЛАМИ

Пусть имеются наборы ${\left(x_i\right)}_{i=\mathrm{0,}\mathrm{...,}n-1}$ и ${\left(y_i\right)}_{i=\mathrm{0,}\mathrm{...,}n-1}$ значений объясняющей $X$ и зависимой $Y$ переменных. Точки $\left(x_i,y_i\right)$, $i=\mathrm{0,}\mathrm{...,}n-1$ координатной плоскости $x\times y$ образуют корреляционное поле $\mathbf{K}$.

Будем предполагать, что последовательность ${\left(x_i\right)}_{i=\mathrm{0,}\mathrm{...,}n-1}$ строго возрастает, т. е. $x_0<x_1<\mathrm{...}<x_{n-1}$. При необходимости последовательности ${\left(x_i\right)}_{i=\mathrm{0,}\mathrm{...,}n-1}$ и ${\left(y_i\right)}_{i=\mathrm{0,}\mathrm{...,}n-1}$ можно перенумеровать так, чтобы условие возрастания выполнялось.

Для точек $u\text{'}$ и $u\text{'}\text{'}$ таких, что $\left\{u^{\text{'}},u^{\text{'}\text{'}}\right\}\cap \left\{x_i:i=\mathrm{0,}\mathrm{...,}n-1\right\}=\varnothing$ и $u^{\text{'}}<u^{\text{'}\text{'}}$, множество $K_{u^{\text{'}};u^{\text{'}\text{'}}}$ точек $\left(x_i,y_i\right)$, $i=\mathrm{0,}\mathrm{...,}n-1$ корреляционного поля $\mathbf{K}$, удовлетворяющих условию $u^{\text{'}}<x_i<u^{\text{'}\text{'}},$ будем называть полосовым участком, а точки $u^{\text{'}}$ и $u\text{'}\text{'}$ его узлами, левым и правым соответственно. Случай бесконечных узлов не исключается.

Узлы полосового участка $K_{u^{\text{'}};u^{\text{'}\text{'}}}$ определяются неоднозначно. Если $u^{\text{'}}<x_0$, то в качестве левого узла можно взять любое число интервала $\left(-\infty ;x_0\right)$ или $-\infty$. Если $x_0<u^{\text{'}}<x_{n-1}$, то найдется такое натуральное число $k$, что $0<k<n-1$ и $x_{k-1}<u^{\text{'}}<x_k$. В этом случае каждое число интервала $\left(x_{k-1};x_k\right)$ может быть взято в качестве левого узла. Если $x_{n-1}<u^{\text{'}\text{'}}$, то правым узлом может быть любое число интервала $\left(x_{n-1};+\infty \right)$ или $+\infty$. Если $x_0<u^{\text{'}\text{'}}<x_{n-1}$, то найдется такое натуральное число $l$, что $0<l\le n-1$ и $x_{l-1}<u^{\text{'}\text{'}}<x_l$. Тогда каждое число интервала $\left(x_{l-1};x_l\right)$ может быть взято в качестве правого узла. Вопрос выбора хороших узлов, т. е. узлов, удовлетворяющих нужным условиям (например, оптимальности), остается открытым. В общем случае, когда на узлы не накладываются условия, левым узлом будем брать $-\infty$ или середину интервала $\left(x_{k-1};x_k\right)$, а правым узлом $+\infty$ или середину интервала $\left(x_{l-1};x_l\right)$.

С геометрической точки зрения для конечных узлов $u\text{'}$ и $u\text{'}\text{'}$ прямые $x=u^{\text{'}}$ и $x=u^{\text{'}\text{'}}$ вырезают из корреляционного поля полосовой участок. Эти прямые будем называть узловыми.

Полосовым разбиением корреляционного поля будем называть конечную последовательность непустых попарно непересекающихся полосовых участков, объединение которых совпадает с корреляционным полем. При этом для каждого полосового участка этой последовательности, кроме последнего, его правый узел совпадает с левым узлом следующего участка данной последовательности. Элементы полосового разбиения будем называть частичными участками разбиения. Таким образом, последовательность $\mathbf{R}={\left(K_{{u^{\text{'}}}_j;{u^{\text{'}\text{'}}}_j}^{\left(j\right)}\right)}_{j=\mathrm{1,}\mathrm{...,}m}$ полосовых участков является полосовым разбиением, если:

1. $K_{{u^{\text{'}}}_j;{u^{\text{'}\text{'}}}_j}^{\left(j\right)}\ne \varnothing$ $\forall j=\mathrm{1,}\mathrm{...,}m$;
2. $K_{{u^{\text{'}}}_j;{u^{\text{'}\text{'}}}_j}^{\left(j\right)}\cap K_{{u^{\text{'}}}_l;{u^{\text{'}\text{'}}}_l}^{\left(l\right)}=\varnothing$ $\forall j=\mathrm{1,}\mathrm{...,}m$; $\forall l=\mathrm{1,}\mathrm{...,}m;$ $j\ne l$;
3. ${u^{\text{'}\text{'}}}_j={u^{\text{'}}}_{j+1}$ $\forall j=\mathrm{1,}\mathrm{...,}m-1.$

Если число частичных участков полосового разбиения равно $m$, то их попарно различные узлы, расположенные в порядке возрастания, обозначим как

$a$, $u^{\left(0\right)}$, $u^{\left(1\right)}$, …, $u^{(m-2)}$, $b$. (1)

Причем в общем случае, $a=-\infty$ и $b=+\infty$. Тогда $K_{a;u^{\left(0\right)}}$ — первый частичный участок разбиения; $K_{u^{\left(0\right)};u^{\left(1\right)}}$ — второй, …; $K_{u^{(m-2)};b}$ — участок m. Найдется такая возрастающая последовательность ${\left(s_j\right)}_{j=\mathrm{1,}\mathrm{...,}m-1}$ первых n — 1 натуральных чисел, что $x_{s_{j+1}-1}<u^{\left(j\right)}<x_{s_{j+1}}$ для каждого $j=\mathrm{0,}\mathrm{...,}m-1$. Для удобства введем обозначения: $s_0=0$, $s_m=n$. Тогда:

$K_{a;u^{\left(0\right)}}=\left\{\left(x_i,y_i\right):i=s_0,\mathrm{...,}{s}_1-1\right\}$;

$K_{u^{(j-1)};u^{\left(j\right)}}=\left\{\left(x_i,y_i\right):i=s_j,\mathrm{...,}{s}_{j+1}-1\right\}$, $j=\mathrm{1,}\mathrm{...,}m-2$;

$K_{u^{(m-2)};b}=\left\{\left(x_i,y_i\right):i=s_{m-1},\mathrm{...,}{s}_m-1\right\}$.

Геометрически прямые $x=u^{\left(0\right)}$, $x=u^{\left(1\right)}$, …, $x=u^{(m-2)}$ разрезают корреляционное поле на полосовые участки, составляющие его полосовое разбиение.

Полосовой участок $K_{u^{\text{'}};u^{\text{'}\text{'}}}$ будем называть линейным, если, рассматривая его как отдельное корреляционное поле, можно эффективно применять линейную аппроксимацию. Полосовое разбиение корреляционного поля будем называть линейным, если каждый его частичный участок является линейным.

На данный момент нет четкого алгоритма построения линейного разбиения корреляционного поля. Можно предложить следующие правила: во-первых, линейность участка определяется визуально, а во-вторых, число таких линейных участков должно быть наименьшим. Поэтому первичной обработкой исходных данных является визуальный анализ корреляционного поля, результатом которого будет последовательность ${\left(l_j\right)}_{j=\mathrm{0,}\mathrm{...,}m-2}$ номеров правых концов интервалов $(x_{l_j-1};x_{l_j})$, $j=\mathrm{0,}\mathrm{...,}m-2$ для определения узлов $u^{\left(j\right)}\in (x_{l_j-1};x_{l_j})$, $j=\mathrm{0,}\mathrm{...,}m-2$ линейного разбиения R.

Для каждого частичного участка линейного разбиения R на отрезке с концами в узлах построим линейную функцию таким образом, чтобы для любых двух соседних частичных участков данного разбиения в их общем узле значения линейных функций совпадали. В результате на интервале $\left(a;b\right)$ будет построена кусочно-линейная функция $\phi$ с фиксированными узлами (1):

$\phi \left(x\right)=\left\{\begin{array}{l}{k}^{\left(0\right)}x+l^{\left(0\right)},x<u^{\left(0\right)},\\ \mathrm{...}\\ k^{\left(j\right)}x+l^{\left(j\right)},u^{(j-1)}<x<u^{\left(j\right)},j=\mathrm{1,}\mathrm{...,}m-\mathrm{2,}\\ \mathrm{...}\\ k^{(m-1)}x+l^{(m-1)},u^{(m-2)}<x,\end{array}\right.$ (2)

где

$k^{\left(j\right)}{u}^{\left(j\right)}+l^{\left(j\right)}=k^{(j+1)}{u}^{\left(j\right)}+l^{(j+1)}$, $j=\mathrm{0,}\mathrm{...,}m-2$. (3)

Значения кусочно-линейной функции $\phi$ в узлах $y^{\left(j\right)}=\phi \left(u^{\left(j\right)}\right)=k^{\left(j\right)}{u}^{\left(j\right)}+l^{\left(j\right)}$ называют узловыми значениями.

Использовав элементарные преобразования, с учетом условий (3), легко получить выражения для коэффициентов $l^{\left(0\right)}$, $k^{\left(j\right)}$, $l^{\left(j\right)}$, $j=\mathrm{1,}\mathrm{...,}m-2$, $l^{(m-1)}$ через узлы и узловые значения. В таком случае кусочно-линейная функция $\phi$ примет вид

$\phi \left(x\right)=\left\{\begin{array}{l}{y}^{\left(0\right)}+k^{\left(0\right)}\left(x-u^{\left(0\right)}\right),x\le u^{\left(0\right)},\\ \mathrm{...}\\ \frac{y^{\left(i\right)}}{u^{\left(i\right)}-u^{(i-1)}}\left(x-u^{(i-1)}\right)-\frac{y^{(i-1)}}{u^{\left(i\right)}-u^{(i-1)}}\left(x-u^{\left(i\right)}\right),u^{(i-1)}\le x\le u^{\left(i\right)},i=\mathrm{1,}\mathrm{...,}m-2\\ \mathrm{...}\\ y^{(m-2)}+k^{(m-1)}\left(x-u^{(m-2)}\right),u^{(m-2)}\le x.\end{array}\right.$ (4)

Геометрически (рис. 1) левый кусок изображается лучом с началом в точке $\left(u^{\left(0\right)};y^{\left(0\right)}\right)$, правый кусок — лучом с началом в точке $\left(u^{(m-2)};y^{(m-2)}\right)$, промежуточный кусок i — отрезком с концами в точках $\left(u^{(i-1)};y^{(i-1)}\right)$ и $\left(u^{\left(i\right)};y^{\left(i\right)}\right)$ для каждого $i=\mathrm{1,}\mathrm{...,}m-2$. Эти отрезки являются звеньями ломаной линии.

Рис. 1. Кусочно-линейная функция, определенная на линейном разбиении корреляционного поля

Из (4) видно, что кусочно-линейная функция $\phi$ зависит от m + 1 параметров. Из них m — 1 узловых значений $y^{\left(0\right)}$, …, $y^{(m-2)}$, перемещающих вершины ломаной вдоль узловых прямых $x=u^{\left(0\right)}$, …, $x=u^{(m-2)}$; $k^{\left(0\right)}$ — угловой коэффициент левого луча, поворачивающий его вокруг начала в точке $\left(u^{\left(0\right)},y^{\left(0\right)}\right)$; $k^{(m-1)}$ — угловой коэффициент правого луча, позволяющий поворачивать его вокруг начала $\left(u^{(m-2)},y^{(m-2)}\right)$.

Таким образом, мы построили класс кусочно-линейных функций с фиксированными узлами. Именно в этом классе, методом наименьших квадратов, мы будем искать аппроксимирующую функцию.

2. МЕТОД НАИМЕНЬШИХ КВАДРАТОВ В КЛАССЕ КУСОЧНО-ЛИНЕЙНЫХ НЕПРЕРЫВНЫХ ФУНКЦИЙ С ФИКСИРОВАННЫМИ УЗЛАМИ

Сумма квадратов отклонений точек корреляционного поля $K$ от кусочно-линейной непрерывной функции (4) имеет вид:

$S=\sum _{i=0}^{n-1}{\left(y_i-\phi \left(x_i\right)\right)}^2=\sum _{j=0}^{m-1}\sum _{i=s_j}^{s_{j+1}-1}{\left(y_i-\phi \left(x_i\right)\right)}^2=\sum _{i=s_0}^{s_1-1}{\left(y_i-\phi \left(x_i\right)\right)}^2+\sum _{j=1}^{m-2}\sum _{i=s_j}^{s_{j+1}-1}{\left(y_i-\phi \left(x_i\right)\right)}^2+\sum _{i=s_{m-1}}^{s_m-1}{\left(y_i-\phi \left(x_i\right)\right)}^2=$

$=\sum _{i=s_0}^{s_1-1}{\left(y_i-y^{\left(0\right)}-k^{\left(0\right)}\left(x_i-u^{\left(0\right)}\right)\right)}^2+\sum _{j=1}^{m-2}\sum _{i=s_j}^{s_{j+1}-1}{\left(y_i-\frac{y^{\left(j\right)}}{u^{\left(j\right)}-u^{(j-1)}}\left(x_i-u^{(j-1)}\right)+\frac{y^{(j-1)}}{u^{\left(j\right)}-u^{(j-1)}}\left(x_i-u^{\left(j\right)}\right)\right)}^2+$

$+\sum _{i=s_{m-1}}^{s_m-1}{\left(y_i-y^{(m-2)}-k^{(m-1)}\left(x_i-u^{(m-2)}\right)\right)}^2$.

Обозначим $x_i^{\left(j\right)}=x_i-u^{\left(j\right)}$, $j=\mathrm{0,}\mathrm{...,}m-2$, тогда

$\begin{array}{c}S=\sum _{i=s_0}^{s_1-1}{\left(y_i-y^{\left(0\right)}-k^{\left(0\right)}{x}_i^{\left(0\right)}\right)}^2+\sum _{j=1}^{m-2}\sum _{i=s_j}^{s_{j+1}-1}{\left(y_i-\frac{y^{\left(j\right)}}{u^{\left(j\right)}-u^{(j-1)}}{x}_i^{(j-1)}+\frac{y^{(j-1)}}{u^{\left(j\right)}-u^{(j-1)}}{x}_i^{\left(j\right)}\right)}^2+\\ +\sum _{i=s_{m-1}}^{s_m-1}{\left(y_i-y^{(m-2)}-k^{(m-1)}{x}_i^{(m-2)}\right)}^2.\end{array}$ (5)

Согласно методу наименьших квадратов оценки параметров функции регрессии, представленной кусочно-линейной непрерывной функцией с фиксированными узлами, находятся как решение оптимизационной модели $S=S\left(k^{\left(0\right)},k^{(m-1)},y^{\left(0\right)},y^{\left(1\right)}\mathrm{,...,}{y}^{(m-2)}\right)\to \min$.

Необходимые условия локального экстремума функции от нескольких переменных имеют вид:

$\frac{\partial S}{\partial k^{\left(0\right)}}=0$, $\frac{\partial S}{\partial y^{\left(0\right)}}=0$, …, $\frac{\partial S}{\partial k^{(m-1)}}=0$. (6)

Введем в рассмотрение величины групповых сумм:

$S_p\left({x^{\left(p\right)}}^2\right)=\sum _{i=s_p}^{s_{p+1}-1}{x_i^{\left(p\right)}}^2$, $S_p\left(x^{\left(p\right)}y\right)=\sum _{i=s_p}^{s_{p+1}-1}{x}_i^{\left(p\right)}{y}_i$, $S_{p+1}\left(x^{\left(p\right)}y\right)=\sum _{i=s_{p+1}}^{s_{p+2}-1}{x}_i^{\left(p\right)}{y}_i$, $S_{p+1}\left({x^{\left(p\right)}}^2\right)=\sum _{i=s_{p+1}}^{s_{p+2}-1}{x_i^{\left(p\right)}}^2$,

$S_p\left(x^{(p-1)}{x}^{\left(p\right)}\right)=\sum _{i=s_p}^{s_{p+1}-1}{x}_i^{(p-1)}{x}_i^{\left(p\right)}$, $p=\mathrm{0,}\mathrm{...,}m-2$;

$S_0\left(x^{\left(0\right)}\right)=\sum _{i=s_0}^{s_1-1}{x}_i^{\left(0\right)}$, $S_0\left(y\right)=\sum _{i=s_0}^{s_1-1}{y}_i$, $S_{m-1}\left(x^{(m-2)}\right)=\sum _{i=s_{m-1}}^{s_m-1}{x}_i^{(m-2)}$, $S_{m-1}\left(y\right)=\sum _{i=s_{m-1}}^{s_m-1}{y}_i$.

Тогда частные производные от $S$, с учетом (5), примут следующий вид.

Частная производная по $k^{\left(0\right)}$:

$\frac{\partial S}{\partial k^{\left(0\right)}}=-2\sum _{i=s_0}^{s_1-1}\left(y_i-y^{\left(0\right)}-k^{\left(0\right)}{x}_i^{\left(0\right)}\right)x_i^{\left(0\right)}$.

Отсюда

$\mathrm{0,5}\frac{\partial S}{\partial k^{\left(0\right)}}=-S_0\left(x^{\left(0\right)}y\right)+y^{\left(0\right)}{S}_0\left(x^{\left(0\right)}\right)+k^{\left(0\right)}{S}_0\left({x^{\left(0\right)}}^2\right).$ (7)

Частная производная по $y^{\left(0\right)}$:

$\frac{\partial S}{\partial y^{\left(0\right)}}=-2\sum _{i=s_0}^{s_1-1}\left(y_i-y^{\left(0\right)}-k^{\left(0\right)}{x}_i^{\left(0\right)}\right)+2\sum _{i=s_1}^{s_2-1}\left(y_i-\frac{y^{\left(1\right)}}{u^{\left(1\right)}-u^{\left(0\right)}}{x}_i^{\left(0\right)}+\frac{y^{\left(0\right)}}{u^{\left(1\right)}-u^{\left(0\right)}}{x}_i^{\left(1\right)}\right)\frac{1}{u^{\left(1\right)}-u^{\left(0\right)}}{x}_i^{\left(1\right)}.$

Тогда

$\mathrm{0,5}\frac{\partial S}{\partial y^{\left(0\right)}}=k^{\left(0\right)}{S}_0\left(x^{\left(0\right)}\right)+y^{\left(0\right)}\left(n^{\left(0\right)}+\frac{S_1\left({x^{\left(1\right)}}^2\right)}{{\left(u^{\left(1\right)}-u^{\left(0\right)}\right)}^2}\right)-y^{\left(1\right)}\frac{S_1\left(x^{\left(0\right)}{x}^{\left(1\right)}\right)}{{\left(u^{\left(1\right)}-u^{\left(0\right)}\right)}^2}-S_0\left(y\right)+\frac{S_1\left(x^{\left(1\right)}y\right)}{u^{\left(1\right)}-u^{\left(0\right)}}.$ (8)

Для произвольного $p=\mathrm{1,}\mathrm{...,}m-3$ находим частную производную по $y^{\left(p\right)}$:

$\begin{array}{c}\frac{\partial S}{\partial y^{\left(p\right)}}=-2\sum _{i=s_p}^{s_{p+1}-1}\left(y_i-\frac{y^{\left(p\right)}}{u^{\left(p\right)}-u^{(p-1)}}{x}_i^{(p-1)}+\frac{y^{(p-1)}}{u^{\left(p\right)}-u^{(p-1)}}{x}_i^{\left(p\right)}\right)\frac{1}{u^{\left(p\right)}-u^{(p-1)}}{x}_i^{(p-1)}+\\ +2\sum _{i=s_{p+1}}^{s_{p+2}-1}\left(y_i-\frac{y^{(p+1)}}{u^{(p+1)}-u^{\left(p\right)}}{x}_i^{\left(p\right)}+\frac{y^{\left(p\right)}}{u^{(p+1)}-u^{\left(p\right)}}{x}_i^{(p+1)}\right)\frac{1}{u^{(p+1)}-u^{\left(p\right)}}{x}_i^{(p+1)}.\end{array}$

Тогда

$\begin{array}{c}\mathrm{0,5}\frac{\partial S}{\partial y^{\left(p\right)}}=-y^{(p-1)}\frac{S_p\left(x^{(p-1)}{x}^{\left(p\right)}\right)}{{\left(u^{\left(p\right)}-u^{(p-1)}\right)}^2}+y^{\left(p\right)}\left(\frac{S_p\left({x^{(p-1)}}^2\right)}{{\left(u^{\left(p\right)}-u^{(p-1)}\right)}^2}+\frac{S_{p+1}\left({x^{(p+1)}}^2\right)}{{\left(u^{(p+1)}-u^{\left(p\right)}\right)}^2}\right)-\\ -y^{(p+1)}\frac{S_{p+1}\left(x^{\left(p\right)}{x}^{(p+1)}\right)}{{\left(u^{(p+1)}-u^{\left(p\right)}\right)}^2}+\frac{S_{p+1}\left(x^{(p+1)}y\right)}{u^{(p+1)}-u^{\left(p\right)}}-\frac{S_p\left(x^{(p-1)}y\right)}{u^{\left(p\right)}-u^{(p-1)}}.\end{array}$ (9)

Частная производная по $y^{(m-1)}$:

$\frac{\partial S}{\partial y^{(m-2)}}=-2\sum _{i=s_{m-2}}^{s_{m-1}-1}\left(y_i-\frac{y^{(m-2)}{x}_i^{(m-3)}}{u^{(m-2)}-u^{(m-3)}}+\frac{y^{(m-3)}{x}_i^{(m-2)}}{u^{(m-2)}-u^{(m-3)}}\right)\frac{x_i^{(m-3)}}{u^{(m-2)}-u^{(m-3)}}-2\sum _{i=s_{m-1}}^{s_m-1}\left(y_i-y^{(m-2)}-k^{(m-1)}{x}_i^{(m-2)}\right).$

Тогда

$\begin{array}{c}\mathrm{0,5}\frac{\partial S}{\partial y^{(m-2)}}=-y^{(m-3)}\frac{S_{m-2}\left(x^{(m-3)}{x}^{(m-2)}\right)}{{\left(u^{(m-2)}-u^{(m-3)}\right)}^2}+y^{(m-2)}\left(\frac{S_{m-2}\left({x^{(m-3)}}^2\right)}{{\left(u^{(m-2)}-u^{(m-3)}\right)}^2}+n^{(m-1)}\right)+\\ +k^{(m-1)}{S}_{m-1}\left(x^{(m-2)}\right)-\frac{S_{m-2}\left(x^{(m-3)}y\right)}{u^{(m-2)}-u^{(m-3)}}-S_{m-1}\left(y\right).\end{array}$ (10)

Частная производная по $k^{(m-1)}$:

$\frac{\partial S}{\partial k^{(m-1)}}=-2\sum _{i=s_{m-1}}^{s_m-1}\left(x_i^{(m-2)}{y}_i-y^{(m-2)}{x}_i^{(m-2)}-k^{(m-1)}{x_i^{(m-2)}}^2\right).$

Тогда

$\mathrm{0,5}\frac{\partial S}{\partial k^{(m-1)}}=y^{(m-2)}{S}_{m-1}\left(x^{(m-2)}\right)+k^{(m-1)}{S}_{m-1}\left({x^{(m-2)}}^2\right)-S_{m-1}\left(x^{(m-2)}y\right).$ (11)

Подставляя правые части равенств (7)–(11) в систему (6), получим систему уравнений:

$\begin{array}{l}{k}^{\left(0\right)}{S}_0\left({x^{\left(0\right)}}^2\right)+y^{\left(0\right)}{S}_0\left(x^{\left(0\right)}\right)=S_0\left(x^{\left(0\right)}y\right),\\ k^{\left(0\right)}{S}_0\left(x^{\left(0\right)}\right)+y^{\left(0\right)}\left(n^{\left(0\right)}+\frac{S_1\left({x^{\left(1\right)}}^2\right)}{{\left(u^{\left(1\right)}-u^{\left(0\right)}\right)}^2}\right)-y^{\left(1\right)}\frac{S_1\left(x^{\left(0\right)}{x}^{\left(1\right)}\right)}{{\left(u^{\left(1\right)}-u^{\left(0\right)}\right)}^2}=S_0\left(y\right)-\frac{S_1\left(x^{\left(1\right)}y\right)}{u^{\left(1\right)}-u^{\left(0\right)}},p=\mathrm{1,}\mathrm{...,}m-3;\end{array}$

$\begin{array}{l}-y^{(p-1)}\frac{S_p\left(x^{(p-1)}{x}^{\left(p\right)}\right)}{{\left(u^{\left(p\right)}-u^{(p-1)}\right)}^2}+y^{\left(p\right)}\left(\frac{S_p\left({x^{(p-1)}}^2\right)}{{\left(u^{\left(p\right)}-u^{(p-1)}\right)}^2}+\frac{S_{p+1}\left({x^{(p+1)}}^2\right)}{{\left(u^{(p+1)}-u^{\left(p\right)}\right)}^2}\right)-y^{(p+1)}\frac{S_{p+1}\left(x^{\left(p\right)}{x}^{(p+1)}\right)}{{\left(u^{(p+1)}-u^{\left(p\right)}\right)}^2}=\\ =\frac{S_p\left(x^{(p-1)}y\right)}{u^{\left(p\right)}-u^{(p-1)}}-\frac{S_{p+1}\left(x^{(p+1)}y\right)}{u^{(p+1)}-u^{\left(p\right)}},\end{array}$ (12)

$\begin{array}{l}-y^{(m-3)}\frac{S_{m-2}\left(x^{(m-3)}{x}^{(m-2)}\right)}{{\left(u^{(m-2)}-u^{(m-3)}\right)}^2}+y^{(m-2)}\left(\frac{S_{m-2}\left({x^{(m-3)}}^2\right)}{{\left(u^{(m-2)}-u^{(m-3)}\right)}^2}+n^{(m-1)}\right)+k^{(m-1)}{S}_{m-1}\left(x^{(m-2)}\right)=\end{array}$

$\begin{array}{l}=\frac{S_{m-2}\left(x^{(m-3)}y\right)}{u^{(m-2)}-u^{(m-3)}}+S_{m-1}\left(y\right),\\ y^{(m-2)}{S}_{m-1}\left(x^{(m-2)}\right)+k^{(m-1)}{S}_{m-1}\left({x^{(m-2)}}^2\right)=S_{m-1}\left(x^{(m-2)}y\right).\end{array}$

Из первого уравнения системы (12) находим

$k^{\left(0\right)}=\frac{S_0\left(x^{\left(0\right)}y\right)}{S_0\left({x^{\left(0\right)}}^2\right)}-y^{\left(0\right)}\frac{S_0\left(x^{\left(0\right)}\right)}{S_0\left({x^{\left(0\right)}}^2\right)}$. (13)

Правую часть равенства (13) подставим во второе уравнение системы (12):

$y^{\left(0\right)}\left(n^{\left(0\right)}+\frac{S_1\left({x^{\left(1\right)}}^2\right)}{{\left(u^{\left(1\right)}-u^{\left(0\right)}\right)}^2}-\frac{S_0{\left(x^{\left(0\right)}\right)}^2}{S_0\left({x^{\left(0\right)}}^2\right)}\right)-y^{\left(1\right)}\frac{S_1\left(x^{\left(0\right)}{x}^{\left(1\right)}\right)}{{\left(u^{\left(1\right)}-u^{\left(0\right)}\right)}^2}=S_0\left(y\right)-\frac{S_1\left(x^{\left(1\right)}y\right)}{u^{\left(1\right)}-u^{\left(0\right)}}-\frac{S_0\left(x^{\left(0\right)}y\right)S_0\left(x^{\left(0\right)}\right)}{S_0\left({x^{\left(0\right)}}^2\right)}$.

Фактически мы из второго уравнения системы (12) исключили $k^{\left(0\right)}$. Подобным образом из последнего равенства системы (12) будем иметь

$k^{(m-1)}=\frac{S_{m-1}\left(x^{(m-2)}y\right)}{S_{m-1}\left({x^{(m-2)}}^2\right)}-y^{(m-2)}\frac{S_{m-1}\left(x^{(m-2)}\right)}{S_{m-1}\left({x^{(m-2)}}^2\right)}.$ (14)

С помощью равенства (14) из предпоследнего уравнения системы (12) исключаем $k^{(m-1)}$ и приходим к уравнению:

$\begin{array}{c}-y^{(m-3)}\frac{S_{m-2}\left(x^{(m-3)}{x}^{(m-2)}\right)}{{\left(u^{(m-2)}-u^{(m-3)}\right)}^2}+y^{(m-2)}\left(\frac{S_{m-2}\left({x^{(m-3)}}^2\right)}{{\left(u^{(m-2)}-u^{(m-3)}\right)}^2}+n^{(m-1)}-\frac{S_{m-1}{\left(x^{(m-2)}\right)}^2}{S_{m-1}\left({x^{(m-2)}}^2\right)}\right)=\\ =\frac{S_{m-2}\left(x^{(m-3)}y\right)}{u^{(m-2)}-u^{(m-3)}}+S_{m-1}\left(y\right)-\frac{S_{m-1}\left(x^{(m-2)}y\right)S_{m-1}\left(x^{(m-2)}\right)}{S_{m-1}\left({x^{(m-2)}}^2\right)}.\end{array}$

К двум уравнениям, полученным из второго и предпоследнего уравнения системы (12) исключением $k^{\left(0\right)}$ и $k^{(m-1)}$, добавляем неиспользованные уравнения системы (12). Приходим к системе уравнений:

$\left\{\begin{array}{l}{a}_{\mathrm{0,0}}{y}^{\left(0\right)}+a_{\mathrm{0,1}}{y}^{\left(1\right)}=b_0,\\ a_{p,p-1}{y}^{(p-1)}+a_{p,p}{y}^{\left(p\right)}+a_{p,p+1}{y}^{(p+1)}=b_p,p=\mathrm{1,}\mathrm{...,}m-\mathrm{3,}\\ a_{m-\mathrm{2,}m-3}{y}^{(m-3)}+a_{m-\mathrm{2,}m-2}{y}^{(m-2)}=b_{m-2},\end{array}\right.$ (15)

коэффициенты которой определяются равенствами:

$a_{\mathrm{0,0}}=n^{\left(0\right)}+\frac{S_1\left({x^{\left(1\right)}}^2\right)}{{\left(u^{\left(1\right)}-u^{\left(0\right)}\right)}^2}-\frac{S_0{\left(x^{\left(0\right)}\right)}^2}{S_0\left({x^{\left(0\right)}}^2\right)},a_{\mathrm{0,1}}=-\frac{S_1\left(x^{\left(0\right)}{x}^{\left(1\right)}\right)}{{\left(u^{\left(1\right)}-u^{\left(0\right)}\right)}^2},b_0=S_0\left(y\right)-\frac{S_1\left(x^{\left(1\right)}y\right)}{u^{\left(1\right)}-u^{\left(0\right)}}-\frac{S_0\left(x^{\left(0\right)}y\right)S_0\left(x^{\left(0\right)}\right)}{S_0\left({x^{\left(0\right)}}^2\right)};$ (16)

$\begin{array}{c}{a}_{p,p-1}=-\frac{S_p\left(x^{(p-1)}{x}^{\left(p\right)}\right)}{{\left(u^{\left(p\right)}-u^{(p-1)}\right)}^2},a_{p,p}=\frac{S_p\left({x^{(p-1)}}^2\right)}{{\left(u^{\left(p\right)}-u^{(p-1)}\right)}^2}+\frac{S_{p+1}\left({x^{(p+1)}}^2\right)}{{\left(u^{(p+1)}-u^{\left(p\right)}\right)}^2},a_{p,p+1}=-\frac{S_{p+1}\left(x^{\left(p\right)}{x}^{(p+1)}\right)}{{\left(u^{(p+1)}-u^{\left(p\right)}\right)}^2},\\ b_p=\frac{S_p\left(x^{(p-1)}y\right)}{u^{\left(p\right)}-u^{(p-1)}}-\frac{S_{p+1}\left(x^{(p+1)}y\right)}{u^{(p+1)}-u^{\left(p\right)}},p=\mathrm{1,}\mathrm{...,}m-3;\end{array}$ (17)

$\begin{array}{c}{a}_{m-\mathrm{2,}m-3}=-\frac{S_{m-2}\left(x^{(m-3)}{x}^{(m-2)}\right)}{{\left(u^{(m-2)}-u^{(m-3)}\right)}^2},a_{m-\mathrm{2,}m-2}=\frac{S_{m-2}\left({x^{(m-3)}}^2\right)}{{\left(u^{(m-2)}-u^{(m-3)}\right)}^2}+n^{(m-1)}-\frac{S_{m-1}{\left(x^{(m-2)}\right)}^2}{S_{m-1}\left({x^{(m-2)}}^2\right)},\\ b_{m-2}=\frac{S_{m-2}\left(x^{(m-3)}y\right)}{u^{(m-2)}-u^{(m-3)}}+S_{m-1}\left(y\right)-\frac{S_{m-1}\left(x^{(m-2)}y\right)S_{m-1}\left(x^{(m-2)}\right)}{S_{m-1}\left({x^{(m-2)}}^2\right)},p=m-2.\end{array}$ (18)

Именно точечные оценки параметров $y^{\left(0\right)}$, …, $y^{(m-2)}$ находятся как решения системы уравнений (15).

Система уравнений (15) называется трехдиагональной. Известен метод прогонки решения такой системы (см., например, (Ильин, Кузнецов, 1985, с. 50)). Для полноты изложения приведем его в разд. 3.

3. МЕТОД ПРОГОНКИ РЕШЕНИЯ ТРЕХДИАГОНАЛЬНОЙ СИСТЕМЫ УРАВНЕНИЙ

Главные миноры матрицы системы (15) обладают свойствами:

$M_0=a_{\mathrm{0,0}}$, $M_1=\left|\begin{array}{cc}{a}_{\mathrm{0,0}}& a_{\mathrm{0,1}}\\ a_{\mathrm{1,0}}& a_{\mathrm{1,1}}\end{array}\right|=a_{\mathrm{0,0}}{a}_{\mathrm{1,1}}-a_{\mathrm{1,0}}{a}_{\mathrm{0,1}}=a_{\mathrm{1,1}}{M}_0-a_{\mathrm{1,0}}{a}_{\mathrm{0,1}}$. (19)

Раскладывая по последней строке, будем иметь

$M_p=\left|\begin{array}{ccccccc}{a}_{\mathrm{0,0}}& a_{\mathrm{0,1}}& 0& \cdots & 0& 0& 0\\ a_{\mathrm{1,0}}& a_{\mathrm{1,1}}& a_{\mathrm{1,2}}& \cdots & 0& 0& 0\\ 0& a_{\mathrm{2,1}}& a_{\mathrm{2,2}}& \cdots & 0& 0& 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& 0& 0& \cdots & a_{p-\mathrm{2,}p-2}& a_{p-\mathrm{2,}p-1}& 0\\ 0& 0& 0& \cdots & a_{p-\mathrm{1,}p-2}& a_{p-\mathrm{1,}p-1}& a_{p-\mathrm{1,}p}\\ 0& 0& 0& \cdots & 0& a_{p,p-1}& a_{p,p}\end{array}\right|=a_{p,p}{M}_{p-1}-a_{p,p-1}{a}_{p-\mathrm{1,}p}{M}_{p-2}.$

Вспомогательные определители получаются из главных миноров заменой последнего столбца столбцом свободных членов:

$M_{0;b}=b_0$, $M_{1;b}=\left|\begin{array}{cc}{a}_{\mathrm{0,0}}& b_0\\ a_{\mathrm{1,0}}& b_1\end{array}\right|=a_{\mathrm{0,0}}{b}_1-a_{\mathrm{1,0}}{b}_0=b_1{M}_0-a_{\mathrm{1,0}}{M}_{0;b}$. (20)

Раскладываем по последней строке

$M_{p;b}=\left|\begin{array}{ccccccc}{a}_{\mathrm{0,0}}& a_{\mathrm{0,1}}& 0& \cdots & 0& 0& b_0\\ a_{\mathrm{1,0}}& a_{\mathrm{1,1}}& a_{\mathrm{1,2}}& \cdots & 0& 0& b_1\\ 0& a_{\mathrm{2,1}}& a_{\mathrm{2,2}}& \cdots & 0& 0& b_2\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& 0& 0& \cdots & a_{p-\mathrm{2,}p-2}& a_{p-\mathrm{2,}p-1}& b_{p-2}\\ 0& 0& 0& \cdots & a_{p-\mathrm{1,}p-2}& a_{p-\mathrm{1,}p-1}& b_{p-1}\\ 0& 0& 0& \cdots & 0& a_{p,p-1}& b_p\end{array}\right|=b_p{M}_{p-1}-a_{p,p-1}{M}_{p-1;b}$.

Из первого уравнения системы (15) получим

$y^{\left(0\right)}=\frac{b_0}{a_{\mathrm{0,0}}}-\frac{a_{\mathrm{0,1}}}{a_{\mathrm{0,0}}}{y}^{\left(1\right)}$. (21)

Методом математической индукции доказывается, что выполняются равенства

$y^{\left(p\right)}=\frac{M_{p;b}}{M_p}-\frac{M_{p-1}{a}_{p,p+1}}{M_p}{y}^{(p+1)},p=\mathrm{1,}\mathrm{...,}m-3;$ (22)

$y^{(m-3)}=\frac{M_{m-3;b}}{M_{m-3}}-\frac{M_{m-4}{a}_{m-\mathrm{3,}m-2}}{M_{m-3}}{y}^{(m-2)},p=m-3.$ (23)

Уравнение (23) и последнее уравнение системы (15) приводят к системе двух линейных уравнений с переменными $y^{(m-3)}$ и $y^{(m-2)}$.

Методом Крамера решения систем линейных уравнений приходим к равенству

$y^{(m-2)}=M_{m-2;b}/M_{m-2}$. (24)

Подставляя правую часть (24) в (23), находим . Применяя к равенству (22) индукцию спуска, получим все значения $y^{\left(p\right)}$ от $p=m-4$ до $p=1$ включительно. Подставляя $y^{\left(1\right)}$ в (21), находим $y^{\left(0\right)}$, что завершает решение системы (15).

Для решения системы (12) необходимо выполнить еще два шага. Значение $y^{(m-2)}$ подставить в равенство (14) и определить $k^{(m-1)}$, а $y^{\left(0\right)}$ подставить в равенство (13) и найти $k^{\left(0\right)}$. Решив систему (12), мы найдем точечные оценки параметров аппроксимирующей функции (4).

Равенства (21)–(24) совместно с (13) и (14) составляют каркас алгоритма определения точечных оценок параметров аппроксимирующей функции. Остановимся на этом подробнее.

4. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ АППРОКСИМИРУЮЩЕЙ ФУНКЦИИ В КЛАССЕ КУСОЧНО-ЛИНЕЙНЫХ ФУНКЦИЙ С ФИКСИРОВАННЫМИ УЗЛАМИ

Входными данными являются массивы значений ${\left(x_i\right)}_{i=\mathrm{0,}\mathrm{...,}n-1}$, ${\left(y_i\right)}_{i=\mathrm{0,}\mathrm{...,}n-1}$ и ${\left(u^{\left(j\right)}\right)}_{j=\mathrm{0,}\mathrm{...,}m-2}$. По ним вычисляются групповые суммы, которые сохраняются в специальных переменных и массивах. Затем находим коэффициенты и свободные члены системы (15) и сохраняем их значения в отдельном массиве. После этого рассчитываем значения главных миноров и вспомогательных определителей и сохраняем их в двух массивах. Применяя равенства (21)–(24) к элементам двух последних массивов, мы найдем $y^{\left(p\right)}$ для $p=\mathrm{0,}\mathrm{...,}m-2$. В заключение равенствами (13) и (14) получаем $k^{\left(0\right)}$, $k^{(m-1)}$. Алгоритм возвращает массив значений параметров $y^{\left(p\right)}$, $p=\mathrm{0,}\mathrm{...,}m-2$, $k^{\left(0\right)}$, $k^{(m-1)}$.

Данный алгоритм реализован на языке Python в виде функции Gluing_m_Linear_Pieces().

Вычисление значения аппроксимирующей функции в точке реализуется с помощью функции mLinearGluing(), а на векторе — с помощью функции mLinearGluingV().

Напомним, что последовательность узлов ${\left(u^{\left(j\right)}\right)}_{j=\mathrm{0,}\mathrm{...,}m-2}$ построена в разд. 1 по номерам $l_j$, $j=\mathrm{0,}\mathrm{...,}m-2$ правых концов интервалов $(x_{l_j-1};x_{l_j})$. Функция IndexKnot() принимает массивы ${\left(l_j\right)}_{j=\mathrm{0,}\mathrm{...,}m-2}$ и ${\left(x_i\right)}_{i=\mathrm{0,}\mathrm{...,}n-1}$ и возвращает массив узлов ${\left(u^{\left(j\right)}\right)}_{j=\mathrm{0,}\mathrm{...,}m-2}$ как середины $u^{\left(j\right)}=\mathrm{0,5}(x_{l_j-1}+x_{l_j})$, $j=\mathrm{0,}\mathrm{...,}m-2$ соответствующих интервалов.

Также организована функция mLinearGluingDeter() для вычисления коэффициента детерминации $R^2$.

5. ПРИМЕР

Работу алгоритма из разд. 4 проиллюстрируем на следующем примере. В табл. 1 приведены значения ${\left(x_i\right)}_{i=\mathrm{0,}\mathrm{...,}45}$, ${\left(y_i\right)}_{i=\mathrm{0,}\mathrm{...,}45}$ соответствующих значений объясняющей $X$ и зависимой $Y$ переменных. Согласно алгоритму, составляем корреляционное поле (рис. 2).

Таблица 1. Значения соответствующих значений объясняющей X и зависимой Y переменных