Aproksymacja

W matematyce termin aproksymacja (ang. approximation) oznacza coś przybliżonego, co zwykle posiada prostszą postać i rozwiązuje pewien problem z pewnym przybliżeniem, ale za to prościej. Interpolacja natomiast oznacza szukanie wartości pośrednich, pomiędzy wartościami dokładnie znanymi (tzw. węzłami interpolacji). W interpolacji funkcja interpolująca przechodzi dokładnie przez węzły interpolacji. W aproksymacji staramy się dopasować funkcję aproksymującą tak, aby przechodziła możliwie najbliżej węzłów aproksymacji - ale nie musi przechodzić przez nie dokładnie. Np. mamy trzy punkty i chcemy dobrać taką prostą, aby przechodziła jak najbliżej tych punktów:

Innym przykładem jest sytuacja, gdy mamy serię pomiarów i chcemy dopasować do nich pewną znaną funkcję tak, aby błędy tego dopasowania były jak najmniejsze.

Liniowa aproksymacja średniokwadratowa (ang. least-squares linear approximation) polega na dopasowaniu funkcji liniowej \({f(x) = ax + b}\) do zbioru  \(n\) punktów pomiarowych \((x_0,y_0)\), \((x_1,y_1)\), ..., \({(x_{n - 1},y_{n - 1})}\) tak, aby suma kwadratów odchyleń punktów od wartości tej funkcji liniowej była jak najmniejsza. Tworzymy funkcję kryterialną:

W celu znalezienia minimum funkcji \(S(a,b)\) wyznaczamy jej pochodne cząstkowe względem \(a\) i \(b\) i przyrównujemy je do zera:

Pochodne cząstkowe obliczamy:

Dzielimy oba równania obustronnie przez \((-2)\):

Wyrażenia wewnątrz sum wymnażamy:

Rozbijamy na pojedyncze sumy i wyciągamy stałe \(a\) oraz \(b\) przed znak sumowania:

Przenosimy niewiadome na lewą stronę, a wyrazy z \(y_i\) na prawą stronę, tworząc ostateczny układ równań normalnych dla prostej:

Układ rozwiązujemy względem \(a\) i \(b\), otrzymując:

We wzorach występują sumy o tych samych zakresach. Zastąpmy je symbolami:

Przed obliczeniem współczynnika \(a\) należy sprawdzić, czy mianownik ułamka jest równy zero. W takim przypadku obliczenia muszą być przerwane, bo oznacza to, iż wszystkie punkty danych wejściowych leżą na prostej pionowej, a takiej nie da się zapisać równaniem \({y = ax + b}\).

Otrzymujemy współczynniki prostej, która aproksymuje dane wejściowe linią prostą. Prosta ta pozwala wyznaczać punkty leżące poza zbiorem danych, co umożliwia nam aproksymacja. Przy interpolacji mogliśmy wyznaczać jedynie punkty pomiędzy danymi (słowo "inter" po łacinie znaczy "między", "pomiędzy"). Musisz jednak pamiętać, iż wyznaczona prosta aproksymacyjna jedynie przybliża prostą bazową, a błąd przybliżenia rośnie wraz z oddalaniem się od zbioru danych wejściowych.

Przykładowa implementacja

W celu przetestowania algorytmu aproksymacji średniokwadratowej poniższy program wykonuje następujące operacje:

• Określa:
  • \(N\) - liczbę punktów (np. 100)
  • przedział \([Xmin; Xmax]\)
  • współczynniki \(aa\) i \(bb\) prostej bazowej \({y = aax + bb}\)
• Wyznacza losowo w tym przedziale \(N\) punktów \(x\) i dla każdego z nich liczy \({y = aax + bb}\). Do wyniku dodaje losową wartość od -1 do 1, co tworzy odchyłkę i daje w wyniku punkty leżące ponad i pod prostą bazową.
• Wg wzorów aproksymacji średniokwadratowej wylicza współczynniki \(a\) i \(b\) dla prostej aproksymującej.
• Wyświetla wynik porównawczy: prostą bazową i prostą aproksymującą.

Wszystko jest opisane w komentarzach.

