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

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

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

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

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

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

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

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

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

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

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

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

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

Любая линейная функция может быть записана уравнением

{\displaystyle{y = a \cdot x + b}}

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

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

{\displaystyle{ F(a,b)= \sum_{i=1}^n (y_i — (a \cdot x_i + b))^2 \rarr min }}

Решение поставленной задачи сводится к нахождению экстремума указанной функции двух переменных. С этой целью находим частные производные функции функции по коэффициентам {a} и {b} и приравниваем их к нулю.

{\displaystyle{ \begin{cases} \displaystyle{\frac {\partial F(a,b)}{\partial a} = -2 \sum_{i=1}^n (y_i — (a \cdot x_i + b)) \cdot x_i = 0} \\ \\ \displaystyle{\frac {\partial F(a,b)}{\partial b} = -2 \sum_{i=1}^n (y_i — (a \cdot x_i + b)) = 0} \\ \end{cases} }}

Решаем полученную систему уравнений

{\displaystyle{ \begin{cases} \displaystyle{a \cdot \sum_{i=1}^n x_i^2 + b \cdot \sum_{i=1}^n x_i = \sum_{i=1}^n x_i \cdot y_i} \\ \\ \displaystyle{a \cdot \sum_{i=1}^n x_i + \sum_{i=1}^n b = \sum_{i=1}^n y_i}\\ \end{cases} }}

{\displaystyle{ \begin{cases} \displaystyle{a \cdot \sum_{i=1}^n x_i^2 + b \cdot \sum_{i=1}^n x_i = \sum_{i=1}^n x_i \cdot y_i} \\ \\ \displaystyle{a \cdot \sum_{i=1}^n x_i + n \cdot b = \sum_{i=1}^n y_i}\\ \end{cases} }}

Определяем значения коэффициентов

{\displaystyle{ \begin{cases} a = {\displaystyle{ \frac {n \cdot \displaystyle{\sum_{i=1}^n x_i \cdot y_i}\ -\ \displaystyle{\sum_{i=1}^n x_i \cdot \sum_{i=1}^n y_i} } {n \cdot \displaystyle{\sum_{i=1}^n x_i^2}\ -\ \displaystyle{ \left( \sum_{i=1}^n x_i \right)^2}} }} \\ \\ b = {\displaystyle{ \frac {\displaystyle{\sum_{i=1}^n y_i}\ -\ a \cdot \displaystyle{\sum_{i=1}^n x_i}} {n} }}\\ \end{cases} }}

Для вычисления коэффициентов необходимо найти следующие составляющие:

{\displaystyle{ sumx = \displaystyle{\sum_{i=1}^n x_i} }}

{\displaystyle{ sumy = \displaystyle{\sum_{i=1}^n y_i} }}

{\displaystyle{ sumx2 = \displaystyle{\sum_{i=1}^n x_i^2} }}

{\displaystyle{ sumxy = \displaystyle{\sum_{i=1}^n x_i \cdot y_i} }}

Тогда значения коэффициентов будут определены как

{\displaystyle{ \begin{cases} a = {\displaystyle{ \frac {n \cdot sumxy \ -\ sumx \cdot sumy} {n \cdot sumx2 \ — \ \displaystyle{ \left( sumx \right)^2}} }} \\ \\ b = {\displaystyle{ \frac {sumy \ -\ a \cdot sumx} {n}} }\\ \end{cases} }}

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

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

y = 8 \cdot x — 3

Функция getData() служит для формирования двух массивов для координат x и y согласно уравнению для указанного количества точек.

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

