МЕТОД ПОСТРОЕНИЯ НЕЭЛЕМЕНТАРНЫХ ЛИНЕЙНЫХ РЕГРЕССИЙ НА ОСНОВЕ АППАРАТА МАТЕМАТИЧЕСКОГО ПРОГРАММИРОВАНИЯ

М атематические проблемы управления

УДК 519.862.6 DOI: http://doi.org/10.25728/pu.2022.4.1

МЕТОД ПОСТРОЕНИЯ НЕЭЛЕМЕНТАРНЫХ ЛИНЕЙНЫХ РЕГРЕССИЙ НА ОСНОВЕ АППАРАТА МАТЕМАТИЧЕСКОГО ПРОГРАММИРОВАНИЯ

М.П. Базилевский

Аннотация. Рассматривается проблема построения неэлементарных линейных регрессий, состоящих из объясняющих переменных и всевозможных комбинаций их пар, преобразованных с помощью бинарных операций минимум и максимум. Задача построения таких моделей формализована в виде задачи частично-булевого линейного программирования. Регулируя в ней ограничения на бинарные переменные, можно контролировать структурную спецификацию неэлементарной линейной регрессии, а именно количество входящих в нее регрессоров, их типы и состав объясняющих переменных. При этом оценки параметров модели находятся приближенно с помощью метода наименьших квадратов. К достоинствам сформулированной задачи относится то, что число ее ограничений не зависит от объема выборки, а знаки оценок при объясняющих переменных согласуются со знаками коэффициентов их корреляции с зависимой переменной. Показано, как на начальном этапе отсекать регрессоры, чтобы сократить время решения задачи и сделать модель вполне интерпретируемой. Построена неэлементарная линейная регрессия для моделирования железнодорожных грузоперевозок в Иркутской области и дана ее интерпретация.

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

ВВЕДЕНИЕ

При проведении регрессионного анализа [1, 2] на основе экономических данных особое внимание уделяется построению производственных функций (ПФ), представляющих собой математические зависимости между объемами выпуска продукции и факторами производства. Теории, методам и применению ПФ целиком посвящена монография [3], выпущенная еще в 1986 г. В ней рассмотрены следующие ПФ: линейная, многорежимная, Кобба -Дугласа, Леонтьева, Аллена, CES (Constant Elasticity of Substitution - с постоянной эластичностью замены факторов), LES (Linear Elasticity of Substitution - с линейной эластичностью замены факторов), Солоу. В настоящее время появляются и исследуются новые модификации ПФ, которые активно применяются в эконометрических исследованиях [4-6]. В данной статье речь пойдет о построении неэлементарных регрессионных моделей, специфицированных на основе известной ПФ Леонтьева:

У = т1п{а1х11, а2х12,..., ах} + £,., г = 1, п,

где п - объем выборки; I - количество объясняющих переменных; у., г = 1, п, - значения объясняемой переменной у ; х^ , г = 1, п, } = 1,1, - значения объясняющих переменных х1, х2, ..., х1; а ^

, ] = 1,1, - неизвестные параметры; вг, г = 1, п, -ошибки аппроксимации. С позиции теории ПФ переменная у в уравнении (1) трактуется как объем выпуска продукции, а х1, ..., хп - как показатели факторов производства.

Отметим, что в монографии [3] выделена еще и «параллельная» функция Леонтьева

у,. = 1шп{апх;1, а12х;2,...,а1;х,7} + ...

+тт{а,1 х,, ак2х2,..., ад,} + В, г = 1, п,

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

порциями факторов, использующих общие ресурсы. Для двух факторов производства х1 и х2 «параллельная» функция Леонтьева называется функцией линейного программирования.

В монографии [7] отмечается, что для нахождения оценок параметров ПФ Леонтьева (1) можно применять методы негладкой оптимизации [8-10], которые, как правило, являются труднореализуемыми. Поэтому в работе [7] задача точного оценивания ПФ (1) с помощью метода наименьших модулей (МНМ) сведена к задаче частично-булевого линейного программирования (ЧБЛП). Вместе с тем в монографии [7] предложен способ приближенного оценивания ПФ Леонтьева, основанный на переборе оценок из предварительно сформированной области определения.

В статье [11] предложена функция, противоположная по смыслу ПФ (1):

индексов переменных; а , 1 = 0,1 + Сг2, 1.,

у.1 = тах |а1 хп, а2х.12,

а1хи } + 8,-

I = 1, п,

а в статье [12] рассмотрен симбиоз функций (1) и (2):

У1 = тп а2xi2,, аХ } + +тах {Р1Х1, Р2 X 2 ,, Рл } + 8.', i =1, п

Задачи точного оценивания параметров регрессий (2) и (3) с помощью МНМ сведены в статьях [11] и [12] к соответствующим задачам ЧБЛП. Повышенное внимание к построению регрессионных моделей с использованием аппарата математического программирования в современной научной литературе (см., например, работы [13-15]) объясняется тем, что за последние годы была существенно развита технология решения задач ЧБЛП.

Данная статья посвящена оцениванию специфицированных на основе ПФ Леонтьева регрессионных моделей с помощью метода наименьших квадратов (МНК) [1, 2]. Впервые такая задача была сформулирована автором в работе [16] для регрессии (1) с двумя объясняющими переменными. А в статье [17] была предложена неэлементарная линейная регрессия (НЛР) вида

Уi =ао+Е

а зха +

1=1

+

"-1 _ Еа 1+1 ттКЦд, Х АЦ,2} +8г- , i = 1

1=1

где ц ;1 и ц 12, 1 = 1, С1 - элементы первого и второго столбца индексной матрицы М^ , содержащей по строкам всевозможные комбинации пар

1 = 1, С1 , - неизвестные параметры. Считается, что значения всех переменных в уравнении (4) строго положительны.

Как видно, НЛР относится к классу нелинейных по параметрам моделей. Но если придать всем

параметрам 11, 1 = 1, С2, определенные значения,

то регрессия становится линейной и нахождение

оценок ее параметров а 1, 1 = 0,1 + С2, с помощью

МНК не вызывает трудностей. В статье [17] установлено, что оптимальные с точки зрения МНК

