Iterative convex estimation of linear regression models under data stochastic heterogeneity

Full Text

Введение

Линейная регрессия представляет собой одну из наиболее распространенных вероятностных моделей, используемых для анализа зависимостей в различных областях знания [1]. Данная модель описывает линейную взаимосвязь между условным математическим ожиданием зависимой переменной $Y$ и набором объясняющих переменных $X_1$, $X_2$, $\dots$, $X_m$ на основе экспериментальных данных $(1, x_{i2}, \dots, x_{im}; y_i)$, где $i = 1, 2, \dots, n$. Формально это выражается соотношением [2]
\[\begin{equation}
    \operatorname{E}(Y|X_2=x_2,\dots,X_m=x_m) = \alpha_1 + \sum_{k=2}^m \alpha_k x_k.
\tag{1}
\end{equation}\]

Таким образом, модель (1) позволяет аппроксимировать величины, которые могут быть представлены в виде
\[
    \mathbf{y} = \mathbf{X}\boldsymbol{\alpha} + \boldsymbol{\varepsilon},
\]
где
\[
\mathbf{y} = \begin{pmatrix}
    y_1 \\
    y_2 \\
    \vdots \\
    y_n
\end{pmatrix}, \;\;
\mathbf{X} = \begin{pmatrix}
    1 & x_{12} & \cdots & x_{1m} \\
    1 & x_{22} & \cdots & x_{2m} \\
    \vdots & \vdots & \ddots & \vdots \\
    1 & x_{n2} & \cdots & x_{nm}
\end{pmatrix} = \begin{pmatrix}
    \mathbf{x}^{\top}_1 \\
    \mathbf{x}^{\top}_2 \\
    \vdots \\
    \mathbf{x}^{\top}_n
\end{pmatrix},\;\;
\boldsymbol{\alpha} = \begin{pmatrix}
    \alpha_1 \\
    \alpha_2 \\
    \vdots \\
    \alpha_m
\end{pmatrix},\;\;
\boldsymbol{\varepsilon} = \begin{pmatrix}
    \varepsilon_1 \\
    \varepsilon_2 \\
    \vdots \\
    \varepsilon_n
\end{pmatrix}.
\]
Здесь $\mathbf{y}$ — вектор наблюдений зависимой переменной; $\mathbf{X} = \{ x_{ij} \}_{n \times m}$ — матрица регрессоров, содержащая значения объясняющих переменных $X_1$, $X_2$, $\dots$, $X_m$, причем $\mathbf{x}^{\top}_i = (1, x_{i2}, \dots, x_{im})$; $\boldsymbol{\alpha}$ —  вектор неизвестных параметров модели; $\boldsymbol{\varepsilon}$ — вектор случайных ошибок наблюдений. 

Задача оценивания коэффициентов модели (1) в общем случае формулируется как задача оптимизации:
\[\begin{equation}
Q(\mathbf{a}) = \sum_{i=1}^{n} \rho\biggl( \biggl| y_i - \sum_{j=1}^{m} a_j x_{ij} \biggr| \biggr) \to \min_{\mathbf{a} \in \mathbb{R}^m},
\tag{2}
\end{equation}\]
где $\rho(\cdot)$ — функция потерь, определенная на положительной полуоси, удовлетворяющая условиям $\rho(0)=0$ и $\rho'(x) > 0$ для всех $x > 0$.

Классический линейный регрессионный анализ опирается на систему предпосылок, которые определяют требования к объясняющим переменным $X_i$, параметрам модели $\boldsymbol{\alpha}$ и случайным отклонениям $\boldsymbol{\varepsilon}$, а также позволяют исследовать статистические свойства оценок коэффициентов регрессии [3]:

• $1^\circ$) на вектор параметров $\boldsymbol{\alpha}$ не наложено ограничений, то есть $\mathcal{A} = \mathbb{R}^m$, где $\mathcal{A}$ "--- множество допустимых (априорных) значений параметров;
• $2^\circ$) матрица  $\mathbf{X}$ детерминирована, т.е. элементы $x_{ij}$ не являются случайными величинами;
• $3^\circ$) матрица $\mathbf{X}$ имеет полный столбцовый ранг: $\operatorname{rank}(\mathbf{X}) = m < n$;
• $4^\circ$) вектор ошибок $\boldsymbol{\varepsilon}$ случаен, следовательно, вектор наблюдений $\mathbf{y}$ также является случайным вектором;
• $5^\circ$) ошибки имеют нулевое математическое ожидание: $\operatorname{E}(\varepsilon_i) = 0$ для всех $i = 1, 2,\dots,n$, или в векторной форме $\operatorname{E}(\boldsymbol{\varepsilon}) = \mathbf{0}$;
• $6^\circ$) ошибки гомоскедастичны и некоррелированы: $\operatorname{cov}(\varepsilon_i,\varepsilon_k) = 0$ при $i \neq k$, $\operatorname{E}(\varepsilon_i^2) = \sigma^2_\varepsilon$ для всех $i,k = 1, 2,\dots,n$, где $\sigma^2_\varepsilon$ — постоянная дисперсия ошибок;
• $7^\circ$) распределение ошибок $\boldsymbol{\varepsilon}$ является нормальным (гауссовым).

При выполнении указанных предпосылок метод наименьших квадратов (МНК) обеспечивает получение эффективных оценок $\boldsymbol{\tilde{\alpha}}$ коэффициентов регрессии [4]. В рамках МНК в (2) функция потерь имеет вид $\rho(x)=\rho_{ls}(x)=x^2$. Статистические свойства МНК-оценок параметров модели и уравнения регрессии в целом хорошо исследованы. Однако в условиях стохастической неоднородности данных возможны нарушения исходных предположений, приводящие к потере устойчивости МНК-оценок.

Во-первых, эмпирические распределения ошибок часто демонстрируют более вытянутые хвосты по сравнению с нормальным распределением, что снижает эффективность МНК-оценок [5–7]. Особенно критичной для МНК является ситуация наличия выбросов — редких, часто асимметричных наблюдений, а также случаев, когда в модель включены лаговые эндогенные переменные [8–10]. Подобные аномалии могут возникать, например, при развитии процессов деградации в механических системах [11].