// Liniowa aproksymacja
// średniokwadratowa
// (C)2026 mgr Jerzy Wałaszek
// Metody numeryczne 0078
//---------------------------

#include <iostream>
#include <windows.h>
#include <iomanip>
#include <cmath>
#include <cstdlib>
#include <ctime>

using namespace std;

// Definicje stałych
//------------------

// Liczba danych
const int N = 100;
// Przedział dla danych
const double Xmin = -10;
const double Xmax =  10;
// Współczynniki prostej bazowej
const double aa =  1;
const double bb = -1;
// Przyrównanie do zera
const double EPS = 1e-12;

// Tablice
//--------

// Współrzędne danych
double x[N];
double y[N];

// Zwraca wartość pseudolosową
// w przedziale [-1;1]
//----------------------------
double random()
{
  return -1.0 + 2.0 * (double)rand() /
         RAND_MAX;
}

// Wyznacza losowe dane
//----------------------
void set_nodes()
{
  int i;
  for(i = 0; i < N; i++)
  {
    x[i] = Xmin + (Xmax - Xmin) *
           (double) rand() / RAND_MAX;
    y[i] = aa * x[i] + bb + random();
  }
}

// Program główny
//---------------
int main()
{
  int i;
  double a,b,Sx,Sy,Sxy,Sx2,M;

  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(4)
       << fixed;

  cout << "Liniowa aproksymacja "
          "średniokwadratowa\n"
          "---------------------"
          "-----------------\n\n";

  // Inicjujemy generator
  // pseudolosowy
  srand(time(nullptr));

  // Generujemy dane
  set_nodes();

  // Obliczamy współczynniki
  // prostej aproksymującej
  Sx = Sy = Sxy = Sx2 = 0;
  for(i = 0; i < N; i++)
  {
    Sx  += x[i];
    Sy  += y[i];
    Sxy += x[i] * y[i];
    Sx2 += x[i] * x[i];
  }

  // Liczymy mianownik
  M = N * Sx2 - Sx * Sx;

  // Testujemy mianownik
  if(fabs(M) < EPS)
    cout << "Błąd, dane w pionie";
  else
  {
    // liczymy a i b
    a = (N * Sxy - Sx * Sy) / M;
    b = (Sy - a * Sx) / N;

    // Wyświetlamy wynik porównawczy
    cout << "prosta bazowa       : "
            "y = " << aa
         << " * x + " << bb
         << endl << endl
         << "prosta aproksymowana: "
            "y = " << a
         << " * x + " << b
         << endl;
  }

  cout << endl;
  system("pause");
  return 0;
}

Liniowa aproksymacja średniokwadratowa
--------------------------------------

prosta bazowa       : y = 1.0000 * x + -1.0000

prosta aproksymowana: y = 1.0010 * x + -1.0218

Zauważ, iż po wyliczeniu z chmury danych prostej aproksymowanej \({y = ax + b}\) możemy przy jej pomocy wyliczać wartości leżące poza przedziałem \([Xmin,Xmax]\). To jest podstawowa różnica w stosunku do interpolacji. Im dalej od tego przedziału tym bardziej rośnie błąd aproksymacji.

# Liniowa aproksymacja
# średniokwadratowa
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0078
#---------------------------

from random import uniform

# Definicje stałych
#------------------

# Liczba danych
n = 100
# Przedział dla danych
xmin = -10
xmax =  10
# Współczynniki prostej bazowej
aa =  1.0
bb = -1.0
# Przyrównanie do zera
eps = 1e-12

# Tablice
#--------

# Współrzędne danych
x = [0.0] * n
y = [0.0] * n

# Wyznacza losowe dane
#----------------------
def set_nodes():
    for i in range(n):
        x[i] = uniform(xmin, xmax)
        y[i] = (aa * x[i] + bb
               + uniform(-1,1))

# Program główny
#---------------

print("Liniowa aproksymacja "
      "średniokwadratowa\n"
      "---------------------"
      "-----------------\n")

# Generujemy dane
set_nodes()