оценки параметров 1 , 1 = 1, С2, для НЛР принад-

лежат промежуткам

фтп, ), 1=й=

где

= тах

Я (11 = тт

'Ц11

1 =

^, ^,, ^ |. Точки х 1 =ят!п и

[ х1ц 12 х2,Ц 12 хП,Ц 12 ]

X1 = X т^ нельзя использовать из-за возникновения

совершенной мультиколлинеарности переменных.

Благодаря этим свойствам, в работе [17] был предложен способ приближенного МНК-оценивания НЛР (4), основанный на переборе значений параметров 1 ^, 1 = 1, С2, из промежутков

(5)-

К сожалению, при построении НЛР (4) с ростом числа объясняющих переменных I существенно возрастает общее количество ее регрессо-ров. Поэтому появляется необходимость в решении задачи отбора некоторого числа наиболее «информативных» регрессоров (ОИР) [7]. Специально для этого в работе [18] были разработаны две стратегии. Каждая из них предполагает формирование по указанному алгоритму множества альтернативных вариантов регрессий, для каждой из которых реализуется описанный в статье [17] способ приближенного МНК-оценивания, а затем выбирается модель с наименьшей величиной суммы квадратов остатков. Главным недостатком предложенного в работе [18] подхода к построению НЛР является то, что он основан на методе полного перебора всех возможных альтернатив, поэтому на решение задачи ОИР может уходить слишком много времени. Более перспективным выглядит следующий подход с использованием аппарата ЧБЛП.

п

В статье [19] задача ОИР при оценивании линейной регрессии с помощью МНК была сведена к задаче ЧБЛП. При этом нерешенной оставалась проблема выбора большого положительного числа М, влияющего как на скорость, так и на само решение задачи, пока не появилась работа [20]. В ней сформулирована задача ЧБЛП, которая приводит к построению линейной регрессии с заданным числом объясняющих переменных, в которой знаки МНК-оценок согласованы со знаками коэффициентов корреляции между переменными у и х^,

1 = 1,1. В ходе вычислительных экспериментов был подтвержден сделанный в статье [21] вывод, что такая задача с ограничениями на знаки коэффициентов решается на порядок быстрее, чем без них. Главной целью данной работы является сведение задачи построения НЛР к эффективно решаемой задаче ЧБЛП, рассмотренной в статье [20].

1. МЕТОД ПОСТРОЕНИЯ НЕЭЛЕМЕНТАРНЫХ _ЛИНЕЙНЫХ РЕГРЕССИЙ_

Уравнение НЛР (4) содержит только одну бинарную операцию - минимум. Здесь и далее под бинарной операцией минимум (максимум) понимается математическая операция, принимающая два аргумента и возвращающая их минимум (максимум). Дополним эту регрессионную модель ре-грессорами с бинарной операцией максимум: I с2

у, =а0 + £ах + £а1+1 ,, ^А».2} +

1=1

1=1

+ £а,+,+а тах{х,ц(1, ^2 } + В, г = 1 п,

1=1

Общее число регрессоров в уравнении (6) стало гораздо больше, чем в уравнении (4), и составляет I + 2С,.

Уравнение НЛР вида (6) вводится впервые, поэтому ставится задача формализовать процесс построения этой модели в виде задачи ЧБЛП. Это можно сделать следующим образом.

Для каждого параметра X,, 1 = 1, С2, из уравнения (6) определим промежутки значений по формулам (5). Затем равномерно разобьем каждый из этих промежутков р точками и сформируем

матрицу Л = (Х*,к ), 1 = 1, С,2 , к = 1, р , элемент 1*,к которой показывает к -е значение параметра X, для 1 -й пары переменных. Заменив в уравнении (6) неизвестные параметры X, на известные элементы матрицы Л, получим:

I С р . .

£а х+££а-кт1п {х>д, 2}

у. =а0 + £а х +

1=1 1=1 к=1

С1 р , ч _

££а+ктах{х,»^ }+в, г=1 n,

1=1 к =1

где а-к, 1 = 1, С1 , к = 1, р , - неизвестные параметры для регрессоров с бинарной операцией минимум, а а +к, 1 = 1, С,2 , к = 1, р , - неизвестные

параметры для регрессоров с бинарной операцией максимум. В модели (7) общее число регрессоров еще больше, чем в модели (6), и равно I + 2рС,2. Например, если общее число переменных I = 100, а число разбиений р = 10, то регрессия (7) будет иметь 99 100 регрессоров.

Сделаем в уравнении (7) замену

1 = т1п {*,, Ь}, 4 = тах {х,1, А.>;,„ },

г = 1, п, 1 = 1, С2 , к = 1, р . Получим модель множественной линейной регрессии:

1 г

у. =а0 +£а Л + ££0^ +

1=1 1=1 к=1

Г-2

С1 р _

+в, г=1, п

1=1 к=1

Далее, как это сделано в статье [19], сведем задачу ОИР для линейной регрессии (8), оцениваемой с помощью МНК, к задаче ЧБЛП. Для этого предварительно проведем нормирование всех переменных из уравнения (8) по известному правилу, вычитая из каждого значения переменной ее среднее арифметическое и деля результат на стандартное отклонение.

Составим для модели (8) уравнение стандартизованной регрессии

Щ

1 У

1=1

1=1 к =1

1=1 к=1

где щ - нормированная переменная у ; , } = 1,1, - нормированные переменные х^, 1 = 1,1; И -к, И+,

1 = 1, С,2 , к = 1, р, - нормированные переменные

- +

1к, ¿д,

1-к, , 1 = 1, С2, к = 1, р ; в;, 1 = 1, /, и в-к, Р+к:

1 = 1, С2 , к = 1, р, - неизвестные стандартизованные коэффициенты; £г, г = 1, п, - новые ошибки аппроксимации.

п

Для модели (9) МНК-оценки находятся по формуле

Р ^ XX ' \ 5 Л

Яхх Я х; Я + х;

где яхх = Я ; х Я ;; Я + ;;

Я + V ; х Я + ;; Я + + ;;

- корреляционная

блочная матрица размера (I + 2рСг2)*(/ + 2рС?), составленная из блоков