Во-вторых, в экономических, медицинских и других исследованиях гетерогенность данных часто приводит к появлению наблюдений, которые могут быть классифицированы как выбросы относительно основной выборки [12–14].

В-третьих, существенную проблему представляет гетероскедастичность — непостоянство дисперсии ошибок наблюдений [3, 15].

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

Если в (2) задать $\rho(x) = \rho_{ld}(x) = |x|$, то получим метод наименьших модулей (МНМ) [16, 17]. Для данного метода разработаны эффективные вычислительные алгоритмы [18–17]. 
Однако при одновременном воздействии нескольких факторов неоднородности, упомянутых выше, МНМ-оценки также теряют устойчивость [9].

В случае, когда в (2) функция потерь $\rho(x) = \rho_{gld}(x)$ удовлетворяет условию 
\[
\frac{d^2 \rho_{gld}(x)}{dx^2} < 0  \text{  для всех  }  x > 0,
\]
приходим к обобщенному методу наименьших модулей (ОМНМ) [21].

Применение вогнутых ($\rho'(x) > 0$, $\rho''(x) < 0$ $\forall x > 0$) и выпукло-вогнутых ($\rho'(x) > 0$ $\forall x > 0$ с точкой перегиба $c > 0$: $\rho''(x) > 0$ при $x < c$ и $\rho''(x) < 0$ при $x > c$) функций потерь приводит к появлению множества локальных минимумов в задаче (2), что существенно увеличивает вычислительную сложность процедуры оценивания. Одним из подходов к решению этой проблемы является сведение исходной невыпуклой задачи робастного оценивания к последовательности задач минимизации взвешенным методом наименьших квадратов (ВМНК) [3, 22] или взвешенным методом наименьших модулей (ВМНМ) [23, 24].

Данный подход был впервые использован для реализации МНМ через ВМНК и известен как метод Вейсфельда [25], или  метод вариационно-взвешенных квадратических приближений [16]. Следует отметить, что в работе [25] метод был эвристически применен к задачам размещения, тогда как в [16] он получил строгое теоретическое обоснование и детальное исследование.

Задачи минимизации в рамках ВМНК и ВМНМ характеризуются целевыми функциями вида
\[\begin{equation}
    V_*(\mathbf{a}) = \sum_{i=1}^{n} p_i \rho_*\biggl( \biggl| y_i - \sum_{j=1}^{m} a_j x_{ij} \biggr| \biggr) \to \min_{\mathbf{a} \in \mathbb{R}^m},
\tag{3}
\end{equation}\]
где $\rho_{*}(\cdot)$ представляет функцию потерь ($\rho_*(x) = \rho_{ls}(x) = x^2$ для МНК, $\rho_*(x) = \rho_{ld}(x) = |x|$ для МНМ), а $p_i \geqslant 0$ — весовые коэффициенты, отражающие точность соответствующих измерений. Эти коэффициенты призваны компенсировать расхождения между эмпирическим распределением ошибок $f_n(x)$ и их теоретически оптимальным распределением:

• нормальным распределением $f_\varepsilon(x) = \frac{1}{\sigma_\varepsilon \sqrt{2\pi}} \exp\bigl(-\frac{x^2}{2\sigma_\varepsilon^2}\bigr)$ для МНК;
• распределением Лапласа $f_\varepsilon(x) = \frac{1}{\sigma_\varepsilon \sqrt{2}} \exp\bigl(-\frac{\sqrt{2}|x|}{\sigma_\varepsilon}\bigr)$ для МНМ.

В силу выпуклости, непрерывности и ограниченности снизу целевых функций в (3) они обладают единственными точками глобального минимума.

Метод вариационно-взвешенных абсолютных приближений, предложенный в работах [23, 24], имеет ряд ограничений:

1. исследован только для случая вогнутых функций потерь $\rho(x)$;
2. не изучена область его применимости в зависимости от свойств $\rho(x)$;
3. отсутствует анализ вычислительной сложности итерационного процесса, требующего решения последовательности выпуклых задач вида (3) вместо исходной задачи (2).

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

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

1. Методы исследования

Рассмотрим итерационные процедуры взвешивания квадратов и модулей невязок.

1.1. Использование ВМНК. В условиях стохастической неоднородности данных при наличии асимметричных отклонений распределения ошибок от нормального закона целесообразно рассматривать функции потерь, удовлетворяющие следующим условиям [4]:

1. существует конечный положительный предел: 
      \[{0 < \lim\limits_{x \to +0}  \frac{\rho(x)}{x^2} = c_1 < \infty};\]
2. функция имеет горизонтальную асимптоту: \[0 < \lim\limits_{x \to +\infty} \rho(x) = d_1 < \infty.\]

Следует отметить, что реализация метода остается возможной и при отсутствии горизонтальной асимптоты у функции потерь при условии выполнения соотношения $\lim\limits_{x \to +\infty}  {\rho(x)}/{x^2} = 0$.

Итерационная процедура заключается в решении последовательности задач минимизации
\[\begin{equation}
    V_{ls}^{(k)}(\mathbf{a}) = \sum_{i=1}^{n} p_i^{(k)} \biggl( \biggl| y_i - \sum_{j=1}^{m} a_j x_{ij} \biggr| \biggr)  \to \min_{\mathbf{a} \in \mathbb{R}^m},
\tag{4}
\end{equation}\]
где весовые коэффициенты определяются так:
\[
p_i^{(k)} = \frac{\rho ( |e_i^{(k-1)} | )}{ (e_i^{(k-1)} )^2}, 
\quad
e_i^{(k-1)} = y_i - \sum_{j=1}^m a_j^{(k-1)} x_{ij}.
\]
Здесь $e_i^{(k-1)}$ — ошибки на $(k-1)$-й итерации; $a_j^{(k-1)}$ — ВМНК-оценки параметров на $(k-1)$-й итерации; $a_j^{(0)}$ — начальные МНК-оценки; $p_i^{(0)} = 1$ для всех $i = 1, 2,\dots,n$. В случае стремления ошибок $e_i^{(k-1)}$ к нулю весовые коэффициенты $p_i^{(k)}$ вычисляются с помощью предельного перехода в нулевой точке.

