Линейная аппроксимация

Алгоритмизация / Линейная аппроксимация

 

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

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

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

В случае если аппроксимация проводится на непрерывном множестве точек (отрезке), аппроксимация называется непрерывной или интегральной. Примером такой аппроксимации может служить разложение функции в ряд Тейлора, то есть замена некоторой функции степенным многочленом.

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

Пусть задан дискретный набор точек, называемых узлами интерполяции, а также значения функции в этих точках. Требуется построить функцию g(x), проходящую наиболее близко ко всем заданным узлам. Таким образом, критерием близости функции является g(xi)=yi.

В качестве функции g(x) обычно выбирается полином, который называют интерполяционным полиномом.

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

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

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

Аппроксимация линейной функцией

Любая линейная функция может быть записана уравнением
Уравнение прямой
Аппроксимация заключается в отыскании коэффициентов a и b уравнения таких, чтобы все экспериментальные точки лежали наиболее близко к аппроксимирующей прямой.

С этой целью чаще всего используется метод наименьших квадратов (МНК), суть которого заключается в следующем: сумма квадратов отклонений значения точки от аппроксимирующей точки принимает минимальное значение:
Метод наименьших квадратов
Решение поставленной задачи сводится к нахождению экстремума указанной функции двух переменных. С этой целью находим частные производные функции функции по коэффициентам a и b и приравниваем их к нулю.
Частные производные МНК
Решаем полученную систему уравнений
Частные производные МНК
Определяем значения коэффициентов
Коэффициенты линейной аппроксимации по МНК
Для вычисления коэффициентов необходимо найти следующие составляющие:
МНК
Тогда значения коэффициентов будут определены как
Коэффициенты линейной аппроксимации

Пример реализации

Для примера реализации воспользуемся набором значений, полученных в соответствии с уравнением прямой

y = 8 · x — 3

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

Реализация на Си

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
// Задание начального набора значений
double ** getData(int n) {
  double **f;
  f = new double*[2];
  f[0] = new double[n];
  f[1] = new double[n];
  for (int i = 0; i<n; i++) {
    f[0][i] = (double)i;
    f[1][i] = 8 * (double)i - 3;
    // Добавление случайной составляющей
    f[1][i] = 8*(double)i - 3 + ((rand()%100)-50)*0.05;
  }
  return f;
}
// Вычисление коэффициентов аппроксимирующей прямой
void getApprox(double **x, double *a, double *b, int n) {
  double sumx = 0;
  double sumy = 0;
  double sumx2 = 0;
  double sumxy = 0;
  for (int i = 0; i<n; i++) {
    sumx += x[0][i];
    sumy += x[1][i];
    sumx2 += x[0][i] * x[0][i];
    sumxy += x[0][i] * x[1][i];
  }
  *a = (n*sumxy - (sumx*sumy)) / (n*sumx2 - sumx*sumx);
  *b = (sumy - *a*sumx) / n;
  return;
}
int main() {
  double **x, a, b;
  int n;
  system("chcp 1251");
  system("cls");
  printf("Введите количество точек: ");
  scanf("%d", &n);
  x = getData(n);
  for (int i = 0; i<n; i++)
    printf("%5.1lf - %7.3lf\n", x[0][i], x[1][i]);
  getApprox(x, &a, &b, n);
  printf("a = %lf\nb = %lf", a, b);
  getchar(); getchar();
  return 0;
}

Результат выполнения
Запуск без случайной составляющей
Реализация линейной аппроксимации по МНК
Запуск со случайной составляющей
Реализация линейной аппроксимации по МНК

Построение графика функции

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

 
Реализация на Си

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#include <windows.h>
const int NUM = 70; // количество точек
LONG WINAPI WndProc(HWNDUINTWPARAMLPARAM);
double **x; // массив данных
      // Определение коэффициентов линейной аппроксимации по МНК