Кхх = ( 1), 1 = й, к = 1П;

Ях; =(гхг-), * = , 1 = й?, к = 1Гр;

Я + =

хг

я _ =

(, *=11, 1=1,С2, к=1 р; (), 1=1С', к=1, р, *=1,1;

я.-.- =(г г-), *1 = *2 =1, p, 1 = 1, См к =1 р;

; ; у г*1*2 ;к1

Я-.+ =(г. .+ ), *1 = С1 *2 = 1 р, 1 = 1, См к =1 р;

Я + =

; х

\ ;-]кх, ) '

1 = 1, С, , к = 1, р , * = 1,1;

Я + =( г+ _

Я;;+ =( г;+ ;+ ; ; \ к

*1 = С1, *2 =1 р> 1 = 1, С,2, к =1, р;

*1 = С1, *2 =1 P, 1 = 1, С,, к =1, р;

Яж = (Яух Я _ Я^ +) - корреляционный блочный вектор размера (, + 2рС2 )х 1, составленный

из блоков Яух =( ^ ) , 1 = 1, I ; Я ; =( Гуг_к ) ,

1 = ГС2, к = 1Гр; я^+=(Гу;.), 1 = Сс2, к = 1Гр

Коэффициент детерминации модели (9) находится по формуле

Я Гух1 Р; + Е Е^ Р- + ЕЕ^ Р+к

1=1

1=1 к =1

1 =1 к=1

Тогда с использованием формул (10) и (11) сформулируем задачу ОИР для линейной регрессии (8):

Я ^ тах,

1 С р

-(1 1 ) М <Е'1 Рк + Е Е 1 Р-

+

+ Е Е^Ик Р+к - Гух, < С1 -51 )М , 1 = 1, I,

-(1 -5Д)М<ЕР* + ЕЕГ. _ Р*1

;1к 12

+

и

е ъ

; 1к 12

*1 =1 *2 =1

р+ - г,;- < (1 -5к )м,

у;]к

1 = 1С, к = 1, р, 1 с р -(1 -5+к)М <2/.Р* +ЕЕ-, Р-1,

+

г

^Егг!^Р+1*2 -Гу1 < (1 -5+к)М

*1 =1 *2 =1

1 = 1, С,2, к = 1, р,

-5 М <Р . <5М , 1 = 1,1:

-5.М <Р-к <5кМ , 1 = 1, С2, к = 1, р,

-5+кМ < Р+к < 5+кМ , 1 = 1, С, , к = 1, р ,

5. е {0,1}, 1 = 1,1.

5-к е {0,1} , 1 = 1, С2, к = 1, р ,

51 е {0,1} , 1 = 1, С2 , к = 1, р ,

, с} р с} р

Е +Е Е 5-к+Е Е 5+к =

1=1 1=1 к =1 1=1 к =1

где т - заданное число регрессоров; 5;., 1 = 1,!, -булевы переменные, заданные по правилу

Г1, если 1-я переменная входит в регрессию,

5/ = 1

10, в противном случае;

5 к, 1 = 1, С2, к = 1, р, - булевы переменные, заданные по правилу

1, если 1-я бинарная операция минимум с к -м преобразованием входит в регрессию, 0, в противном случае;

5-к =

5 +к, 1 = 1, С2 , к = 1, р, - булевы переменные, заданные по правилу

к=1

*=1 к=1

,=1

,=1

* =1 я, =1

81 =

1, если 1-м бинарная операция

максимум с к-м преобразованием

входит в регрессию,

0, в противном случае;

М - большое положительное число.

Достоинством задачи ЧБЛП (12)-(22) является то, что число ее ограничений не зависит от объема выборки п .

В задаче ЧБЛП (12)-(22) стратегия построения НЛР регулируется ограничениями на бинарные переменные. Рассмотрим следующие стратегии.

Стратегия 1. Отбор т регрессоров в линейной регрессии (7).

Для этого просто нужно решить задачу (12)-(22). В этом случае итоговая модель может содержать несколько регрессоров с одинаковой бинарной операцией и с одинаковой парой переменных, но с разными значениями параметра X .

Стратегия 2. Приближенное оценивание НЛР (6) с помощью МНК (без отбора регрессоров). Для этого нужно решить задачу с целевой

функцией (12), ограничениями (13)-(21) и

р р

£8 -к = 1, £ 8; = 1, 1=1, с

2

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

Стратегия 3. Отбор т регрессоров в НЛР (6).

Для этого нужно решить задачу с целевой функцией (12), ограничениями (13)-(22) и

£8 -к *1, £8*1, 1 = I С2.

Заметим, что регулируя ограничения на бинарные переменные, можно контролировать тип входящих в НЛР (6) регрессоров. Так, например, если добавить в задачу (12)-(22) ограничения

С! р С1 р

££8-к = 0, ££8 ]к = 0, то получим задачу ОИР

1=1 к=1 1=1 к=1

для линейной регрессии. Если добавить ограниче-

I С, р

ния £81 = 0, ££8+к = 0, то получим задачу 1=1 1=1 к=1 ОИР для регрессии только с бинарными операция-

ми

минимум, а если ограничения £8  = 0,

1=1

££8-к = 0 - то задачу ОИР для регрессии только

1=1 к=1

с бинарными операциями максимум.

Помимо этого можно контролировать состав входящих в модель переменных. Для этого введем

бинарную матрицу V = {у у}, г = 1,, + 2 рС, ,

1 = 1,,, в которой

1, если 1-я переменная входит в г -й регрессор модели (7), 0, в противном случае.

уи =

Тогда интеграция в задачу (12)-(22) линейных ограничений

£ уа51 +£ £ ^+р( 1 -1)6-/к +

1=1

1 =1 к=1

+££ у,,+рс2+р( н)5+к *1,г=1, 1=1 к =1

позволяет строить НЛР с т регрессорами, в которую каждая объясняющая переменная входит не более одного раза. В этом случае автоматически выполняются условия (23).

К сожалению, для задачи (12)-(22) не до конца ясно, как задавать большие числа М . Для решения этой проблемы поступим так, как предложено в работе [20]. Заменим ограничения (13)-(18) на следующие:

' С р