Оценки параметров в (4) вычисляются по формуле
\[\begin{equation}
\mathbf{a} = (\mathbf{U}^{\top} \mathbf{U})^{-1} (\mathbf{U}^{\top} \mathbf{v}),
\tag{5}
\end{equation}\]
где $\mathbf{U} = \mathbf{W}^{-1} \mathbf{X}$, $\mathbf{v} = \mathbf{W}^{-1} \mathbf{y}$, $\mathbf{W} = \operatorname{diag} (p_1^{(k)},\dots,p_n^{(k)} )$.

Рассмотрим три классические функции потерь, предложенные Эндрюсом, Мешалкиным и Рамсеем [4]. Для унификации анализа введем единый параметр $\lambda$, характеризующий их свойства:

• $\lambda \to 0$ все три функции стремятся к квадратичной $\rho(x) = x^2$;
• при $\lambda \neq 0$ функции обладают горизонтальной асимптотой.

Графическая иллюстрация поведения функций потерь при значениях параметра $\lambda=0.1$ и 0.5 представлена на рис. 1.

Рис. 1. Графики функций потерь $\rho(x)$ при $\lambda=0.1$ и $0.5$
[Figure 1. Plots of loss functions $\rho(x)$ at $\lambda=0.1$, and $0.5$]

Функция потерь Эндрюса задается кусочно-определенным выражением
\[\begin{equation}
\rho_A(x) =
    \begin{cases} 
(2\lambda)^{-1} \bigl[1 - \cos (x \sqrt{2\lambda} )\bigr], & |x| < \pi/\sqrt{2\lambda}, \\
        \lambda^{-1}, & |x| \geqslant \pi/\sqrt{2\lambda}.
    \end{cases}
\tag{6}
\end{equation}\]
Соответствующие весовые коэффициенты вычисляются по формуле
\[
    p_i^{(k)} = \frac{1 - \cos (e_i^{(k-1)} \sqrt{2\lambda} )}{2\lambda  (e_i^{(k-1)})^2}.
\]
Для корректной программной реализации алгоритма необходимо отдельно рассматривать случай нулевых невязок. Предельное значение весового коэффициента при $e_i^{(k-1)} \to 0$ определяется следующим образом:
\[\begin{equation}
\lim_{x \to 0} \frac{1 - \cos (x \sqrt{2\lambda} )}{2\lambda x^2} =   \frac{1}{2}.
\end{equation}\]
Таким образом, в вычислительной реализации следует полагать $p_i^{(k)} = 1/2$ при $e_i^{(k-1)} = 0$.

Функция потерь Мешалкина определяется выражением
\[\begin{equation}
    \rho_M(x) = \frac{1 - \exp(-{\lambda x^2}/{2})}{\lambda}.
\tag{7}
\end{equation}\]

Соответствующие весовые коэффициенты вычисляются так:
\[
p_i^{(k)} = \frac{1 - \exp(-{\lambda (e_i^{(k-1)})^2}/{2})}{\lambda (e_i^{(k-1)})^2}.
\]
Предельное значение при нулевых невязках определяется следующим образом:
\[\begin{equation}
    \lim_{x \to 0} \frac{1 - \exp (- {\lambda x^2}/{2})}{\lambda x^2} = \frac{1}{2}.
\end{equation}\]

Функция потерь Рамсея задается формулой
\[\begin{equation}
\rho_R(x) = \frac{1 - (1 + \sqrt{\lambda}|x|)\exp (-\sqrt{\lambda}|x|)}{\lambda}.
\tag{8}
\end{equation}\]
Весовые коэффициенты для данной функции следующие:
\[
p_i^{(k)} = \frac{1 -  (1 + \sqrt{\lambda} |e_i^{(k-1)} | )\exp (-\sqrt{\lambda} |e_i^{(k-1)} |)}{2 \lambda (e_i^{(k-1)})^2}.
\]
Предельное значение:
\[
\lim_{x \to 0} \frac{1 -  ( 1 + \sqrt{\lambda} |x|  ) \exp ( -\sqrt{\lambda} |x|  )}{\lambda x^2} = 
  \frac{1}{2}.
\]

Таким образом, алгоритм итерационно пересчитывает значения весовых коэффициентов до тех пор, пока разница между весами текущей и предыдущей итераций не станет минимальной. Экспериментальные исследования показали, что при достижении точности порядка $10^{-12}$ скорость сходимости метода существенно снижается. Практически оптимальным компромиссом между точностью и вычислительной эффективностью является выполнение двух условий: максимальная допустимая разность весов должна составлять $10^{-8}$,  а максимальное число итераций ограничивается 50 циклами.

Для обеспечения единого подхода к анализу различных функций потерь определим  в уравнениях (6)–(8) значения параметра $\lambda$, соответствующие точке перегиба при $x=3\sigma$ (рис. 2). Численные расчеты при $\sigma=1$ дают следующие значения:

• $\lambda \approx 0.137$ для функции Эндрюса (6);
• $\lambda \approx 0.111$ для функций Мешалкина (7) и Рамсея (8). 

Рис. 2. Графики вторых производных функций потерь в точках перегиба $x =3 \sigma$
[Figure 2. Graphs of second derivatives of loss functions at inflection points $x =3 \sigma$]

Отметим, что решение задачи регрессионного оценивания методом взвешенных наименьших квадратов остается возможным даже при условии $\lim\limits_{x \to +0} {\rho(x)}/{x^2} = 0$. Однако в этом случае распределение ошибок будет иметь более короткие хвосты по сравнению с нормальным распределением, асимптотически приближаясь к равномерному распределению. Следует подчеркнуть, что такая ситуация является скорее теоретической, поскольку на практике встречается крайне редко.

В противоположном случае, когда $\lim\limits_{x \to +0}  {\rho(x)}/{x^2} = \infty$, возникает принципиально иная ситуация: весовые коэффициенты, соответствующие малым ошибкам, неограниченно возрастают. Это приводит к неустойчивости оценок параметров регрессии, делая метод неприменимым на практике.

1.2. Использование ВМНМ. Здесь возможны два варианта.

