Metode najmanjih kvadrata

Često se pri raznim mjerenjima pojavljuje veliki broj podataka. Da bi lakše protumačili te podatke pribjegava se grafičkom predstavljanju tih podataka. Na primjer pri izvodjenju nekog eksperimenta dobijeni su podaci $$(x_0,y_0),\,(x_1,y_1),\dots,(x_7,y_7).$$

Ovakve podatke moguće je predstaviti grafikom, na primjer kao na sljedećoj slici.

Podaci iz prethodne tabele predstavljeni grafički

Na prvi pogled grafik linearne funkcije najbolje bi odgovara ovim eksperimentalnim podacima. Sljedeći korak bi bio da odredimo linearnu funkciju

$$y=ax+b,$$

odnosno koeficijente $a$ i $b.$ Postavlja se pitanje: Koja je to prava koja je najbliže prolazi predstavljenim tačkama $(x_0,y_0)-(x_7,y_7)?$ ili u opštem slučaju tačkama $(x_0,y_0)-(x_n,y_n)?$ Pretpostavimo da su $a$ i $b$ tražene vrijednosti. Neka je $(x_k,y_k)$ jedna od datih tačaka i neka kroz ovu tačku ne prolazi prava $y=ax+b,$ tada je $y_k\neq ax_k+b,$ tj. $y_k-a x_k-b\neq 0.$

Da bi izračunali $a$ i $b$ posmatrajmo funkciju $$\varphi(a,b)=\sum_{k=0}^{n}(ax_k+b-y_k)^2. \:(\star)$$

Ako greška ima normalnu statističku raspodjelu, onda minimizacijom funkcije $\varphi$ dobijamo najbolje procjena za $a$ i $b.$ Računanje parametara $a$ i $b$ na ovakav način, odnosno linearne funkcije kojom ćemo aproksimirati naziva se $l_2$ aproksimacija.

Napomena. Norma $l_p$ definisana je na sljedeći način $$ || x ||_p=\left\{ \sum_{i=1}^{n}|x_i|^p \right\}^{1/p},\:1\leqslant p<\infty,$$ za neki vektor $x=(x_1,\ldots,x_n)^T.$

Odredimo sada minimum funkcije $\varphi.$ Potrebni uslovi za minimum funkcije $\varphi$ su

$$ \frac{\partial \varphi}{\partial a}=0,\quad \frac{\partial \varphi}{\partial b}=0,$$

računajući prethodne uslove dobijamo

\begin{align*} \sum_{k=0}^{n}{2(ax_k+b-y_k)x_k}=&0\\ \sum_{k=0}^{n}{2(ax_k+b-y_k)}=&0 , \end{align*}

Prethodne jednačine su po nepoznatim $a$ i $b$ i možemo ih napisati u obliku

\begin{align} \left(\sum_{k=0}^{n}{x^2_k}\right)a+ \left(\sum_{k=0}^{n}{x_k}\right)b=&\sum_{k=0}^{n}{x_k y_k}\\ \left(\sum_{k=0}^{n}{x_k}\right)a+ (n+1)b=&\sum_{k=0}^{n}{y_k}, \end{align}

gdje je $\sum_{k=0}^{n}{1}=n+1.$

Rješenje prethodnog sistema je

\begin{align*} a=& \frac{1}{d}\left[(n+1)\sum_{k=0}^{n}{x_k y_k}-\left(\sum_{k=0}^{n}{x_k}\right)\left(\sum_{k=0}^{n} {y_k}\right) \right] \\ b=&\frac{1}{d}\left[\left(\sum_{k=0}^{n} {x^2_k}\right)\left(\sum_{k=0}^{n} {y_k}\right) -\left(\sum_{k=0}^{n}{x_k}\right)\left(\sum_{k=0}^{n} {x_ky_k}\right) \right], \end{align*}

gdje je $$d=(n+1)\sum_{k=0}^{n}{x^2_k}-\left(\sum_{k=0}^{n}{x_k^2 }\right)^2.$$