void getApprox(double **m, double *a, double *b, int n) {
  double sumx = 0;
  double sumy = 0;
  double sumx2 = 0;
  double sumxy = 0;
  for (int i = 0; i<n; i++) {
    sumx += m[0][i];
    sumy += m[1][i];
    sumx2 += m[0][i] * m[0][i];
    sumxy += m[0][i] * m[1][i];
  }
  *a = (n*sumxy - (sumx*sumy)) / (n*sumx2 - sumx*sumx);
  *b = (sumy - *a*sumx) / n;
  return;
}
// Задание исходных данных для графика
// (двумерный массив, может содержать несколько рядов данных)
double ** getData(int n) {
  double **f;
  double a, b;
  f = new double*[3];
  f[0] = new double[n];
  f[1] = new double[n];
  f[2] = new double[n];
  for (int i = 0; i<n; i++) {
    double x = (double)i * 0.1;
    f[0][i] = x;
    f[1][i] = 8 * x - 3 + ((rand() % 100) - 50)*0.05;
  }
  getApprox(f, &a, &b, n); // аппроксимация
  for (int i = 0; i<n; i++) {
    double x = (double)i * 0.1;
    f[2][i] = a*x + b;
  }
  return f;
}
// Функция рисования графика
void DrawGraph(HDC hdc, RECT rectClient, double **x, int n, int numrow = 1) {
  double OffsetY, OffsetX;
  double MAX_X = 0;
  double MAX_Y = 0;
  double ScaleX, ScaleY;
  double min, max;
  int height, width;
  int X, Y;    // координаты в окне (в px)
  HPEN hpen;
  height = rectClient.bottom - rectClient.top;
  width = rectClient.right - rectClient.left;
  // Область допустимых значений X
  min = x[0][0];
  max = x[0][0];
  for (int i = 0; i<n; i++) {
    if (x[0][i] < min)
      min = x[0][i];
    if (x[0][i] > max)
      max = x[0][i];
  }
  double temp = max - min;
  MAX_X = max - min;
  OffsetX = min*width / MAX_X;  // смещение X
  ScaleX = (double)width / MAX_X; // масштабный коэффициент X
                  // Область допустимых значений Y
  min = x[1][0];
  max = x[1][0];
  for (int i = 0; i<n; i++) {
    for (int j = 1; j <= numrow; j++) {
      if (x[j][i] < min)
        min = x[j][i];
      if (x[j][i] > max)
        max = x[j][i];
    }
  }
  MAX_Y = max - min;
  OffsetY = max*height / (MAX_Y);  // смещение Y
  ScaleY = (double)height / MAX_Y; // масштабный коэффициент Y
                   // Отрисовка осей координат
  hpen = CreatePen(PS_SOLID, 0, 0); // черное перо 1px
  SelectObject(hdc, hpen);
  MoveToEx(hdc, 0, OffsetY, 0);  // перемещение в точку (0;OffsetY)
  LineTo(hdc, width, OffsetY);   // рисование горизонтальной оси
  MoveToEx(hdc, OffsetX, 0, 0);  // перемещение в точку (OffsetX;0)
  LineTo(hdc, OffsetX, height);  // рисование вертикальной оси
  DeleteObject(hpen);  // удаление черного пера
             // Отрисовка графика функции
  int color = 0xFF; // красное перо для первого ряда данных
  for (int j = 1; j <= numrow; j++) {
    hpen = CreatePen(PS_SOLID, 2, color); // формирование пера 2px
    SelectObject(hdc, hpen);
    X = (int)(OffsetX + x[0][0] * ScaleX);    // координаты начальной точки графика
    Y = (int)(OffsetY - x[j][0] * ScaleY);
    MoveToEx(hdc, X, Y, 0);  // перемещение в начальную точку
    for (int i = 0; i<n; i++) {
      X = OffsetX + x[0][i] * ScaleX;
      Y = OffsetY - x[j][i] * ScaleY;
      LineTo(hdc, X, Y);
    }
    color = color << 8; // изменение цвета пера для следующего ряда
    DeleteObject(hpen); // удаление текущего пера
  }
}
// Главная функция
int  WINAPI  WinMain(HINSTANCE  hInstance,
  HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow) {
  HWND hwnd;
  MSG msg;
  WNDCLASS w;
  x = getData(NUM); // задание исходных данных
  memset(&w, 0, sizeof(WNDCLASS));
  w.style = CS_HREDRAW | CS_VREDRAW;
  w.lpfnWndProc = WndProc;
  w.hInstance = hInstance;
  w.hbrBackground = CreateSolidBrush(0x00FFFFFF);
  w.lpszClassName = "My Class";
  RegisterClass(&w);
  hwnd = CreateWindow("My Class", "График функции",
    WS_OVERLAPPEDWINDOW, 500, 300, 500, 380, NULLNULL,
    hInstance, NULL);
  ShowWindow(hwnd, nCmdShow);
  UpdateWindow(hwnd);
  while (GetMessage(&msg, NULL, 0, 0)) {
    TranslateMessage(&msg);
    DispatchMessage(&msg);
  }
  return msg.wParam;
}
// Оконная функция
LONG WINAPI WndProc(HWND hwnd, UINT Message,
  WPARAM wparam, LPARAM lparam) {
  HDC hdc;
  PAINTSTRUCT ps;
  switch (Message) {
  case WM_PAINT:
    hdc = BeginPaint(hwnd, &ps);
    DrawGraph(hdc, ps.rcPaint, x, NUM, 2); // построение графика
    EndPaint(hwnd, &ps);
    break;
  case WM_DESTROY:
    PostQuitMessage(0);
    break;
  default:
    return DefWindowProc(hwnd, Message, wparam, lparam);
  }
  return 0;
}