Первый вариант представляет собой итерационную реализацию ОМНМ с функциями потерь $\rho_{gld}(x)$, удовлетворяющими в нуле условию
\[
0 < \lim\limits_{x \to +0} \frac{\rho(x)} {x} = c_2 < \infty.
\]
В качестве конкретных примеров таких функций потерь рассмотрим $\operatorname{arctg}(|x|)$ и $\ln(|x| + 1)$.  Для обеих этих функций выполняется условие
\[\begin{equation}
\lim_{x \to +0} \frac{\rho_{gld}(x)}{|x|} = 1.
\end{equation}\]

Второй вариант предполагает использование ВМНМ для итерационного решения задачи (2) с функциями потерь, обладающими следующими свойствами:
\[\begin{equation}
\lim_{x \to +0} \frac{\rho(x)}{x} = 0, \quad \lim_{x \to \infty} \frac{\rho(x)}{x} < \infty.
\end{equation}\]
Конкретные примеры таких функций потерь $\rho_{ls}(x)$ приведены в (6)–(8).

Важно отметить, что задачу (2) с функциями потерь, для которых
\[\begin{equation}
\lim_{x \to +0} \frac{\rho(x)}{x} = \infty,
\end{equation}\]
невозможно свести к последовательности задач ВМНМ-минимизации. Как и в п. 1.1, при малых ошибках весовые коэффициенты в таком случае неограниченно возрастают, что делает невозможным получение устойчивых оценок. В частности, это относится к функциям потерь вида $\rho(x) = |x|^\beta$ при $0 < \beta < 1$.

Вместо исходной задачи (1) итерационно решается последовательность задач ВМНМ:
\[\begin{equation}
V_{ld}^{(k)}(\mathbf{a}) = \sum_{i=1}^{n} p_i^{(k)} 
\biggl| y_i - \sum_{j=1}^{m} a_j x_{ij} \biggr| \to \min_{\mathbf{a} \in \mathbb{R}^m},
\tag{9}
\end{equation}\]
где весовые коэффициенты $p_i^{(k)}$ определяются так:
\[
    p_i^{(k)} = \frac{\rho(|e_i^{(k-1)}|)}{|e_i^{(k-1)}|},
\]
а ошибки $e_i^{(k-1)}$ на $(k-1)$-й итерации вычисляются по формуле
\[
    e_i^{(k-1)} = y_i - \sum_{j=1}^{m} a_j^{(k-1)} x_{ij}.
\]
Здесь $a_j^{(k-1)}$ — оценки параметров на $(k-1)$-й итерации, $a_j^{(0)}$ — МНМ"=оценки, $p_i^{(0)} = 1$ для всех $i = 1, \dots, n$.

Аналогично предыдущему случаю (см. п. 1.1), регрессионное оценивание с помощью ВМНМ остается корректным и при $e_i^{(k-1)} \to 0$ благодаря предельному переходу в нуле.

В работах [23, 24] предложен альтернативный подход, где весовые коэффициенты определяются так:
\[
    p_i^{(k)} = \rho'(|e_i^{(k-1)}| ).
\]
Однако следует отметить важное отличие: в общем случае $\rho'(x) \neq \rho(x)/x$. Поэтому
\[
\lim_{k \to \infty} \Bigl( p_i^{(k)} - \frac{\rho (|e_i^{(k-1)}| )}{|e_i^{(k-1)}|} \Bigr) \neq 0,
\]
а последовательность оценок (9) реализует иную функцию потерь. Например, при использовании $\rho(x) = \operatorname{arctg}(|x|)$ фактически минимизируется функция $\tilde{\rho}(x) =  {|x|}/({x^2+1})$.

Для решения задачи ВМНМ-оценивания применим алгоритмы покоординатного [26] и модифицированного градиентного [20] спусков по узловым прямым. Прежде чем описать принцип работы алгоритмов, введем необходимые определения.

Пусть $ \boldsymbol{\Omega} = \{\Omega_1, \ldots, \Omega_n\}$ — множество всех $n$ гиперплоскостей \(\Omega_i\):
\[
\Omega_i: {y_i - \langle\mathbf{a}, \mathbf{x}_i\rangle = 0},\quad  i = 1, \dots, n,
\]
где каждая гиперплоскость $\Omega_i$ соответствует $i$-му наблюдению выборки.

Узловой точкой называется точка пересечения $m$ линейно независимых гиперплоскостей:
\[\begin{equation}
\mathbf{u} = \bigcap_{s \in M} \Omega_s, \quad M = \{k_1, \dots, k_m\}, \quad k_1 < k_2 < \dots < k_m, \quad k_l \in \{1, \dots, n\}.
\end{equation}\]
Множество всех узловых точек обозначим через $U$.

Узловой прямой называется прямая, являющаяся пересечением $(m-1)$ линейно независимых гиперплоскостей:
\[\begin{equation}
l_{(k_1, \dots, k_{m-1})} = \bigcap_{i \in \{k_1, \dots, k_{m-1}\}} \Omega_i, \quad k_l \in \{1, \dots, n\}.
\end{equation}\] 

Алгоритм покоординатного спуска реализует следующий итерационный процесс.

1.  На нулевой итерации выбирается произвольная узловая точка $\mathbf{u}^{(0)}$, полученная решением СЛАУ, составленной из $m$ произвольных наблюдений.
2. Для текущей узловой точки $\mathbf{u}^{(k)}$ анализируются все инцидентные ей узловые прямые.
3. Находится узловая прямая, на которой целевая функция $V_{ld}^{(k)}(\mathbf{a})$ достигает минимального значения.
4. Осуществляется переход в новую узловую точку $\mathbf{u}^{(k+1)}$ вдоль найденной прямой.
5. Процесс продолжается до достижения точки $\mathbf{u}^{(*)}$, которая является локальным минимумом на множестве всех инцидентных ей узловых прямых.

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