Функция getApprox() рассчитывает коэффициенты аппроксимирующей прямой. И поскольку функция не может вернуть сразу 2 значения коэффициентов, используется передача «местодержателей» (указателей) для a и b.

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

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
#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
// Задание начального набора значений
void getData(int n, double *x, double *y) 
{
    for (int i = 0; i < n; i++) 
    {
        x[i] = (double)i;
        y[i] = 8 * (double)i - 3;
        // Добавление случайной составляющей
        //y[i] = 8 * (double)i - 3 + ((rand() % 100) - 50) * 0.05;
    }
    return;
}
// Вычисление коэффициентов аппроксимирующей прямой
void getApprox(double* x, double *y, int n, double* a, double* b) 
{
    double sumx = 0;
    double sumy = 0;
    double sumx2 = 0;
    double sumxy = 0;
    for (int i = 0; i < n; i++) 
    {
        sumx += x[i];
        sumy += y[i];
        sumx2 += x[i] * x[i];
        sumxy += x[i] * y[i];
    }
    *a = (n * sumxy - (sumx * sumy)) / (n * sumx2 - sumx * sumx);
    *b = (sumy - *a * sumx) / n;
    return;
}
int main() 
{
    double* x, *y, a, b;
    int n;
    system("chcp 1251");
    system("cls");
    printf("Введите количество точек: ");
    scanf("%d", &n);
    x = new double[n];
    y = new double[n];
    getData(n, x, y);
    for (int i = 0; i < n; i++)
        printf("%5.1lf - %7.3lf\n", x[i], y[i]);
    getApprox(x, y, n, &a, &b);
    printf("a = %lf\nb = %lf", a, b);
    delete[] x;
    delete[] y;
    getchar(); getchar();
    return 0;
}

Результат выполнения: Запуск без случайной составляющей

Запуск со случайной составляющей

Реализация линейной аппроксимации по МНК

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

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

{\displaystyle{ \displaystyle{\frac {\partial F(a,b)}{\partial a} = -2 \sum_{i=1}^n (y_i — (a \cdot x_i + b)) \cdot x_i = 0} }}

{\displaystyle{ \displaystyle{a \cdot \sum_{i=1}^n x_i^2 + b \cdot \sum_{i=1}^n x_i = \sum_{i=1}^n x_i \cdot y_i} }}

{\displaystyle{ a = {\displaystyle{ \frac {\displaystyle{\sum_{i=1}^n x_i \cdot y_i}\ -\ b \cdot \displaystyle{\sum_{i=1}^n x_i} } {\displaystyle{\sum_{i=1}^n x_i^2}\ } }} }}

В этом случае текст программы для поиска коэффициента угла наклона аппроксимирующей прямой будет следующий (имя функции 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
47
48
49
50
51
#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
// Задание начального набора значений
void getData(int n, double *x, double *y) 
{
    for (int i = 0; i < n; i++) 
    {
        x[i] = (double)i;
        y[i] = 8 * (double)i - 3;
        // Добавление случайной составляющей
        //y[i] = 8 * (double)i - 3 + ((rand() % 100) - 50) * 0.05;
    }
    return;
}
// Вычисление коэффициентов аппроксимирующей прямой
void getApproxA(double* x, double *y, int n, double* a, double *b) 
{
    double sumx = 0;
    double sumx2 = 0;
    double sumxy = 0;
    for (int i = 0; i < n; i++) 
    {
        sumx += x[i];
        sumx2 += x[i] * x[i];
        sumxy += x[i] * y[i];
    }
    *a = (sumxy - *b * sumx) / sumx2;
    return;
}
int main() 
{
    double* x, * y;
    double a, b = 0; // заданная точка пересечения с осью y
    int n;
    system("chcp 1251");
    system("cls");
    printf("Введите количество точек: ");
    scanf("%d", &n);
    x = new double[n];
    y = new double[n];
    getData(n, x, y);
    for (int i = 0; i < n; i++)
        printf("%5.1lf - %7.3lf\n", x[i], y[i]);
    getApproxA(x, y, n, &a, &b);
    printf("a = %lf\nb = %lf", a, b);
    delete[] x;
    delete[] y;
    getchar(); getchar();
    return 0;
}

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

Прокрутить вверх