Pienimmän neliösumman sovitus

V2 19.3. 2002

Aluksi luentoesimerkki

KRE esim. 1 s. 915

> restart:

Warning, the name changecoords has been redefined

> with(LinearAlgebra):alias(Tr=Transpose):

> xd:=[-1.3,-0.1,0.2,1.3]; yd:=[0.103,1.099,0.808,1.897];

> n:=nops(xd);

> plot([seq([xd[i],yd[i]],i=1..n)],view=[-2..2,0..2],style=point,symbol=circle,symbolsize=17,axes=box);

> datakuva:=%:

> C:=<<1,1,1,1>|<op(xd)>>;

> M:=Tr(C).C;

> B:=Tr(C).Vector(yd);

> ab:=LinearSolve(M,B);

> suora:=ab[1]+ab[2]*x;

> display(plot(suora,x=-2..2,color=blue),datakuva);

Muistellaan aluksi interpolaatiotehtävää.

Esim:

> xd:=[$1..10];yd:=[1.3,3.5,4.2,5.0,7.0,8.8,10.1,12.5,13.0,15.6];

> plot([seq([xd[i],yd[i]],i=1..10)],style=point,symbol=cross,symbolsize=18);

> p:=interp(xd,yd,x);

> display(%%,plot(p,x=0..10,color=blue));

Interpolaatio pakottaa polynomin usein voimakkaaseen heilahteluun, jolloin "metsä katoaa näkyvistä puilta".

Luovutaan vaatimuksesta, että polynomin tulisi kulkea datapisteiden kautta. Valitaan alempiasteinen polynomi, joka sopivassa mielessä approksimoi parhaiten annettua dataa.

Kriteeriksi valitaan usein funktiomallin (tässä polynomi) ja datapisteissä laskettujen erotuksten neliösumman minimointi.

Tätä kautta menetelmä liittyy funktioiden ääriarvoteemaan.

Polynomi-PNS-tehtävässä on annettu datapisteitä ( ) (paljon) ja tehtävänä on sovittaa esim. polynomi

annettuun dataan niin, että minimoituu.

Yleensä datapisteiden lkm =n on paljon suurempi kuin polynomin asteluku m (tässä 3).

Jos A on matriisi , jonka alkiot ovat , i=1..n, j=0..3, saadaan b-kertoimien määräämiseksi ylimääräytyvä yhtälösysteemi A b = y, missä siis A on n x 3 - matriisi.

PNS-ratkaisuun päästään kertomalla tämä systeemi A:n transpoosilla, jolloin ratkaistavaksi tulee

systeemi

(1) (A^T)A b = A^T y (missä A^T tarkoittaa A:n transpoosia)

(Tämä johdetaan joko osittaisderivoimalla tai ortogonaaliprojektioon perustuvalla geometrisluonteisella

argumentilla.)

Jos tuntemattomia kertoimia on yhtä monta kuin datapisteitä (polynomin asteluku = n-1), niin systeemillä

A b = y on yksikäsitteinen ratkaisu, mikäli x-pisteet ovat erillisiä. Tällöin on kyse interpolaatiopolynomista.

Yllä oleva ratkaisu voidaan Maplella tehdä rakentamalla matriisi A ja ratkaisemalla systeemi (1) .

Maplen linalg :ssa on myös leastsqrs - funktio, jolle annetaan pelkkä A - matriisi ja datapisteet, joka

siis ratkaisee ylimääräytyvän systeemin A b = y PNS-mielessä.

Tämän lisäksi Maplen stats -pakkauksessa on vielä automaattisempi funktio fit[leastsquare] , joka tekee

kaiken puolestamme (siis ei tarvitse muodostaa edes kerroinmatriisia).

Suosimme kahta ensin mainittua tapaa ensinnäkin opin kannalta, mutta myös siksi, että viimemainittu on syntaksiltaan outo ilmestys.

>

> restart:with(LinearAlgebra):with(linalg):with(plots):#with(stats):

Warning, the name changecoords has been redefined

Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value

Warning, the protected names norm and trace have been redefined and unprotected

> alias(Tr=Transpose):

Esim. 1

Muodosta PNS-suora annettuun pisteistöön.

> Plista:=[[2.8,30],[2.9,26],[3,33],[3.1,31],[3.2,33],[3.2,35],[3.2,37],[3.3,36],[3.4,33]];

> Matrix(Plista);

> plot(Plista);

> gp:=plot(Plista,style=POINT,symbol=circle):

> xd:=map(lista->lista[1],Plista);
n:=nops(xd):yd:=map(lista->lista[2],Plista);

>

> A:=<Vector([1$n])| Vector(xd)>;

Jos haluttaisiin jatkaa sarakkeittain latomista (ja siis korkeampiasteisten polynomien sovitusta), voitaisiin edetä tähän tapaan: (Alempana on kuitenkin yleispätevämpiä keinoja.)

> <Vector([1$n])| Vector(xd) | Vector(map(t->t^2,xd))>;

Tehdään aluksi "teoriaa" mukaillen:

>

> kertoimet:=LinearSolve(Tr(A).A,Tr(A).Vector(yd));

Suoremmin saadaan näin:

> leastsqrs(A,Vector(yd));

> suora:=innerprod(kertoimet,<1,x>);

> gs:=plot(suora,x=2.8..3.5,color=blue):

>

> display({gs,gp});

Esim. 2

Sovita a) 2. asteen PNS-polynomi ja b) interpolaatiopolynomi

> P:=[[1.4,7400],[1.8,7500],[2.4,7600],[3,7500],[4,7200]];

> xd:=map(lis->lis[1],P);
yd:=map(lis->lis[2],P);

Jos tuo näyttää oudolta, voitaisiin tehdä funktiot poimieka:=x->x[1]; poimitoka:=x->x[2]; ja suorittaa xd:=map(poimieka,P);yd:=map(poimitoka,P);

Tässä on kätevä ja yleispätevä tapa muodostaa matriisi.

> A:=matrix(5,3,(i,j)->xd[i]^(j-1));

Toisaalta yleissivistykseen kuuluu myös tietää, että kyseessä on ns. Vandermonden matriisi. Siihenkin

on linalg:ssa valmis funktio.

> Matrix(vandermonde(xd)):%[1..5,1..3]; # Vanhat tietorakenteet muutettava uusiksi (Matrix), jotta kätevät poimimispiirteet yms. toimisivat.

Käytetään nyt sitä keskimmäistä "harmaata laatikkoa"

> kertoimet:=leastsqrs(A,yd);

> poly2:=innerprod(kertoimet,[1,x,x^2]);

> display(plot(P,style=point,symbol=circle,color=black),plot(poly2,x=1.3..4.2));

>

b)

Kyseessä on interpolaatioply, jonka asteluku = datapisteiden lkm -1 = 4. No siihen löytyy taas ihan

valmiskin funktio:

(Jos ratkaistaisiin LSQ-tehtävänä, matriisina olisi koko Vandermonde (transposilla kertominen voitaisiin "supistaa pois").)

> intpol:=interp(xd,yd,x);

> display(plot(intpol,x=1.3..4.2,color=blue),plot(P,style=point,symbol=circle,color=black),plot(poly2,x=1.3..4.2,color=gray));

>

SVD:n käyttö LSQ:ssa