Для каждой узловой точки $\mathbf{u} = (u_1,\ldots,u_m)$ производная по направлению узловой прямой вычисляется по формуле
\[
\frac{\partial V_{ld}^{(k)}(\mathbf{u})}{\partial \mathbf{l}_{(k_1, \dots, k_{m-1})}} = 
\sum_{i=1}^n \biggl( \bigl( c_1 + x_{i2} c_2 + \dots + x_{im} c_m \bigr) \cdot 
\operatorname{sign}\biggl( \sum_{j=1}^m u_j x_{ij} - y_i \biggr) \biggr),
\]
где $\mathbf{l}_{(k_1, \dots , k_{m-1})} = (c_1,c_2,\dots,c_m)$ — направляющий вектор узловой прямой $l_{(k_1, \dots  , k_{m-1})}$.

Критерий минимума формулируется следующим образом: если односторонние производные
\[
\frac{\partial V_{ld}^{(k)}(\mathbf{u})}{\partial \mathbf{l}_{(k_1, \dots, k_{m-1})}} \Bigr|_{\mathbf{u}_-} \quad \text{и} \quad 
\frac{\partial V_{ld}^{(k)}(\mathbf{u})}{\partial \mathbf{l}_{(k_1, \dots, k_{m-1})}} \Bigr|_{\mathbf{u}_+}
\]
в узловой точке $\mathbf{u}$ имеют разные знаки, то эта точка является минимумом функции $V_{ld}^{(k)} (\mathbf{a})$ на узловой прямой $l_{(k_1, \dots  , k_{m-1})}$.

Для повышения быстродействия ОМНМ-оценивания регрессии с вогнутыми функциями потерь $\rho_{gld} (x)$ вместо полного перебора всех узловых точек [21] можно использовать приближенные алгоритмы, такие как обобщенный спуск — итерационный алгоритм оценивания регрессии на основе ОМНМ [27]. В этом алгоритме уменьшается исследуемая окрестность точного решения и допускается определенная доля ошибок с заданной вероятностью (например, не более 0.05), т.е. сохраняется статистическая незначимость потери точности. 

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

2. Обсуждение результатов

Наличие выбросов в выборке может существенно ухудшить качество и надежность модели, что особенно критично для линейного регрессионного анализа. Для комплексного исследования этого эффекта мы проводим серию вычислительных экспериментов с различными типами засорений данных, используя метод статистических испытаний Монте—Карло [28]. В экспериментах генерируются выборки с числом параметров модели $m=2,\dots,5$.

Генерация засоренных случайных ошибок $\varepsilon$ выполняется по схеме Тьюки—Хьюбера [29, 30] с вероятностью засорения $\gamma=0.1$:
\[
F_\varepsilon(x) = (1 - \gamma)F_0(x) + \gamma F_H(x),
\]
где $F_0(x)$ — функция распределения основных ошибок (используется стандартное нормальное распределение); $F_H(x)$ — функция распределения засорений (рассматриваются двустороннее и одностороннее распределения Коши и Лапласа).

Нормальное распределение выбрано как эталон с быстро убывающими хвостами, распределение Коши — как крайний случай с вытянутыми хвостами, максимально подверженный влиянию выбросов, а распределение Лапласа занимает промежуточное положение между ними.

Вычислительные эксперименты проводились на персональном компьютере Dell G5 5587 со следующими характеристиками: 6-ядерный процессор Intel Core i7-8750H; тактовая частота до 4.1 ГГц. Программная реализация алгоритмов выполнена на языке C++ в среде разработки Microsoft Visual  Studio 2019.

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

2.1. Взвешенный метод наименьших квадратов. В результате проведения 100 вычислительных экспериментов при $n=50$, $100$, $\dots$, $300$ и $m=2$, $3,\dots,5$ было выявлено, что итерационный алгоритм ВМНК вне зависимости от функции потерь имеет сильную зависимость от вида засорений. Так, число итераций при добавлении односторонних или двусторонних засорений Коши возросло в среднем в 1.5—2 раза.

Сравнение по критериям Стьюдента и Манна--Уитни среднего времени работы алгоритма ВМНК для данных без засорения и с засорением с распределением Коши для функций потерь Эндрюса и Мешалкина отклонило нулевую гипотезу об однородности. Для функции Рамсея количество итераций и время работы алгоритма возросло в 1.2—2.5 раза, что связано с более медленным достижением асимптоты (рис. 1), это привело к выполнению нулевой гипотезы вне зависимости от вида засоряющего распределения.

Ввиду отсутствия значимой разницы между временем работы алгоритма при односторонних и двусторонних засорениях далее условно разделим их на 3 типа: распределения Гаусса, распределения с засорениями Коши и распределения с засорениями Лапласа. Для распределений Коши и Лапласа взяты усредненные значения вне зависимости от симметричности отклонений. Зависимости для среднего времени вычисления 1000 экспериментов от размерности выборки с коэффициентом детерминации $R^2=0.98$ представлены в табл. 1 и на рис. 3.

Таблица 1. Зависимости среднего времени вычисления 1000 вычислительных экспериментов от размерности выборки $n$, порядка модели $m$ и распределения засорений ошибок [Dependencies of the average computation time for 1000 simulation runs  on the sample size $n$, model order $m$, and error contamination distribution]

Loss Function	Gaussian	Laplace	Cauchy
Andrews	$0.00005 \cdot n^{2.93}$	$0.0001 \cdot n^{2.77} m^{0.16}$	$0.0002 \cdot n^{2.75} m^{0.15}$
Meshalkin	$0.00005 \cdot n^{2.8} m^{0.08}$	$0.0001 \cdot n^{2.78} m^{0.16}$	$0.0003 \cdot n^{2.71} m^{0.14}$
Ramsay	$0.0002 \cdot n^{2.82} m^{0.08}$	$0.0002 \cdot n^{2.78} m^{0.15}$	$0.0008 \cdot n^{2.54} m^{0.16}$

Из полученных зависимостей выделим три важных момента:

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

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

Рис. 3. Зависимость десятичного логарифма среднего времени обработки 1000 вычислительных экспериментов (сек.) от размера выборки $n$ при $m=3$: (a) для незасоренного распределения ошибок; (b) для распределения с засорениями по Коши
Figure 3. Decimal logarithm of the average processing time for 1000 simulation runs (seconds) as a function of sample size $n$ at $m=3$: (a) error distribution without contamination; (b) Cauchy-contaminated distribution]