-(1 -51)М1 *£ 1 Рк +££1 р«

к =1

+

«=1 к =1

+££ 1 Р+ * а- 5 )м+1 = 1,,

«=1 к=1 1

I

-(1-5 -к) м-и *£ гг р« + £ £ Г р-1^2 +

5=1 «1 =1 «2 =1

г'2 С р

££гг. р;2-г*

+££Г+

«1 =1 «2 =1

у1}к

1 = 1, С, к=1, р,

-(1- 5+к ) ^ р5 + £ £ Г ¿к Р-52 +

5=1 «1 =1 «2 =1

г'2 С р

£ К* - Гуг+к *

у?1к

1 = 1, С, , к = 1, р ,

0 *Р 1 *5 М р 1, 1е Ъ, 5М *Р1 *0, 1еЪ,

0 *Р-к *5кМр1к, 1к е ^

к=1

к=1

Г<2

р

к=1

к=1

5 jkMß-k <ß jk < ^ j k e Jß-,

0 <ßjk <5jkMk, j, k e Jßj,

5iMk <ßjk < 0, j, k e Jßj,

где Jр и Jр - индексные множества, построенные из множества {1,2,..., l}, элементы которых удовлетворяют условиям ryx > 0 и ryx < 0 соответственно; Jp+_ и J__ - индексные множества, построенные из множества {{1, 2},..., {1, p}, {2,1},..., {2, p},..., {Cl 1},..., {Cl p}} , элементы

которых удовлетворяют условиям r _ > 0 и

yzjk

r _ < 0; Jn++ и J~ - индексные множества, по-

yzjk Р Р

строенные из множества {{1,2},..., {1, p}, {2,1},..., {2, p},..., {Cl 1},...,{Cl p}}, элементы

которых удовлетворяют условиям r + > 0 и

yzjk

<m, = 1/rj, ^; jk =1 /г*, Mj =1 /r, jk' j =1' C'' k = p •

Для нахождения чисел M~ в ограничениях

нужно решить серию из l задач линейного программирования с целевыми функциями Uj ^ min при ограничениях

0 <ß j < Mру, j е Jp+,

Mßj < ß j < 0 , j е J_ ,

о <ß j, < мр, j, к e j;,

M- <ß -k < J, k e Jß-' 0<ßjk <, j k e Jßj,

Mß)k <ß)k <j k e J-,

Zj ßk jZZj ßk +

к=1 »=1 к=1 J

Г1

Cl p _

ZZj &k-ry, = uj,j=1,l,

s=1 k=1 J

l C p

Zr - ß»+ZZr- - ß-» +

t—' k ' » ¿—¡¿—i z-,„Z -k 1 »1»2

+ZZr+ - ßj» -r - = j, j=1, с2,k=1, p,

i—ti—t za„z k 1 »1»2 yz-k l

Z^ ß» jZ Z ßv2

»1 =1 »2 =1

ßt1»2-ryzj = utk, j =1, C2,k =1, P,

»1 =1 »2 =1 l

zjk 2 yzjk

с с

Z r, ß j +Z Z ßjk +Z Z <1

yzjk jk yzjk

j=1 j=1 k=1 J j=1 k =1 J

Для нахождения чисел M+ нужно решить се-

uj

рию из l задач линейного программирования с целевыми функциями Uj ^ max при ограничениях

(34)-(43). Аналогично находятся числа M - , M j ,

ujk ujk

M- , M ++ решением серии из pCf задач линей-

ujk ujk

ного программирования с целевыми функциями

ujk ^ min , u-k ^ max, u +jk ^ min , u +jk ^ max

соответственно при ограничениях (34)-(43).

Таким образом, решение задачи ЧБЛП с целевой функцией (12) и ограничениями (19)-(22), (25)-(33) приводит к построению линейной регрессии (7) с m регрессорами, в которой знаки оценок ß-параметров согласованы со знаками соответствующих коэффициентов корреляции ре-грессоров с переменной y , т. е. справедливы неравенства ß jryXj > 0, j =1, l; ß-k r^ > 0, ßjjk r^ > 0,

7 = 1, С2 , к = 1, р . Стратегия построения НЛР в этой задаче по-прежнему регулируется, например, ограничениями (23) и (24) на бинарные переменные.

В работах [20, 21] экспериментально установлено, что задача ЧБЛП (12), (19)-(22), (25)-(33) решается на порядок быстрее, чем задача (12)-(22). Кроме того, из-за согласованности знаков оценок Р-параметров и знаков соответствующих коэффициентов корреляции для полученной регрессии становятся справедливыми формулы для абсолютных вкладов переменных в общую детерминацию

R2:

Cj = Vyxj ßj. j = 1 l. CZj^C = ryzjkßjk,

»1 =1 »2 =1

ч = Р7к, 7=1, сг,к = 1, р,

по которым можно судить о степени влияния каждого регрессора на переменную у .

Сделаем два важных замечания относительно решения задачи (12), (19)-(22), (25)-(33).

Замечание 1. Поскольку в результате решения задачи знаки оценок Р-параметров согласуются со знаками соответствующих коэффициентов корре-

»=1

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

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

Замечание 2. Пусть, например, модель (8) имеет при параметре а- регрессор г- = т1п{х1, 8х2}. Тогда при переходе к кусочно-заданному представлению при параметре а- будет либо переменная х , либо 8х . Если окажется, что г > 0, то

1' 2 у*п

оценка параметра а1 1 гарантированно будет положительной и переменные х1 и 8х2 будут влиять на у со знаком «плюс». В этом случае оба коэффициента корреляции гух1 и гух2 должны быть положительны, иначе возникает проблема с интерпретацией модели. А если окажется, что г _ < 0, то

у^11

оценка параметра а1 1 гарантированно будет отрицательной и переменные х1 и 8х2 будут влиять на у со знаком «минус». В таком случае оба коэффициента корреляции гх и гухдолжны быть отрицательны. Из всего этого следует, что после согласования с экспертами знаков коэффициентов корреляции гух , 1 = 1,,, необходимо сформировать

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

1

^2

или

1 = 1, С, , к = 1, р ,

у уг+к ух»11 ух»- '

или

угк