Zadatak. Izračunati za podatke $(0.5,2.25),(1,3.7),\,(2,4.1),\,(2.5,4.3),\,(3,5);$ linearnu funkciju $y=ax+b $ metodom najmanjih kvadrata.

In [1]:
from IPython import get_ipython;  # briše varijable, funkcije ...
get_ipython().magic('reset -sf')  #

import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import pyplot
from matplotlib import rc

pyplot.style.use('tableau-colorblind10')
newparams={'figure.figsize':(15,8),'font.size': 13,'axes.grid': True,'lines.markersize':10, 'lines.linewidth':2,'figure.dpi':300}
pyplot.rcParams.update(newparams)

xdata=[.5,1,2,2.5,3]
ydata=[2.25,3.7,4.1,4.3,5]

#xdata=np.linspace(1,15,200)
#def f1(x):
#    return 4*np.sin(3*x)*np.cos(14*x)+x
#ydata=f1(xdata)
broj=len(xdata)

def lin(x):
    xzbir=0; yzbir=0; xyzbir=0; xkvadrat=0; d=0; a=0; b=0;   
    for i in range(broj):
        xzbir=xzbir+xdata[i]
        yzbir=yzbir+ydata[i]
        xyzbir=xyzbir+xdata[i]*ydata[i]
        xkvadrat=xkvadrat+(xdata[i])**2
    
    d=(broj)*xkvadrat-xzbir**2    
    a=((broj)*xyzbir-xzbir*yzbir)/d
    b=(xkvadrat*yzbir-xzbir*xyzbir)/d
    print('Koeficijenti su a=%f i'%a, 'b=%f.'%b)
    return a*x+b

iks=np.linspace(np.min(xdata)-.25,np.max(xdata)+.25,1000)
fig=pyplot.figure()
pyplot.plot(xdata,ydata,'o',markersize=10)
pyplot.plot(iks,lin(iks),linewidth=3)
pyplot.legend(('Podaci','Linearna funkcija'))
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/najmanji2a.pdf')
pyplot.close(fig)
pyplot.show()
Koeficijenti su a=0.917442 i b=2.218605.
Konstrukcija linearne funkcije za uzorak od 5 tačaka

Konstrukcija linearne funkcije za uzorak od 200 tačaka

Generalni slučaj

Umjesto računanja linearne funkcije $y=ax+b,$ metoda najmanjih kvadrata može se poopštiti na računanje funkcije $$y=\sum_{j=0}^{m}{c_jg_j(x)}$$ kojom ćemo predstaviti podatke koji su nam na raspolaganju. Funkcije $g_0,\ldots,g_m$ su poznate, fiksirane i linearno nezavisne. Koeficijente $c_0,\ldots,c_m$ ponovo treba izračunati koristeći prethodno izloženi princip najmanjih kvadrata. Drugim riječima, ponovo su poznate kordinate tačaka $(x_k,y_k),\,i=0,\ldots,n;$ koje smo dobili mjerenjem, eksperimentalno ili na neki drugi način, i cilj nam je izračunati koeficijente $c_0,\ldots,c_m;$ tako da funkcija
\begin{equation} \varphi(c_0,\ldots,c_m)=\sum_{k=0}^{n}{\left[ \sum_{j=0}^{m}{c_jg_j(x_k)-y_k}\right]^2 } \label{funkcija1} \end{equation} ima što je moguće manju vrijednost. Potrebni uslovi za dobijanje minimuma funkcije $\varphi$ su $$ \frac{\partial \varphi}{\partial c_i}=0,\:i=0,\ldots,m$$ gdje su $$\frac{\partial \varphi}{\partial c_i}=\sum_{k=0}^{n}{2\left[\sum_{j=0}^{m}{c_jg_j(x_k)-y_k} \right]}g_i(x_k),\:i=0,\ldots,m.$$ Na kraju dobijamo $$ \sum_{j=0}^{m}{\left[ \sum_{k=0}^{n}{g_i(x_k)g_j(x_k)}\right]c_j}= \sum_{k=0}^{n}{y_kg_i(x_k)},\:i=0,\ldots,m.$$

Zadatak. Izračunati za podatke $(-1.5,2.3),\, (-1,1.1),\,(-0.5,0.24),\,(0,0.01),\,(0.5,0.26),\,(1,1.2),\,(1.5,2.25)$ kvadratnu funkciju $y=c_2x^2+c_1x+c_0$ metodom najmanjih kvadrata.

In [8]:
from IPython import get_ipython;  # briše varijable, funkcije ...
get_ipython().magic('reset -sf')  #

import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import rc

pyplot.style.use('tableau-colorblind10')
newparams={'figure.figsize':(15,8),'font.size': 13,'axes.grid': True,'lines.markersize':10, 'lines.linewidth':2,'figure.dpi':300}
pyplot.rcParams.update(newparams)

xdata=np.array([-1.5,-1,-.5,0,.5 ,1,1.5])
ydata=np.array([2.3,1.1,.24,.01,.26,1.2,2.25])

brojTacaka=len(xdata)
brojBaznih=3

matricaBaznih=np.zeros((brojBaznih,brojTacaka))
matricaSistema=np.zeros((brojBaznih,brojBaznih))
vektorSlob=np.zeros(brojBaznih)

def bazne(x,stepen):
    return x**stepen

for i in range(brojBaznih):
    for j in range(brojTacaka):
        matricaBaznih[i,j]=bazne(xdata[j],i)
        
for i in range(brojBaznih):
    for j in range(brojBaznih):
        for k in range(brojTacaka):
            matricaSistema[i,j]=matricaSistema[i,j]+matricaBaznih[i,k]*matricaBaznih[j,k]

for i in range(brojBaznih):
    for k in range(brojTacaka):
        vektorSlob[i]=vektorSlob[i]+ydata[k]*matricaBaznih[i,k]      
        
koeficijenti=linalg.solve(matricaSistema,vektorSlob)        

pol=np.polynomial.Polynomial(koeficijenti)

iks=np.linspace(min(xdata),max(xdata),250)
y=pol(iks)
fig=pyplot.subplot()
pyplot.plot(xdata,ydata,'o',markersize=10)
pyplot.plot(iks,y)  
pyplot.legend(('Podaci','polinom $p_{10}$'),loc='upper center');
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/najmanji4d.png')
#pyplot.xlabel(' Broj tačaka je %i'%len(xdata));
#pyplot.show()
pyplot.close()
Konstrukcija polinoma $p_2(x)=c_2x²+c_1x+c_0$ za uzorak od 7 tačaka
Konstrukcija polinoma $p_2(x)=c_2x²+c_1x+c_0$ za uzorak od 200 tačaka
Konstrukcija polinoma $p_3(x)=c_3x^3+c_2x^2+c_1x+c_0$ za uzorak od 200 tačaka
Konstrukcija polinoma $p_4(x)=c_4x^4+c_3x^3+c_2x^2+c_1x+c_0$ za uzorak od 200 tačaka
Konstrukcija polinoma $p_5(x)=c_5x^5+c_4x^4+c_3x^3+c_2x^2+c_1x+c_0$ za uzorak od 200 tačaka
Konstrukcija polinoma $p_{10}$ za uzorak od 500 tačaka

Napomena 1. Kao bazne funkcije umjesto $1,x,\ldots,x^n$ mogu se koristiti i druge funkcije.


Napomena 2. Ovdje je predstavljena najjednostavija metoda. Druge metode i detaljnija analiza može se vidjeti na primjer u [J97] ili [TB97].

Literatura

[CD04] W. Cheney, D. Kincaid, Numerical Mathematics and Computing, Brooks/Cole-Thompson Learning, USA, 2004. ISBN: 053489937.

[J97] J.W. Demmel, Applied Numerical Linear Algebra, SIAM, USA, 1997. ISBN: 0-89871-389-7.

[TB97] L.N.Trefethen, D.Bau III, Numerical Linear Algebra, SIAM, USA, 1997. ISBN: 0-89871-361-7.