Оценим вычислительную сложность ВМНК при вычислении коэффициентов регрессии по формуле (5). Основные вычислительные операции и их сложность:

• формирование матриц $\mathbf{U}$ и векторов $\mathbf{v}$: $O(n^2 m)$ и $O(n^2)$ операций соответственно;
• вычисление матричного произведения $\mathbf{U}^{\top} \mathbf{U}$: $O(n m^2)$ операций;
• вычисление обратных матриц: $O(m^3)$ операций;
• умножение на вектор $\mathbf{v}$: $O(n m)$ операций;
• определение вектора весов $\mathbf{p}$: $O(n m)$ операций;
• диагонализация вектора весов: $O(n)$ операций.

Таким образом, общая вычислительная сложность алгоритма составляет $O(n^2 m k)$ при $n \gg m$, где $k$ — количество итераций алгоритма.

Число итераций существенно зависит от способа засорения выборки, что свидетельствует о влиянии вида распределения и не позволяет использовать средние значения для однозначного определения соотношений для каждой функции потерь. Как показали исследования (см. табл. 1), наименьшее число итераций наблюдается для незасоренных выборок, а наибольшее — для выборок с засорениями по Коши. Это позволяет определить примерную область значений для числа операций в нотации Big O.

Экспериментальные данные демонстрируют слабую тенденцию роста числа итераций с увеличением размера выборки. В частности, для функции потерь Эндрюса получены следующие эмпирические зависимости:

• $k \approx 9.01 \cdot n^{0.01} m^{0.08}$ ($R^2=0.74$) при нормальном распределении;
• $k \approx 10.56 \cdot n^{0.08} m^{0.19}$ ($R^2=0.75$) при засорении по Коши.

Отметим слабую зависимость числа итераций $k$ от размера выборки и ее усиление при удлинении хвостов у распределения ошибок. Получены следующие характерные диапазоны числа итераций $k$:

• $10.5< k <19.8$ для функции Эндрюса; 
• $11.3< k<20.6 $ для функции Мешалкина;
• $25.5<k<27.6$ для функции Рамсея.

Поэтому на практике в качестве $k$ целесообразнее использовать средние арифметические значения.

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

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

Результаты показывают, что при $n=50$ для незасоренных данных и данных с засорением Лапласа СКО всех коэффициентов (кроме свободного члена) составляет 0.01 и уменьшается с ростом объема выборки. Это подтверждает высокую точность модели в описании зависимости между переменными. 
Однако при засорении Коши наблюдаются значительные колебания СКО — от 0.1 до нескольких сотен, что делает применение ВМНК в таких случаях нецелесообразным.

2.2. Взвешенный метод наименьших модулей. Согласно критериям Стьюдента и Манна—Уитни, вероятность неотклонения нулевой гипотезы об однородности выборок для среднего времени расчета 1000 экспериментов и среднего числа итераций 100 экспериментов при использовании функций потерь $\operatorname{arctg}(|x|)$ и $\ln(|x|+1)$ для алгоритмов спуска превышает 0.05. Следовательно, у ВМНМ, в отличие от ВМНК, нет зависимости от вида распределения анализируемой выборки и алгоритм использует схожее количество вычислительных ресурсов как для незасоренных данных, так и для выборок с засорениями Коши и Лапласа.

Зависимости среднего времени расчета 1000 экспериментов для трех алгоритмов (покоординатного, градиентного и обобщенного спусков) с $R^2=0.99$ представлены в табл. 2. В алгоритме обобщенного спуска доля ошибочных решений не превысила установленный порог в 5%, что статистически незначимо.

Таблица 2. Зависимости среднего времени вычисления 1000 вычислительных экспериментов алгоритмов спуска размерности выборки $n$ и порядка модели $m$
[Dependencies of the average computation time for 1000 simulation runs of descent algorithms on the sample size $n$, and model order $m$] 

$\rho(x)$ \ Algorithm	Coordinate Descent	Gradient Descent	Generalised Descent
$\ln(|x|+1)$	$0.0001 \cdot n^{2.11} m^{2.03}$	$0.0007 \cdot n^{1.3} m^{2.27}$	$7.8 \cdot 10^{-7} \cdot n^{1.4} m^{9.4}$
$\operatorname{arctg}(|x|)$	$0.0002 \cdot n^{2.0} m^{2.08}$	$0.0009 \cdot n^{1.27} m^{2.28}$	$3 \cdot 10^{-6} \cdot n^{1.4} m^{8.65}$

По сравнению с ВМНК в ВМНМ наблюдается значительно более сильная зависимость от числа параметров модели $m$, обусловленная экспоненциальным ростом количества узловых точек при увеличении размерности. Этот эффект особенно выражен в алгоритме обобщенного спуска, где требуется уточнение решения полным перебором. В то же время уменьшение влияния размера выборки $n$ компенсирует возросшую зависимость от $m$, так как для построения надежной регрессионной модели необходимо выполнение условия $n \gg m$.

Использование функции потерь $\ln(|x|+1)$ позволило сократить как время работы алгоритмов, так и среднее количество итераций. По сравнению с $\operatorname{arctg}(|x|)$ для градиентного и покоординатного спусков наблюдается выигрыш в 5–10% по времени вычислений, который снижается с увеличением размерности выборки $n$. В случае обобщенного спуска с логарифмической функцией потерь при $m=2$ достигается примерно двукратное ускорение, однако при росте $m$ разница в быстродействии  уменьшается до 1.5–5%.  

Для оценки вычислительной сложности алгоритмов покоординатного и градиентного спусков отметим, что основное изменение в алгоритмах связано с добавлением итерационного пересчета, и после первой итерации процесс спуска практически прекращается. Для оценки дополнительной сложности, вносимой итерациями, рассчитано количество анализируемых узловых точек  (${\approx m^{1.22}n^{0.77}}$) и переходов между узловыми прямыми ($\approx m^{0.26}/n^{0.18}$) на одной итерации. Согласно результатам работы [31], после уточнения сложность алгоритмов составляет $O(m^3 n^2+m^{1.97} n^{1.54} k)$ (для покоординатного спуска) и $O( m^3 n^2 \ln n + C_{\zeta n}^m (m^3 + mn))$ (для обобщенного спуска)  с областью решений, соответствующей 95%-му уровню доверительной вероятности $\zeta_{95}=(7.06 \cdot m^{0.57})/n^{0.93}$.