^11

^12

1 = 1, С, , к = 1, р .

Исключение противоречивых переменных естественным образом уменьшит время решения задачи построения НЛР. Это время можно еще значительно уменьшить, если дополнить выражения (45) и (46) условиями

Г > Г , Г , > Г

уг)к уг)к

> г , 1 = 1, С, , к = 1, р

где г - выбранное из промежутка [0, 1) число. Чем больше число г , тем меньше становится количество переменных и время решения задачи.

2. МОДЕЛИРОВАНИЕ

Для построения НЛР были собраны ежегодные статистические данные за период с 2000 по 2020 г. для зависимой переменной у - отправление грузов железнодорожным транспортом общего пользования в Иркутской области (млн руб.), и шестидесяти двух переменных х1, х2, ..., х62, предположительно влияющих на у . Сначала из этого списка было исключено шесть переменных, у которых значение коэффициента корреляции с у по абсолютной величине не превышало 0,2. Затем значения коэффициентов корреляции для оставшихся пятидесяти шести переменных были переданы двум экспертам, представляющим Управление Восточно-Сибирской железной дороги. Их задачей было исключить те переменные, для которых знаки коэффициентов корреляции с у не соответствуют экономическому смыслу решаемой задачи. В результате работы экспертов осталось восемь факторов:

х2 - процент трудоспособного населения от общей численности;

х3 - численность рабочей силы (тыс. чел.); х5 - численность пенсионеров (тыс. чел.); х8 - число собственных легковых автомобилей на 1000 человек населения (шт.);

х18 - число предприятий и организаций; х20 - кредиторская задолженность организаций (млн руб.);

х22 - производство электроэнергии (млрд кВтч);

х58 - тарифы на грузовые перевозки (железнодорожный транспорт), усл.ед.

Значение переменной х58 за 2001 г. было назначено равным 1000 усл. ед. С его помощью по известным индексам тарифов были найдены оставшиеся значения переменной х58.

Значения коэффициентов корреляции с переменной у для отобранных переменных составляют соответственно гух = 0,785 , г = 0,543,

ух2 ух3

Гух5 = -0,483 :

г = -0,446,

ух8

^ = 0,538 =

Гух20 = ^,204, Гух22 = 0,476, Гух58 = "А465.

Влияние выбранных переменных на переменную y можно обосновать следующим образом:

• рост численнсти рабочей силы x2 и x3, а также числа предприятий x18 и количества производимой электроэнергии x22, приводит к увеличению объемов производимой регионом продукции, что влечет за собой повышение спроса на грузовые перевозки ж/д транспортом, в то время как рост значения переменной x5 тормозит развитие экономики, снижая спрос на перевозки;

• избыток собственных автомобилей x8 у населения снижает спрос как на пассажирские, так и на грузовые перевозки ж/д транспортом;

• рост суммарного объема кредиторской задолженности организаций x20 негативно сказывается на экономике региона, поскольку, например, может повлечь за собой наложение различных штрафных санкций;

• увеличение тарифов на грузовые перевозки x58 естественно снижает спрос на перевозки грузов по железной дороге.

Затем для каждой пары отобранных переменных по формулам (5) были определены промежутки для значений параметров Xj. После этого для

формирования матрицы Л каждый промежуток был равномерно разбит четырьмя точками. В результате удалось сформировать 4С82 = 112 переменных z-k, j = 1, 28, k = 1, 4, преобразованных с помощью бинарной операции минимум, и столько же переменных z-k, j = 1, 28, k = 1, 4, преобразованных с помощью операции максимум. Далее из этих 224 переменных были исключены те, для которых не выполняются условия (45)-(47) при r = 0,2 . Таких переменных оказалось 140. В итоге к построению НЛР мы подошли, имея в наличии 92 переменных, из которых 8 объясняющих и 84 преобразованных с помощью функций min и max.

Построение НЛР осуществлялось на основе решения задачи ЧБЛП с целевой функцией (12) и ограничениями (19)—(21), (25)-(33). Подчеркнем, что ограничение (22) на число входящих в модель регрессоров не ставилось. Для того чтобы в итоговую модель каждая объясняющая переменная входила не более одного раза, были учтены ограничения (24). Для решения задач ЧБЛП использовался решатель LPSolve IDE, а для формирования математических моделей задач для этого решателя была разработана специальная программа в среде программирования Delphi. Сначала с помощью

этой программы были найдены неизвестные числа в ограничениях (25)-(27). Для этого было решено 184 задачи линейного программирования с соответствующими целевыми функциями и линейными ограничениями (34)-(43). Затем с использованием найденных чисел и разработанной программы для решателя LPSolve была сформулирована задача ЧБЛП (12), (19)-(21), (24)-(33), содержащая 284 ограничения, 92 вещественных и 92 бинарных переменных. Решение осуществлялось на персональном компьютере с процессором Intel Core i5 (3.40 ГГц, 4 ядра) и объемом оперативной памяти 8 ГБ. В результате примерно за 30 секунд была построена следующая НЛР:

у = -24,5274 +1,1895 min {х,, 0,00093 3jq Л -

-0,0196min{x5, 0,006754x20}-

-0,0323min{x8, 0,11725x58} +

+ 0,0254max{x3, 23,079x22},

где в скобках под коэффициентами указаны значения ¿-критерия Стьюдента, а над коэффициентами - абсолютные вклады регрессоров в общую детерминацию, найденные по формулам (44). Оказалось, что все регрессоры значимы по критерию Стьюдента для уровня значимости а = 0,05 .

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

Коэффициент детерминации НЛР (48) Я2 = 0,946183, что говорит о высоком качестве построенной модели.

Значения коэффициентов вздутия дисперсии для регрессоров модели (48) не превышают 10, что говорит об отсутствии в ней мультиколлинеарно-сти. Стоит отметить, что контролировать мульти-коллинеарность в сформулированной задаче ЧБЛП пока тоже нельзя.

Таким образом, выполняются все условия, чтобы отнести построенную модель (48) к вполне интерпретируемым.

Модель (48) в кусочно-заданном виде представлена в таблице.

10

CONTROL SCIENCES No.4 • 2022

Уравнения модели (48) для различных диапазонов значений переменных

Уравнение НЛР Диапазоны значений переменных

у = -24,527 + 0,001 1х18 - 0,00013х20 -0,003 8х58 +0,0254х3 х х х х > 0,000933,-^ > 0,00675, > 0,117,^ > 23,08 х18 х20 х58 х22

V = -24,527 + 0,001Ц8 -0,00013х20 -0,0038х58 +0,5857х22 > 0,000933,-^ > 0,00675,> 0,117,-^ < 23,08 х18 х20 х58 х22

V = -24,527 + 0,0011х18 - 0,00013х20 - 0,0323х8 + 0,0254х3 > 0,000933,> 0,00675,< 0,117^-^^ > 23,08 х18 х20 х58 х22

V = -24,527 + 0,00 Щ8 -0,00013x^-0,0323.^ +0,5857л22 х х х х > 0,000933,^ > 0,00675,—^ < 0,117,^ < 23,08 х18 х20 х58 х22

у = -24,527 + 0,001Ц8 -0,0196х5 -0,0038х58 +0,0254.^ х х х х > 0,000933,^ < 0,00675, > 0,117,^ > 23,08 х18 х20 х58 х22

у = -24,527 + 0,0011х18 -0,0196Х5 -0,0038Х58 +0,5857Х22 х х х х > 0,000933,^ < 0,00675, > 0,117,^ < 23,08 х18 х20 х58 х22

V = -24,527 + 0,0011х18 - 0,0196х5 - 0,0323х8 + 0,0254х3 > 0,000933,-^ < 0,00675, < 0,117,-^ > 23,08 х18 х20 х58 х22

у = -24,527 + 0,0011.^ -0,0196х5 -0,0323х8 +0,5857х22 > 0,000933,-^ < 0,00675, < 0,117,-^ < 23,08 х18 х20 х58 х22

у = —24,527 +1,1 895Х2 -0,00013Х20 -0,0038Х58 + 0,0254Х3 х х х х < 0,000933,^ > 0,00675,—^ > 0,117,^ > 23,08 х18 х20 х58 х22

V = -24,527 + 1,1 895,Т2 -0,00013,т20 -0,0038.Т58 +0,5857.т22 х х х х < 0,000933,^ > 0,00675, > 0,117,^ < 23,08 х18 х20 х58 х22

V = —24,527 +1,1 895Л2 -0,000 13Л20 -0,0323.^ +0,0254л3 х х х х < 0,000933,^ > 0,00675, < 0,117,^ > 23,08 х18 х20 х58 х22

V = —24,527 +1,1 895Х2 -0,00013Х20 -0,0323Х8 +0,5857Х22 < 0,000933,-^ > 0,00675, < 0,117,-^ < 23,08 х18 х20 х58 х22

у = -24,527 + 1,1895Х2 -0,0196Х5 -0,0038Х58 +0,0254Х3 < 0,000933,-^ < 0,00675, > 0,117,-^ > 23,08 х18 х20 х58 х22

V = -24,527 + 1Д895х2 -0,0196х5 -0,0038х58 +0,5857х22 х х х х < 0,000933,^ < 0,00675,—^ > 0,117,^ < 23,08 х18 х20 х58 х22

V = -24,527 + 1, 1895л2 -0,0196Л5 -0,0323.^ +0,0254л3 х х х х < 0,000933,^ < 0,00675, < 0,117,^ > 23,08 х18 х20 х58 х22

V =-24,527 + 1,1895.Т2 -0,0196.Т5 -0,0323.Т8 +0,5857.Т22 < 0,000933,-^ < 0,00675,< 0,117,-^ < 23,08 х18 х20 х58 х22

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

^ = 0,000933 , 2 = 0,00675, 2 = 0,117,

3 = 23,08 играют роль точек переключения для

следующих четырех автоматически сформированных показателей:

- отношение процента трудоспособного населения (х2) к числу предприятий и организаций (х18);

- отношение численности пенсионеров (х5) к кредиторской задолженности (х20);

- отношение числа собственных легковых автомобилей (х8) к тарифам на грузовые ж/д перевозки

( х58);

- отношение численности рабочей силы ( х3 ) к объему производства электроэнергии (х22).

Тогда справедлива следующая интерпретация. • Если значение показателя х2/х18 не меньше, чем 0,000933, то на отправление грузов ж/д транспортом влияет число предприятий и организаций х18, а процент трудоспособного населения х2 не влияет. При этом с увеличением числа предприятий и организаций х!8 на одну единицу (при неиз-

менных значениях остальных переменных) отправление грузов у возрастает в среднем на 0,0011 млн руб. А если значение показателя х2/х!8 меньше, чем 0,000933, то на отправление грузов влияет процент трудоспособного населения х2, а число предприятий и организаций х18 не влияет. При этом с увеличением процента трудоспособного населения х2 на одну единицу (при неизменных значениях остальных переменных) отправление грузов у возрастает в среднем на 1,1895 млн руб.

• Если значение показателя х5/х20 не меньше, чем 0,00675, то на отправление грузов ж/д транспортом влияет кредиторская задолженность х20, а численность пенсионеров х5 не влияет. При этом с увеличением кредиторской задолженности х20 на 1 млн руб. (при неизменных значениях остальных переменных) отправление грузов у убывает в среднем на 0,00013 млн руб. А если значение показателя х5/х20 меньше, чем 0,00675, то на отправление грузов влияет численность пенсионеров х5, а кредиторская задолженность х20 не влияет. При этом с увеличением численности пенсионеров х5 на 1 тыс. человек (при неизменных значениях остальных переменных) отправление грузов у убывает в среднем на 0,0196 млн руб.

• Если значение показателя х8/х58 не меньше, чем 0,117, то на отправление грузов ж/д транспортом влияют тарифы на грузовые ж/д перевозки х58, а число собственных легковых автомобилей х8 не влияет. При этом с увеличением тарифов на грузовые ж/д перевозки х58 на одну условную единицу (при неизменных значениях остальных переменных) отправление грузов у убывает в среднем на 0,0038 млн руб. А если значение показателя х8/х58 меньше, чем 0,117, то на отправление грузов влияет число собственных легковых автомобилей х8, а тарифы на грузовые ж/д перевозки х58 не влияют. При этом с увеличением числа собственных легковых автомобилей на 1000 человек населения х8 на 1 шт. (при неизменных значениях остальных переменных) отправление грузов у убывает в среднем на 0,0323 млн руб.

• Если значение показателя х3/х22 не меньше, чем 23,08, то на отправление грузов ж/д транспортом влияет численность рабочей силы х3, а объемы производства электроэнергии х22 не влияют. При этом с увеличением численности рабочей силы х3 на 1 тыс. человек (при неизменных значениях остальных переменных) отправление грузов у возрастает в среднем на 0,0254 млн руб. А если значение показателя х3/х22 меньше, чем 23,08, то на отправление грузов влияют объемы производства электроэнергии х22, а численность рабочей силы х3 не влияет. При этом с увеличением объемов про-

изводства электроэнергии х22 на 1 млрд кВтч (при неизменных значениях остальных переменных) отправление грузов у возрастает в среднем на 0,5857 млн руб.

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

Также отметим, что сама по себе НЛР лучше подходит для моделирования в условиях мульти-коллинеарности, чем традиционная линейная регрессия. Дело в том, что чем больше в НЛР бинарных операций, тем выше число ее степеней свободы по сравнению с линейной регрессией. Это означает, что НЛР позволяет «вместить» в себя больше переменных при меньшем числе регрессо-ров, чем линейная регрессия. Например, НЛР (48) содержит только четыре регрессора, но при этом восемь переменных, поэтому шанс на то, что в ней будет присутствовать мультиколлинеарность, априори ниже, чем шанс столкнуться с мультикол-линеарностью в линейной регрессии со всеми восемью переменными.

ЗАКЛЮЧЕНИЕ

В работе рассмотрена НЛР, в которой помимо бинарной операции минимум используется еще и бинарная операция максимум. Предложен метод построения НЛР, в основе которого лежит решение задачи ЧБЛП. В результате решения этой задачи определяется структурная спецификация НЛР и ее приближенные МНК-оценки. Показано, как с помощью регулирования ограничений на бинарные переменные можно контролировать структурную спецификацию НЛР. Продемонстрировано, каким образом нужно исключать на начальном этапе противоречивые переменные, чтобы уменьшить время решения задачи и гарантировать интерпретируемость НЛР. С помощью предложенного метода построена модель, позволившая выявить новые закономерности функционирования железнодорожного транспорта в Иркутской области,

недоступные при использовании классической линейной регрессии.

Предложенный в работе метод является универсальным и может применяться для построения НЛР в любой предметной области при наличии статистических данных, состоящих только из положительных переменных. Проводимое в его рамках разбиение параметров приводит к формированию задачи ЧБЛП, оптимальное решение которой при достаточно большом числе разбиений дает оценки, практически не отличающиеся от оптимальных МНК-оценок рассматриваемой НЛР. Естественно, что с ростом числа разбиений будет расти и время решения задачи. Тем не менее, как продемонстрировано в работах [20, 21] на примере линейной регрессии, такая задача ЧБЛП решается на порядок быстрее, чем при использовании стандартных переборных процедур. Тестирование скорости построения НЛР по выборкам разных объемов с помощью предложенного метода будет проведено в последующих работах автора.

ЛИТЕРАТУРА

1. Arkes, J. Regression Analysis: A Practical Introduction. -Routledge, 2019. - 362 p.

2. Westfall, P.H., Arias, A.L. Understanding Regression Analysis: A Conditional Distribution Approach. - Chapman and Hall/CRC, 2020. - 514 p.

3. Клейнер Г.Б. Производственные функции: Теория, методы, применение. - М.: Финансы и статистика, 1986. - 239 с. [Kleiner, G.B. Proizvodstvennye funktsii: Teoriya, metody, primenenie. - Moscow: Finansy i statistika, 1986. - 239 s. (In Russian)]

4. Onalan, O., Basegmez, H. Estimation of Economic Growth Using Grey Cobb-Douglas Production Function: An Application for US Economy // Journal of Business Economics and Finance. - 2018. - Vol. 7, no. 2. - P. 178-190.

5. Yankovyi, O., Koval, V., Lazorenko, L., et al. Modeling Sustainable Economic Development Using Production Functions // Studies of Applied Economics. - 2021. - Vol. 39, no. 5.

6. Ishikawa, A. Why Does Production Function Take the Cobb-Douglas Form? // Statistical Properties in Firms' Large-scale Data. - Springer, Singapore, 2021. - P. 113-135.

7. Носков С.И. Технология моделирования объектов с нестабильным функционированием и неопределенностью в данных. - Иркутск: РИЦ ГП «Облинформпечать», 1996. -320 с. [Noskov, S.I. Tekhnologiya modelirovaniya ob"ektov s nestabil'nym funktsionirovaniem i neopredelennost'yu v dannykh. - Irkutsk: RITS GP «Oblinformpechat'», 1996. - 320 s. (In Russian)]

8. Шор Н.З. Методы минимизации недифференцируемых функций и их приложения. - Киев: Наук. думка, 1979. -200 с. [Shor, N.Z. Metody minimizatsii nedifferentsiruemykh funktsii i ikh prilozheniya. - Kiev: Nauk. dumka, 1979. - 200 s. (In Russian)]

9. Scaman, K., Bach, F., Bubeck, S., et al. Optimal Algorithms for Non-smooth Distributed Optimization in Networks // Advances in Neural Information Processing Systems. - 2018. - Vol. 31.

10.Khamaru, K., Wainwright, M. J. Convergence Guaranl Class of Non-convex and Non-smooth Optimization ] // Journal of Machine Learning Research. - 2019. - Vo 154. - P. 1-52.

11.Иванова Н.К., Лебедева С.А., Носков С.И. Иденти параметров некоторых негладких регрессий // Инф онные технологии и проблемы математического м вания сложных систем. - 2016. - № 17. - С. [Ivanova, N.K., Lebedeva, S.A., Noskov, S.I. Identifikat ametrov nekotorykh negladkikh regressii // Informa tekhnologii i problemy matematicheskogo modeli slozhnykh sistem. - 2016. - No. 17. - P. 107-110. sian)]

12.Носков С.И., Хоняков А.А. Программный комплекс ения некоторых типов кусочно-линейных регресси формационные технологии и математическое мод( ние в управлении сложными системами. - 2019. - № С. 47-55. [Noskov, S.I., Khonyakov, A.A. Programmr pleks postroeniya nekotorykh tipov kusochno-linei gressii // Informa-tsionnye tekhnologii i matematiches elirovanie v upravlenii slozhnymi sistemami. - 2019. (4). - P. 47-55. (In Russian)]

13. Park, Y.W., Klabjan, D. Subset Selection for Multipl Regression via Optimization // Journal of Global Opti

- 2020. - Vol. 77. - P. 543-574.

14. Chung, S., Park, Y.W., Cheong, T. A Mathematical 1 ming Approach for Integrated Multiple Linear Regress set Selection and Validation // Pattern Recognition. -Vol. 108. - P. 107565.

15. Bertsimas, D., Li, M.L. Scalable Holistic Linear Regr Operations Research Letters. - 2020. - Vol. 48, no. 3. -208.

16.Базилевский М.П. МНК-оценивание параметров с цированных на основе функций Леонтьева двухфа моделей регрессии // Южно-Сибирский научный ве 2019. - № 2 (26). - С. 66-70. [Bazilevskii, M.P otsenivanie parametrov spetsifitsirovannykh na osnove Leont'eva dvukhfaktornykh modelei regressii // Sibirskii nauchnyi vestnik. - 2019. - No. 2 (26). - P. 6< Russian)]

17.Базилевский М.П. Оценивание линейно-неэлеме: регрессионных моделей с помощью метода наш квадратов // Моделирование, оптимизация и инфор ные технологии. - 2020. - Т. 8. - № 4 (31). [Bazilevs Otsenivanie lineino-neehlementarnykh regressionnykh s pomoshch'yu metoda naimen'shikh kvadratov // M vanie, optimizatsiya i informatsionnye tekhnologii. -Vol. 8, - no. 4 (31). (In Russian)]

18.Базилевский М.П. Отбор информативных опера1 построении линейно-неэлементарных регрессион делей // International Journal of Open Information T gies. - 2021. - Т. 9. - № 5. - С. 30-35. [Bazilevskii, bor informativnykh operatsii pri postroenii neehlementarnykh regressionnykh modelei // Inte Journal of Open Information Technologies. - 2021. -no. 5. - P. 30-35. (In Russian)]

19.Базилевский М.П. Сведение задачи отбора информ регрессоров при оценивании линейной регрессион дели по методу наименьших квадратов к задаче ч; булевого линейного программирования // Модели оптимизация и информационные технологии. - 201

- № 1 (20). - С. 108-117. [Bazilevskii, M.P. Svedenii

ees for a Problems l. 20, no.

фикация ормаци-делиро-107-110. siya par-tsionnye rovaniya (In Rus-

построй // Ин-елирова-3 (4). -yi kom-nykh re-oe mod- No. 3

e Linear mization.

rogram-ion Sub-2020. -

ession // P. 203-

пецифи-кторных стник. -. MNK-funktsii Yuzhno-6-70. (In

нтарных еньших мацион-kii, M.P. modelei odeliro-2020. -

ий при ных мо-echnolo-M.P. Ot-lineino-rnational Vol. 9,

ативных ной мо-астично-рование, 8. - Т. 6. zadachi

otbora informativnykh regressorov pri otsenivanii lineinoi re-gressionnoi modeli po metodu naimen'shikh kvadratov k zadache chastichno-bulevogo lineinogo programmirovaniya // Modelirovanie, optimizatsiya i informatsionnye tekhnologii. -2018. - Vol. 6, no. 1 (20). - P. 108-117. (In Russian)] 20..Базилевский М.П. Способ определения параметра M в задаче частично-булевого линейного программирования для отбора регрессоров в линейной регрессии // Вестник Технологического университета. - 2022. - Т. 25. - № 2. - С. 62-66. [Bazilevskii, M.P. Sposob opredeleniya parametra M v zadache chastichno-bulevogo lineinogo programmirovaniya dlya otbora regressorov v lineinoi regressii // Vestnik Tekhno-logicheskogo universiteta. - 2022. - Vol. 25, no. 2. - P. 62-66. (In Russian)]

21.Konno, H., Yamamoto, R. Choosing the Best Set of Variables in Regression Analysis Using Integer Programming // Journal of Global Optimization. - 2009. - Vol. 44. - P. 273- 282.

Статья представлена к публикации руководителем регионального редсовета М.И. Гераськиным.

Поступила в редакцию 23.04.2022, после доработки 3.08.2022. Принята к публикации 31.08.2022.

Базилевский Михаил Павлович - канд. техн. наук., Иркутский государственный университет путей сообщения, г. Иркутск, Н mik2178@yandex.ru.

A METHOD FOR CONSTRUCTING NONELEMENTARY LINEAR REGRESSIONS BASED ON MATHEMATICAL PROGRAMMING

M.P. Bazilevskiy

Irkutsk State Transport University, Irkutsk, Russia H mik2178@yandex.ru

Abstract. This paper is devoted to constructing nonelementary linear regressions consisting of explanatory variables and all possible combinations of their pairs transformed using binary minimum and maximum operations. Such models are formalized through a 0-1 mixed integer linear programming problem. By adjusting the constraints on binary variables, we control the structural specification of a nonelementary linear regression, namely, the number of regressors, their types, and the composition of explanatory variables. In this case, the model parameters are approximately estimated using the ordinary least squares method. The formulated problem has advantages: the number of constraints does not depend on the sample size, and the signs of the estimates for the explanatory variables are consistent with the signs of their correlation coefficients with the dependent variable. Regressors are eliminated at the initial stage to reduce the time for solving the problem and make the model quite interpretable. A nonelementary linear regression of rail freight in Irkutsk oblast is constructed, and its interpretation is given.

Keywords: nonelementary linear regression, ordinary least squares method, 0-1 mixed integer linear programming problem, subset selection, coefficient of determination, interpretation, rail freight.