# Obliczamy współczynniki
# prostej aproksymującej
sx, sy, sxy, sx2 = 0, 0, 0, 0
for i in range(n):
    sx  += x[i]
    sy  += y[i]
    sxy += x[i] * y[i]
    sx2 += x[i] * x[i]

# Liczymy mianownik
m = n * sx2 - sx * sx

# Testujemy mianownik
if abs(m) < eps:
    print("Błąd, dane w pionie")
else:
    # liczymy a i b
    a = (n * sxy - sx * sy) / m
    b = (sy - a * sx) / n

    # Wyświetlamy wynik porównawczy
    print(f"prosta bazowa       : "
          f"y = {aa:.4f}"
          f" * x + {bb:.4f}\n"
          f"prosta aproksymowana: "
          f"y = {a:.4f}"
          f" * x + {b:.4f}")

print()
input("Naciśnij ENTER...")

Gdy dane wejściowe nie są rozłożone wzdłuż prostej, lecz raczej np. po łuku, wtedy należy użyć aproksymacji parabolicznej, która dopasuje do nich parabolę \({f(x) = ax^2 + bx + c}\) w podobny sposób, jak poprzednio została dopasowana funkcja liniowa \({f(x) = ax + b}\).

Na wejściu mamy zbiór danych w postaci punktów \({(x_i,y_i), i = 0,1, \dots , n - 1}\), do których będziemy dopasowywać odpowiednią parabolę. Definiujemy funkcję kryterialną, która jest sumą kwadratów pionowych odchyłek wyliczanej paraboli od punktów danych:

Szukamy minimum funkcji \({S(a,b,c)}\). Warunkiem koniecznym istnienia minimum jest zerowa wartość pochodnych cząstkowych tej funkcji względem zmiennych \(a\), \(b\) i \(c\).

Obliczamy pochodne cząstkowe, pamiętając o regule różniczkowania funkcji złożonej. Wnętrza nawiasów różniczkujemy odpowiednio względem \(a\): co daje \(-x_{i}^{2}\), względem \(b\): (\(-x_{i}\)) oraz względem \(c\): (\(-1\)):

Wszystkie trzy równania dzielimy obustronnie przez \(-2\), aby pozbyć się stałych czynników przed sumami:

Następnie wymnażamy wyrażenia wewnątrz sum przez czynniki stojące przed nawiasami:

Korzystając z liniowości operacji sumowania, rozbijamy powyższe równania na pojedyncze sumy, wyprowadzając szukane stałe \(a\), \(b\) oraz \(c\) przed znaki sum.

Porządkujemy układ, przenosząc wyrazy zawierające szukane niewiadome na lewą stronę, a sumy z wartościami \(y_i\) na prawą stronę. Otrzymujemy układ trzech równań liniowych z trzema niewiadomymi:

W celu uproszczenia wprowadźmy symbole pomocnicze:

Zapiszmy układ równań z wykorzystaniem tych symboli pomocniczych:

Rozwiązując powyższy układ równań liniowych jedną z opisanych wcześniej metod (np. Cramera, Gaussa, LU, itp.), wyznaczamy bezpośrednio poszukiwane współczynniki wielomianu aproksymującego: \(a\), \(b\) oraz \(c\).

Przykładowa implementacja

W celu przetestowania algorytmu parabolicznej aproksymacji średniokwadratowej poniższy program wykonuje następujące operacje:

• Określa:
  • \(N\) - liczbę punktów (np. 100)
  • przedział \([Xmin; Xmax]\)
  • współczynniki \(aa\), \(bb\) i \(c\) paraboli bazowej \({y = aax^2 + bbx + cc}\)
• Wyznacza losowo w tym przedziale \(N\) punktów \(x\) i dla każdego z nich liczy \({y = aax^2 + bbx + cc}\). Do wyniku dodaje losową wartość od -1 do 1, co tworzy odchyłkę i daje w wyniku punkty leżące ponad i pod parabolą bazową.
• Wg wzorów aproksymacji średniokwadratowej wylicza współczynniki \(a\) , \(b\) i \(c\) dla paraboli aproksymującej.
• Wyświetla wynik porównawczy: parabolę bazową i parabolę aproksymującą.

Wszystko jest opisane w komentarzach.