Результат выполнения
Реализация линейной аппроксимации по МНК (график)

Аппроксимация с фиксированной точкой пересечения с осью y

В случае если в задаче заранее известна точка пересечения искомой прямой с осью y, в решении задачи останется только одна частная производная для вычисления коэффициента a.
Частная производная по a
В этом случае текст программы для поиска коэффициента угла наклона аппроксимирующей прямой будет следующий (имя функции getApprox() заменено на getApproxA() во избежание путаницы).
 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
// Задание начального набора значений
double ** getData(int n) {
  double **f;
  f = new double*[2];
  f[0] = new double[n];
  f[1] = new double[n];
  for (int i = 0; i<n; i++) {
    f[0][i] = (double)i;
    f[1][i] = 8 * (double)i - 3;
    // Добавление случайной составляющей
    //f[1][i] = 8 * (double)i - 3 + ((rand() % 100) - 50)*0.05;
  }
  return f;
}
// Вычисление коэффициентов аппроксимирующей прямой
void getApproxA(double **x, double *a, double b, int n) {
  double sumx = 0;
  double sumx2 = 0;
  double sumxy = 0;
  for (int i = 0; i<n; i++) {
    sumx += x[0][i];
    sumx2 += x[0][i] * x[0][i];
    sumxy += x[0][i] * x[1][i];
  }
  *a = (sumxy - b*sumx) / sumx2;
  return;
}
int main() {
  double **x, a, b;
  int n;
  system("chcp 1251");
  system("cls");
  printf("Введите количество точек: ");
  scanf("%d", &n);
  x = getData(n);
  for (int i = 0; i<n; i++)
    printf("%5.1lf - %7.3lf\n", x[0][i], x[1][i]);
  b = 0;
  getApproxA(x, &a, b, n);
  printf("a = %lf\nb = %lf", a, b);
  getchar(); getchar();
  return 0;
}

Результат выполнения программы поиска коэффициента угла наклона аппроксимирующей прямой при фиксированном значении b=0:
Аппроксимация при фиксированном b


Назад: Алгоритмизация

