Laboratorium 4: Aproksymacja

Wstęp

Podstawowym problemem interpolacji jest to, że stara się przeprowadzić funkcję przez wszystkie dane, które posiadamy (węzły). Ma to sens tylko i wyłącznie wtedy, gdy dane są dokładne i niezaburzone. Gdy jest inaczej — powinniśmy myśleć o jakimś ich uśrednieniu.

Aproksymacja

Aproksymacja to taka metoda „przybliżania" danych, w której zadaną funkcję stara się tak poprowadzić, żeby była jak najbliżej posiadanych punktów.

Osobną kwestią jest ustalenie co to znaczy „jak najbliżej"? Jeżeli mamy zestaw danych pomiarowych (par) $(x_i,y_i),;i=1,2,\dots ,N$, szukamy takiej funkcji $f(x)$ aby: $$Q=\sum _{i=1}^N(f(x_i)-y_i{)}^2\to min!$$

Tak postawione zadanie jest bardzo trudne — minimalizacja polega na wybraniu funkcji takiej, żeby… Znacznie prościej jest rozwiązywać zadanie następujące. Niech $f(x)=g(x,a)$ gdzie $a$ jest wektorem parametrów $a=(a_1,a_2,\dots ,a_M),;M\le N$ $$Q=\sum _{j=1}^N(g(x_j,a)-y_i{)}^2\to min!$$ Teraz zadanie optymalizacji jest łatwiejsze — musimy wybrać wektor liczb. Kolejne uproszczenie polega na rozważaniu zadanie liniowego względem parametrów: $$g(x,a)=\sum _{j=1}^M{a}_j{\phi }_j(x),$$ a zadanie optymalizacji wygląda tak: $$Q=\sum _{i=1}^N{(\sum _{j=1}^M{a}_j{\phi }_j(x_i)-y_i)}^2\to min!$$

Jego rozwiązanie jest stosunkowo proste — wystarczy wyliczyć pochodne cząstkowe $\frac{\partial Q}{\partial a_j}$ i rozwiązać układ równań: $$\frac{\partial Q}{\partial a_j}=0;j=1,2,\dots ,M$$

Zadanie dalej się upraszcza gdy przyjąć, że funkcja ${\phi }_j(x)=x^j$.

Aproksymacja a interpolacja

W przypadku zadania interpolacji żądamy, aby funkcja interpolująca przeszła przez wszystkie punkty (węzły interpolacyjne).

Poniżej przedstawiam zestaw punktów (pomiary temperatury termometrem IR).

Krzywe interpolacyjne mogą wyglądać tak jak na kolejnym rysunku. Czerwonymi kropkami zaznaczone są węzły interpolacji. Zwracam uwagę, że różnica między interpolacją Hermite’a a krzywymi sklejanymi nie jest specjalnie wielka. Niepokojąco natomiast wyglądają różnice pochodnych — pochodna interpolacji splajnami sześciennymi jest gładka.

Aproksymacja — Mathematica

Do realizacji aproksymacji wykorzystać można w Mathematici funkcję Fit. Jej wywołanie jest następujące:

Fit[data, funs, vars]

gdzie data to zestaw danych (par punktów), funs funkcja lub wektor funkcji którymi przybliżamy. Na przykład: $\{1,x,x^2,x^3,x^4,x^5,x^6,x^7,x^8,x^9,x^{10},x^{11},x^{12}\}$, vars — zmienna lub zmienne niezależne.

Powyższy zestaw jednomianów w różnych potęgach można łatwo wygenerować automatycznie: $$\text{funs}=\text{Table}[x^i,i,0,10]$$ $$\{1,x,x^2,x^3,x^4,x^5,x^6,x^7,x^8,x^9,x^{10}\}$$

i dalej: $$\text{funkcja1}=\text{Fit}[\text{dane}[[\text{All},2]],\text{funs},x]$$ w wyniku dostajemy współczynniki wielomianu: $$\begin{array}{rl}& -\text{8.003515809018834}\stackrel{`}{}\text{*}{}^{\wedge }\text{-20}{x}^{10}+\text{1.1265995371896338}\stackrel{`}{}\text{*}{}^{\wedge }\text{-16}{x}^9\\ & -\text{6.645297813196042}\stackrel{`}{}\text{*}{}^{\wedge }\text{-14}{x}^8+\text{2.1252203010076843}\stackrel{`}{}\text{*}{}^{\wedge }\text{-11}{x}^7\\ & -\text{3.981375997713063}\stackrel{`}{}\text{*}{}^{\wedge }\text{-9}{x}^6+\text{4.408253297181539}\stackrel{`}{}\text{*}{}^{\wedge }\text{-7}{x}^5\\ & -0.0000278354x^4+0.00093783x^3-0.0165472x^2+0.22931x-1.59072S\end{array}$$ Korzysta się z otrzymanej funkcji aproksymacyjnej dosyć łatwo, na przykład: $$\text{Plot}[\text{funkcja1},x,0,288,\text{PlotStyle}\to \text{Blue}]$$ albo $$\text{ff1}[\text{x\_}]=\text{funkcja1};$$ i $$\text{Plot}[\text{ff1}[x],\{x,0,288\}]\text{Plot}[\text{ff1}[x],\{x,0,288\}]$$

W przypadku bardziej skomplikowanych zadań wykorzystać można również funkcje:

• FindFit (funkcja aproksymująca nie musi liniowo zależeć od parametrów),
• LinearModelFit (tylko dla modeli liniowych) i chyba najogólniejszą:
• NonlinearModelFit.

Na poniższej ilustracji przykład aproksymacji dobowych zmian temperatury z termometru IR wielomianami stopnia 12 (zielony) i 4 (niebieski)).

Jak widać — cały problem sprowadza się do wyboru odpowiedniej funkcji aproksymacyjnej.

Matlab

Możliwości matlaba w zakresie aproksymacji wydają się być mniejsze. Toolbox Curve Fitting zawiera funkcję o nazwie fit i wywołaniu:

fit(x,y,fitType)

$x$ i $y$ to dane wejściowe. jako fitType podać należy łańcuch znaków określający rodzaj aproksymacji. Możliwości opisuje dokumentacja. Są tam wielomiany do stopnia 9 i parę innych funkcji.

Najprostsze użycie (korzystające z dostarczonych z matlabem danych przykładowych) wyglądać może tak:

load census;
f=fit(cdate,pop,'poly2')

census to jeden ze standardowych, przykładowych danych (w tym wypadku statystycznych) MATLABa zawierających informacje dotyczące liczby ludności USA w latach 1790–1990.

f = 
     Linear model Poly2:
     f(x) = p1*x^2 + p2*x + p3
     Coefficients (with 95% confidence bounds):
       p1 =    0.006541  (0.006124, 0.006958)
       p2 =      -23.51  (-25.09, -21.93)
       p3 =   2.113e+04  (1.964e+04, 2.262e+04)

Funkcja $f()$ to wielomian drugiego stopnia; funkcja podaje wartość współczynników oraz granice ich istotności.

plot(f,cdate,pop)

Aproksymacja a regresja

Wyobraźmy sobie, że mamy $n$ pomiarów $x_i$ jakiegoś parametru i chcemy zaaproksymować je wartością stałą $\overline{x}$. Chcielibyśmy, aby ta stała była jak najbliższa wszystkim pomiarom. Interesuje nas zatem taki problem: $$Q=\sum _{i-1}^n(x_i-\overline{x}{)}^2\to min!$$ czyli szukamy takiej wartości $\overline{x}$, która minimalizuje $Q$. Policzmy więc pierwszą pochodną $\frac{dQ}{d\overline{x}}$ (będziemy przyrównywać ją do zera):

$$\frac{dQ}{d\overline{x}}=\sum _{i=1}^{n}2(x_i-\overline{x})=2\sum _{i=1}^n{x}_i-2n\overline{x}=0$$

zatem

$$\overline{x}=\frac{1}{n}\sum _{i=1}^n{x}_i.$$

Wzór ten przypomina nam znany ze statystyki wzór na średnią.

Nie od rzeczy będzie wspomnieć, że aproksymacja ma bardzo wiele wspólnego ze znaną ze statystyki regresją. W pewnym sensie jest to to samo (choć nie należy mówić tego głośno) — w przypadku regresji jest cała otoczka związana z probabilistyką (w szczególności zakłada się, że $x_i$ są to obserwacje pewnej zmiennej losowej $X$, a to bardzo silne założenie — mówi ono, o tym, że istnieje rozkład prawdopodobieństwa zmiennej losowej $X$). Można w takim przypadku pokazać, że wyliczona wartość $\overline{x}$ ma pewne pożądane właściwości — wraz ze wzrostem $n$, $\overline{x}$ zmierza do wartości średniej rozkładu (jest estymatorem wartości średniej) i, że jest to estymator nieobciążony.

W technice wykorzystuje się średnią do „polepszania" wyników pomiarów. Zakładamy, że wartość $a$ mierzona jest z pewnym addytywnym błędem, czyli: $x_i=a+{\zeta }_i$; zaburzenia ${\zeta }_i$ są niezależnymi realizacjami obserwacji pewnej zmiennej losowej $Z$ o średniej 0. Zatem wyliczając $\sum _{i=1}^n{x}_i$, po dokonaniu odpowiednio wielu pomiarów „odkryć" możemy prawdziwą wartość $a$.

Podobne interpretacje można zaprezentować również dla innych zadań, w których stosujemy aproksymację.

Zadanie do wykonania

Wybrać jakiś przebieg dobowy i przybliżyć go za pomocą jakiejś sprytnej funkcji (która dobrze będzie oddawała istotę zmienności przebiegu.

Uwagi:

1. Dane otrzymane z pomiarów zawierają bardzo duże wartości współrzędnych $x$. Może to stanowić problem podczas aproksymacji. Stanowczo więc zalecam przesunięcie czasu do zera (to znaczy pierwszy pomiar dokonywany jest w chwili 0, a następne co 300 sekund). Można to osiągnąć tak:
   $\text{dane}=\text{Import}[\text{AVERAGE300.dat}];$
   $\text{xmin}=\text{[}Min[\text{dane}[[\text{All},1]]];$
   $\text{dane}[[\text{All},1]]=\text{dane}[[\text{All},1]]-\text{xmin}$
   Można też porównać otrzymane wyniki z wynikami aproksymacji tylko wartości $y$ ($x$ przyjmuje wartości 1,2,…):
   $funkcja1=Fit[dane[[All,2]],funs,x]$

2. Jeżeli chodzi o wybór funkcji — sugeruję zacząć od wielomianów. Teoretycznie, im wyższy stopień wielomianu — tym przybliżenie lepsze. Tylko nie wiadomo czy sensowniejsze. Ambitni mogą wymyślić jakąś funkcję nieliniową lub złożyć z kawałków (patrz tutorial).