Алгоритм градиентного спуска включает три основных этапа: 

• $P_1$) поиск первого приближения в области из $48.3 \cdot m^{0.2}/n^{0.2}$ случайных наблюдений; 
• $P_2$) анализ полной выборки;
• $P_3$) итерационный пересчет. 

Число переходов между узловыми прямыми и  число рассматриваемых узловых точек на этапе $P_1$ составляет $m^{0.8}n^{0.2}$ и $m^{1.5}n^{0.5}$, а на этапе $P_2$ — $mn^{0.15}$ и $m^{1.8}n^{0.5}$.

При  $n \gg m$ вычислительная сложность этапов алгоритма такая:
\[
\begin{aligned}
        P_1 &= O\bigl(m^{1.8}n^{1.2} (m^2 + \ln n + n^{0.3}/m^{0.3} )\bigr) = O(m^{1.8}n^{1.5}), \\
        P_2 &= O\bigl(m^2n^{1.15} (m^2 + \ln n + n^{0.35}/m^{0.2} )\bigr) = O(m^{1.8}n^{1.5}), \\
         P_3 &= O\bigl(m^{1.26}n^{0.82}k (m^2 + \ln n + n^{0.72}/m^{0.29} )\bigr) = O(m^{0.97}n^{1.54}k).
\end{aligned}
\]

Общая вычислительная сложность алгоритма:
\[
    P = P_1 + P_2 + P_3 = O (m^{1.8}n^{1.5} + m^{0.97}n^{1.54}k ).
\]

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

Применение критериев Стьюдента и Манна–Уитни показывает, что хотя вероятность принятия нулевой гипотезы в отдельных случаях снижается до 50%, среднее число итераций остается статистически неразличимым для различных типов засорений. Для используемых функций потерь установлены следующие зависимости $k$ с коэффициентом детерминации $R^2=0.97$:

• $k \approx 1.87 \cdot n^{0.12} m^{0.17}$ (в среднем 4.3 итерации) для $\ln(|x|+1)$;
• $k \approx 2.4 \cdot n^{0.09} m^{0.16}$ (в среднем 4.6 итерации) для $\operatorname{arctg}(|x|)$. 

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

Оценка точности выполнена путем сравнения с решением, полученным методом обобщенного спуска. Согласно данным [27], оценки этого алгоритма отклоняются от результатов полного перебора узловых точек менее чем в 5% случаев, что статистически незначимо. В итоге решения, полученные алгоритмами ВМНМ, отличаются от точного в 100% случаев, однако эти отклонения меньше 0.5%  и уменьшаются с ростом размера выборки. 

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

Анализ изменения времени работы алгоритмов показывает, что введение итерационного пересчета весов привело к увеличению времени работы алгоритмов: 

• в 1.3–1.6 раза для  градиентного спуска; 
• в 1.4–2 раза для покоординатного спуска.

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

2.3. Сравнительный анализ. Сравнительный анализ рассмотренных алгоритмов показывает, что наивысшая точность достигается при использовании обобщенного метода наименьших модулей, что обусловлено его максимальным соответствием результатам полного перебора. Однако данный алгоритм демонстрирует сильную зависимость от числа параметров модели $m$ из-за необходимости уточнения решения полным перебором.

Алгоритмы ВМНК и ВМНМ, согласно анализу среднеквадратических отклонений, обеспечивают высокую точность и согласованность результатов. При этом ВМНК и покоординатный спуск более чувствительны к числу наблюдений $n$, что существенно, поскольку для качественного моделирования требуется выполнение условия $n \gg m$.  Дополнительно отметим, что в ВМНК количество итераций пересчета весовых коэффициентов в 2–6 раз превышает соответствующий показатель ВМНМ.

Для объективной оценки вычислительной трудоемкости обратимся к графикам сравнения времени работы алгоритмов (рис. 4), где представлены оптимальные по быстродействию варианты: ВМНК с функцией потерь Эндрюса (для незасоренных данных) и ВМНМ с логарифмической функцией потерь.

Рис. 4. Десятичный логарифм среднего времени вычислений (сек.) для 1000 вычислительных экспериментов алгоритмами ВМНК, ВМНМ и обобщенного спуска при (a) $m=3$ и (b) $n=150$
[Figure 4. Decimal logarithm of average computation time (seconds) for 1000 simulation runs using WLS, WLAD and generalized descent algorithms for cases: (a) $m=3$, and (b) $n=150$]

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

Сравнение точности алгоритмов ВМНМ и ВМНК по разнице СКО от истинных значений выявило, что  при $n=50$ и $m=2$ разница СКО составляет $0.006$ и уменьшается с ростом $m$ и $n$.

Полученные результаты свидетельствуют, что различия в точности оценок алгоритмов ВМНМ и ВМНК малы и малосущественны в контексте практического применения.

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

• ошибки первого рода («ложное срабатывание») — отвержение верной нулевой гипотезы об отсутствии связи между явлениями или искомого эффекта;
• ошибки второго рода («пропуск цели») — принятие неверной нулевой гипотезы.

МНК преимущественно подвержен ошибкам первого рода из-за чувствительности его оценок к нарушениям предпосылок $1^\circ$–$7^\circ$. A некорректное применение невыпуклых функций потерь может привести к ошибкам второго рода, когда игнорируются реальные изменения в данных.

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

Если нет оснований считать данные стохастически неоднородными, то МНК остается простым и надежным выбором для моделирования, при этом МНК- и МНМ-оценки будут близки. В случаях автокорреляции ошибок и пространственной гетероскедастичности (когда дисперсия ошибок коррелирует с объясняющими переменными) обычно достаточно применять обобщенный МНК, широко используемый в эконометрике [3].