Комментариев к записи: 23

  • Код не совсем рабочий. В случае (n*sumx2 - sumx*sumx)==0.0; Получим в результате a=-inf,b=inf Это дают точки расположенные вдоль вертикальной линии. Как в этом случае д.б. доработан код?

  • Роман Алексеевич
    Спасибо, вы сэкономили мне кучу времени! Правда, я тут же потратил ещё больше, увлекшись вашим сайтом )). Итог - за то же время два результата: готовый велосипед + полезная закладка + положительные эмоции. Мне всегда нравились женщины с логическим-аналитическим складом мышления. Успехов вам!

  • Дмитрий
    А как можно найти погрешность найденного данным способом коэффициента наклона (а)?

    • Елена Вставская
      Если мы знаем точное значение коэффициента угла наклона Ареал., то относительную погрешность можно посчитать как |Аизм.-Ареал.|/Ареал.*100%


  • Елена, спасибо за статью! Если позволите, "недокомплимент" взяла себе в карман, теперь знаю как отвечать безынтересным людям со зрением -10!))))))))

  • Отличный пример и работа, а если известна начальная точка пересечения с Y. то коэффициент a изменится. Как можно это встроить в решение задачи?

  • Ахмеров Роман
    В задачах интерполяции (и линейной и кубической) использовал определитель Ванденмонда как универсальное средство не вдаваясь в суть. Теперь примерно понял. Спасибо. P.S У Петра свои нереализованные подростковые фантазии видимо остались. )))

  • Николай
    День добрый (аль вечер), В интернете нашёл вариант МНК, в котором минимизируется не разница по координате y, а расстояние от точки до прямой. Там делается некое допущение a*a+b*b=1 (дабы конкретизировать коэффициенты a, b, c) и после этого всё сводится к формуле E(a*xi+b*yi+c)^2=>min - таким образом становится очень сильно похоже на то, что делаете Вы. Однако, на этом решение прекращается. Я попробовал Вашим способом - найти производные от данной функции по трём коэффициентам. Получилось три уравнения, но эти уравнения без свободных членов, таким образом их единственным решением получается a=0, b=0, c=0. Вот эти уравнения: a*sumx2+b*sumxy+c*sx=0 a*sumxy+b*sumy2+c*sumy=0 a*sumx+b*sumy+c*n=0 подозреваю, что надо каким-то образом использовать ограничение a*a+b*b=1, но тогда там квадратный корень получается, что окончательно ставит меня в тупик.

    • Елена Вставская
      Что касается расстояния от точки до прямой, то я считаю, что это - дополнительное усложнение алгоритма. Для метода наименьших квадратов и аппроксимации линейной функцией (как он описан в статье) единственным исключительным вариантом является прямая x=a. Но если применять этот метод к реальному набору исходных данных, то для значений x и y всё равно будет какой-либо разброс, и этот алгоритм можно использовать с некоторым допущением.

  • Николай
    Пытаюсь использовать данный подход для распознавания линий на картинке. Не очень хорошо работает на вертикальных линиях (x=5 и т. п.) - из-за равенства уменьшаемого и вычитаемого получается неопределённость 0/0. Не подскажете ли, нельзя ли данный способ как-то переделать, чтобы он работал для уравнения прямой вида A x + B y + C = 0? Ещё есть беда, что при большом количестве точек очень легко выйти за пределы значений 32-битных значений. Пока лечу элементарным прореживанием данных (есть мысль, что можно усреднять группы точек, но пока гложат сомнения). Не силён в математике, поэтому осмелюсь спросить у автора: нельзя ли данный алгоритм переделать, чтобы он работал по добавлению каждой новой точки? Не помню, как эта штука называется, но есть функции-агрегаты, которые так можно легко реализовать (например, sum - каждое новое значение просто прибавляем к результату), есть некоторые, которые можно, но с приворотами (среднее-арифметическое по-моему), а какие-то вообще нельзя (медиану, вроде как). Данный случай к какому относится?

    • Елена Вставская
      Здравствуйте, Николай! Да, для вертикальных прямых этот алгоритм применять затруднительно. Если Вы хотите использовать его для уравнения прямой вида A x + B y + C = 0, то нужно просто составить функционал для этого случая, найти частные производные и вывести значения для коэффициентов. Что касается увеличения количества точек, то мы по сути и храним суммы значений (в различных вариантах), поэтому увеличение количества точек реализовать достаточно просто. Для этого нужно сделать суммы внутри функции переменными типа static (чтобы они сохраняли значение при выходе из функции), а затем с каждой новой точкой обращаться к функции для пересчета.

    • Может стоит просто привести уравнение вида A x + B y + C = 0 к y = B'x +C' ? B' = -A/B C' = -C/B ?



  • Какая-то вы не особо привлекательная, Елена. Неудачно на фото получились или действительно природа красотой обделила?




Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *