Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория

в 6:16, , рубрики: C#, аппроксимация, регрессия

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

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

Небольшая теория  аппроксимации функции  методом наименьших квадратов

Есть выборка данных зависимости значений X_i от Y_i

X_i

Y_i

0

0

1

1

2

3

5

10

Мы знаем что, например,  при X равным 1, Y будет равен 1, но нам нужно узнать значения между узловыми точками или узнать значение Y при X равным 10. В таком случае используют регрессионный анализ. Один из вариантов получения функции по заданным точкам это аппроксимация. В конечном итоге должна быть построена функция Y^*  значения которой при X будут примерно равны значениям Y.

X_i

Y_i

Y_i-Y^*(X_i)

0

0

≈0

1

1

≈0

2

3

≈0

5

10

≈0

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

displaystylesum_{i=0}^{n} (Y_i-Y^*(x_i))^2rightarrow0

Т.е. в результате должны быть найдены коэффициенты функции Y^*(x_i), чтобы разница между точками Yi и значениями функции Y^*(x_i).

Для функции аппроксимирования я использую числовой ряд:

displaystylesum_{i=0}^{n} a_nx^n

На моей практике максимальную повторяемость реальной функции даст n=4. Для простоты следующего примера я буду использовать функцию  n = 1 (линейную функцию kx+b)

Возьмем линейную функцию Y^*(x)=a_1x+a_0, тогда аппроксимируемая функция будет выглядеть

F(x)=sum_{i=0}^n=(Y_i-a_1x_i+a_0)^2

Далее нужно найти коэффициенты a1, a2. Для этого составим систему уравнений из частных производных по данной функции.

Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 23

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

Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 24

Выразим из второго уравнения a_0 :

a_0=frac{(sum_{i=0}^ny_i-a_1sum_{i=0}^nx_i)}{n}

Тогда

a_1=frac{nsum_{i=0}^nx_iy_i-sum_{i=0}^nx_isum_{i=0}^ny_i}{nsum_{i=0}^nx_i^2-(sum_{i=0}^nx_i)^2}

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

Подставим значения из выборки и найдем коэффициенты:

sum_{i=0}^nx_iy_i=57 sum_{i=0}^nx_isum_{i=0}^ny_i=112nsum_{i=0}^nx_i^2=120(sum_{i=0}^nx_i)^2=64 sum_{i=0}^ny_i=14sum_{i=0}^nx_i=8a_1=frac{4*57-112}{120-64}=2.07a_0=frac{14-2.07*8}{4}=-0.64

В конечном итоге мы построили регрессионную функцию, которая проходит через наши точки с минимально возможным отклонением. Обычно линейную функцию используют для определения тренда выборки, т.е. мы знаем, что с ростом значений X, Y будет возрастать в 2.07 раза. Для более “точного повторения” используют функции 4 или 5 степени.

X_i

Y_i

Y_i-Y*(X_i)

0

0

0

1

1

1.43

2

3

3.5

5

10

9.71

Реализация на платформе .NET

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

Для этого понадобится библиотека MathNet.Symbolics, которую можно установить через диспетчера пакетов NUuGet . Данная библиотека будет вычислять частные производные, а так же упрощать выражения.

Вычислять систему уравнений будем с помощью метода Гаусса https://en.wikipedia.org/wiki/Gaussian_elimination , реализацию можно найти в интернете или в моей библиотеке Regression в папке MathService.

Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 41
static void Main()
     {
         //Символьное выражение членов полинома a1 * x + a2 
         var a1 = SymbolicExpression.Variable("a1");
         var a2 = SymbolicExpression.Variable("a2");

         //Аппроксимируемая выборка значений пара X Y
         Dictionary<double, double> XtoY = new()
         {
             { 0, 0 }, { 1, 1 }, { 2, 3 }, { 5, 10 }
         };

         //Построение функции F(x) 
         SymbolicExpression Fx = 0;
         foreach (var item in XtoY)
         {
             var expr = item.Value + (a1 * item.Key + a2);
             Fx += expr * expr;
         }

         //Частные производные в строковом виде
         var diffA1 = Fx.Differentiate(a1).RationalSimplify(a1).ToString().Trim().Split();
         var diffA2 = Fx.Differentiate(a2).RationalSimplify(a2).ToString().Trim().Split();

         List<string[]> polysList = new() { diffA1, diffA2 };
         // Step 1 Построение матрицы Y
         double[,] Y = new double[2, 1];
         for (int i = 0; i < Y.GetUpperBound(0) + 1; i++)
         {
             var yValues = GetArrayParseY(polysList[i]).ToArray();
             for (int j = 0; j < Y.GetUpperBound(1) + 1; j++)
             {
                 Y[i, 0] += yValues[j];
             }
         }

         // Step 2 Построение матрицы X
         double[,] X = new double[2, 2];
         for (int i = 0; i < X.GetUpperBound(0) + 1; i++)
         {
             var xValues = GetArrayParseX(polysList[i]).ToArray();
             for (int j = 0; j < X.GetUpperBound(1) + 1; j++)
             {
                 X[i, j] += xValues[j];
             }
         }
         // Step 3 Решение системы уравнений
         Gaus gaus = new();
         double[] result = gaus.Roots(X, Y);
     }

/// <summary>
 /// Вернуть из строк массив X значений
 /// </summary>
 /// <param name="strings"></param>
 /// <returns></returns>
 private static IEnumerable<double> GetArrayParseX(string[] strings)
 {
     var a = strings.Where(x => !x.Equals("+"));
     var listAr = new double[2];
     int sign = 1;
     foreach (var stringExpression in a)
     {
         if (stringExpression.Equals("-"))
         {
             sign = -1;
             continue;
         }

         if (!stringExpression.Contains("*a")) continue;
         var str = stringExpression.Remove(0, stringExpression!.IndexOf("a", StringComparison.Ordinal)+1);
         var index = int.Parse(str);
         var st1r = stringExpression.Substring(0, stringExpression!.IndexOf("*", StringComparison.Ordinal));
         listAr[index-1] = sign * double.Parse(st1r, CultureInfo.InvariantCulture);
         sign = 1;
     }
     return listAr;
 }
 /// <summary>
 /// Вернуть из строки массив Y значений
 /// </summary>
 /// <param name="strings"></param>
 /// <returns></returns>
 private static IEnumerable<double> GetArrayParseY(string[] strings)
 {
     var a = strings.Where(x => !x.Equals("+"));
     var sign = 1;
     foreach (var stringExpression in a)
     {
         if (stringExpression.Equals("-"))
         {
             sign = -1;
             continue;
         }

         if (stringExpression.Contains("*a")) continue;
         yield return sign * double.Parse(stringExpression, CultureInfo.InvariantCulture);
         sign = 1;
     }
 }

В результате выполнения кода

SymbolicExpression Fx = 0;
         foreach (var item in XtoY)
         {
             var expr = item.Value + (a1 * item.Key + a2);
             Fx += expr * expr;
         }

функция F(x) будет равна
(0.0 + a_0)^2 + (1.0 + 1.0*a_1 + a_0)^2 + (3.0 + 2.0*a_1 + a_0)^2 + (10.0 + 5.0*a_1 + a_0)^2

Далее частные производные функции по а_1 - frac{dF(x)}{da_1} и а_0 - frac{dF(x)}{da_0} вернут значения

frac{dF(x)}{da_1}=114.0 + 60.0*a_1 + 16.0*a_0

frac{dF(x)}{da_0}=28.0 + 16.0*a_1 + 8*a_0

Остается составить матрицу для решения системы уравнений методом Гаусса вида A*X=B

Матрица A

A=begin{matrix}60&16\16&8end{matrix}

Столбец свободных коэффициентов B

B=begin{matrix}114\28end{matrix}

Столбец переменных

X=begin{matrix}a_0\a_1end{matrix}

Для составления матриц из массива строк были применены методы GetArrayParseX и GetArrayParseY

В результате получаем массив коэффициентов полинома:

result[0] = 2.07

result[1] = -0.64

Таким образом можно реализовывать многофакторную регрессию, когда на входе больше, чем одно значение X. Сделать это можно перемножением полиномов друг на друга:

Возьмем полином второй степени вида:

a_0+a_1+a_2x^2+a_3x^3

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

(a_0+a_1x_1+a_2x_1^2+a_3x^3)(a'_0+a'_1x_2+a'_2x^2_2+a'_3x^3_2)=a_0+a_1*x_2+a_2*x_2^2+a_3**x_1+a_4*x_1^2+a_5*x_1*x_2+a_6*x_1^2*x_2+a_7*x_1*x_2^2+a_8*x_1^2*x_2^2+a_9*x_1^3++a_{10}*x_1^3*x_2+a_{11}*x_1^3*x_2^2+a_{12}*x_2^3+a_{13}*x_1*x_2^3+a_{14}*x_2^3*x_1^2+a_{15}*x_2^3**x_1^3

Практическое применение

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

Давление на эталоне:

Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 59

Код давления и температуры при 16 градусах:

код давления при 16 градусах

код давления при 16 градусах
код температуры при 16 градусах

код температуры при 16 градусах

Код давления и температуры при 18 градусах

код давления при 18 градусах

код давления при 18 градусах
код температуры при 18 градусах

код температуры при 18 градусах

Код давления и температуры при 21 градусе

код давления при 21 градусе

код давления при 21 градусе
код температуры при 21 градусе

код температуры при 21 градусе

Код давления и температуры при 23 градусах

код давления при 23 градусах

код давления при 23 градусах
код температуры при 23 градусах

код температуры при 23 градусах

Код давления и температуры при 27 градусах

код давления при 27 градусах

код давления при 27 градусах
код температуры при 27 градусах

код температуры при 27 градусах

Объединим все наши коллекции:

Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 70
Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 71
Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 72

Функция аппроксимирования полиномом третьей степени будет выглядеть:

Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 73
Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 74

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

Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 75
Многофакторное аппроксимирование на платформе .Net. Часть первая. Теория - 76

Приложение RegressionFromExcel

Пользоваться MathCad было не очень удобно:

·        Время расчета коэффициентов могло достигать до 1 минуты

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

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

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

RegressionFromExcel/README.md at main · Georgiy-smr/RegressionFromExcel

В папке Exemple лежит файл TestData, в котором уже есть тестовые данные, а так же посчитанные коэффициенты. Для расчета нужно заполнить столбцы с данными и  выбрать степень полинома в приложении 2 или 3. Далее приложение откроет файл, считает данные и вычислит коэффициенты, которые вставит в отдельные ячейки справа (как в примере TestData).

Выбор степени полинома, по которому будет рассчитаны коэффициенты двухфакторной аппроксимации

Выбор степени полинома, по которому будет рассчитаны коэффициенты двухфакторной аппроксимации

Более подробно о работе приложения будет описано во второй части.

Сравним результаты вычислений MathCad и RegressionFromExcel.

С коэффициентами MathCad

С коэффициентами MathCad
Приложение RegressionFromExcel

Приложение RegressionFromExcel

Мы видим достаточную сходимость погрешностей на уровне 2% от допустимой погрешности на диапазон преобразователя давления.

Использованная литература:

  1. Кобзарь А. И. Прикладная математическая статистика. — М.: Физматлит, 2006.

  2. Методы корреляционного и регрессионного анализа

Автор: mrsmirnovgd

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js