При наличии засорений, не связанных с объясняющими переменными, рекомендуется использовать ВМНМ с функциями потерь, удовлетворяющими условиям $\lim\limits_{x \to +0}  {\rho(x)}/{x} = 0$ и $\lim\limits_{x \to \infty} {\rho(x)}/{x} < \infty$, либо ВМНК с ${\lim\limits_{x \to \infty} \rho(x) < \infty}$. Если засорения носят симметричный характер, подходит метод Хьюбера [30], где ${\lim\limits_{x \to \infty} {\rho(x)}/{x} < \infty}$, практические аспекты  реализации которого рассмотрены в [32].

Когда в модель включаются лаговые значения зависимой переменной, нарушается предпосылка $2^\circ$, так как эти переменные становятся случайными и коррелируют с ошибками $\boldsymbol{\varepsilon}$. В таких случаях при несимметричных засорениях требуются особо устойчивые оценки на основе ВМНМ с функциями потерь $\rho_{gld}(x)$.

Потребность в вогнутых и выпукло-вогнутых функциях потерь возникает при стохастической неоднородности данных в различных областях: экономике и финансах [33, 34], медицине [35], метеорологии [36], геофизике [37], технике [11, 38, 39, 40] и других.

2.5. Анализ динамики безработицы в Свердловской области в период пандемии COVID-19. Согласно данным Федеральной службы государственной статистики, в Свердловской области в период пандемии COVID-19 зафиксирован значительный рост уровня регистрируемой безработицы. 

В развитие исследования [34] нами проанализированы ежемесячные данные за 2022–2023 годы. На рис. 5 представлена временная динамика уровня безработицы в регионе, демонстрирующая более чем пятикратное увеличение показателя в кризисный период (апрель 2020 – ноябрь 2021 гг.).

Рис. 5. Динамика уровня безработицы в Свердловской области за 2016–2023 годы
[Figure 5. Unemployment rate dynamics in Sverdlovsk Region, 2016-2023]

Для количественного анализа динамики этого показателя построена линейная регрессионная модель вида 
\[
Y = a_1X_1 + a_2X_2,
\]
где $Y$ — уровень регистрируемой безработицы (%), $X_1 = 1$ (константа), $X_2$ — порядковая переменная (1–96, соответствуя месяцам с января 2016 г. по декабрь 2023 г.).

Результаты оценивания параметров модели представлены в табл. 3. Особый интерес представляют расчеты по данным, исключающим пандемийный период¹ (последняя строка таблицы), которые соответствуют предпосылкам $1^\circ$–$7^\circ$ и могут рассматриваться как эталонные для всего анализируемого периода. 

Таблица 3. Результаты оценивания параметров регрессионной модели [Regression model parameter estimation results]

Метод [Method]	$a_1$	$a_2$	$R^2$
ВМНМ [WLAD (Weighted Least Absolute Deviations)]	$1.517$	$-0.00897$	$0.182$
МНМ  [LAD (Least Absolute Deviations)]	$1.508$	$-0.00833$	$0.164$
ОМНМ  [GLAD (Generalized Least Absolute Deviations)]	$1.418$	$-0.00806$	$0.213$
ВМНК [WLS (Weighted Least Squares)]	$1.966$	$-0.00216$	$0.047$
МНК [OLS (Ordinary Least Squares)]	$1.669$	$-0.00177$	$0.002$
МНК (без пандемии) [OLS (excl. pandemic)]	$1.493$	$-0.00896$	$0.853$

Анализ результатов позволяет сделать следующие выводы:

• наибольшую устойчивость к выбросам при сохранении точности оценивания демонстрируют ВМНМ и ОМНМ;
• ВМНМ точнее оценивает коэффициенты тренда, но менее устойчив к односторонним выбросам ($R^2 = 0.182$ против $0.213$ у ОМНМ);
• традиционный МНК показал низкую эффективность ($R^2 = 0.002$) и не позволяет построить трендовую модель, хотя его взвешенная версия (ВМНК) дает некоторое улучшение ($R^2 = 0.047$);
• на данных без пандемийного периода ВМНМ и ОМНМ демонстрируют статистически значимые результаты ($R^2 = 0.847$ и $0.830$ соответственно).

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

Заключение

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

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

Когда приоритетом является максимальная точность без жестких ограничений на вычислительные затраты, особенно при комплексных нарушениях предпосылок регрессионного анализа, рекомендуется использовать ОМНМ, реализованный как алгоритм обобщенного спуска, избегая функций потерь с условием $\lim\limits_{x \to +0} {\rho(x)}/{x} = \infty$.

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

Конкурирующие интересы. Авторы заявляют об отсутствии конфликта интересов.
Авторский вклад. Все авторы внесли существенный вклад в исследование: участвовали в разработке методологии, проведении анализа и подготовке рукописи. Окончательная версия статьи была одобрена всеми соавторами перед подачей в печать.
Финансирование. Работа выполнена в соответствии с планом научно-исследовательских работ Института экономики Уральского отделения РАН.
Благодарности. Авторы выражают искреннюю признательность анонимным рецензентам за внимательное прочтение работы и конструктивные замечания, которые способствовали улучшению качества публикации.

---

¹Из выборки были удалены наблюдения за период пандемии.

About the authors

Oleg A. Golovanov

Institute of Economics, Ural Branch of RAS

Email: golovanov.oa@uiec.ru
ORCID iD: 0000-0002-9977-6954
SPIN-code: 4130-8355
Scopus Author ID: 58522704600
https://www.mathnet.ru/rus/person206252

Junior Researcher, Center for Economic Security¹

Russian Federation, 620014, Yekaterinburg, Moskovskaya str., 29

Alexander N. Tyrsin

Ural Federal University named after the first President of Russia B. N. Yeltsin; Institute of Economics, Ural Branch of RAS

Author for correspondence.
Email: at2001@yandex.ru
ORCID iD: 0000-0002-2660-1221
SPIN-code: 1408-1093
Scopus Author ID: 8503427500
ResearcherId: T-5975-2017
https://www.mathnet.ru/rus/person29355

Dr. Tech. Sci., Professor; Leading Researcher; Center for Economic Security¹; Head of Department; Dept. of Applied Mathematics and Mechanics²

Russian Federation, 620002, Yekaterinburg, Mira str., 19; 620014, Yekaterinburg, Moskovskaya str., 29

References

Supplementary files