Bepalen paraboolfunctie met reeks en top

[Marcoe]

Beste forummers,

Mijn probleem is het bepalen van een 2e-graads functie uit een reeks x en y waarden met een bekende top. Het bepalen van de top wordt mathematisch uitgevoerd door het bepalen van de beste regressiecoeffcient.(er zijn mogelijk betere manieren om dit te doen, maar daar kom ik eventueel later nog op terug). Het aantal meetparen is variabel van 3 tot maximaal 9. De randvoorwaarde voor de functie is dat het altijd een bergparabool betreft, de meetparen stellen het gedeelte links van de parabooltop voor. In totaal zijn het altijd 10 meetparen, dat wil dus zeggen dat als de top bepaald is de paren rechts van de top (1 tot maximaal 7 paren) verwijderd worden.

Het bepalen van de functie kan ik indirect doen door per meetpaar en de top de functie te bepalen en vervolgens de beste fit te kiezen met behulp van de regressiecoefficient. Dit is niet helemaal juist en ik vraag me daarom af hoe ik dit wel goed kan doen.
Onderstaand een voorbeeldje met x en y-waarden:

Code: Selecteer alles

0.0025	21.95	
0.0075	22.32	
0.0125	22.67	
0.0175	23.24	
0.0225	23.42	
0.0275	23.43	
0.0325	23.46	
0.0375	23.31	
0.0425	23.62	
0.0475	23.36	

Alvast bedankt.

Mvg, Marco

[arie]

Via lineaire regressie kan je de coëfficiënten a, b en c van je parabool bepalen.
Je krijgt dan de functie
y = a*x^2 + b*x + c
die het beste fit bij je meetwaarden.
Zie bijvoorbeeld http://en.wikipedia.org/wiki/Regression ... near_model

In jouw voorbeeld geldt voor de formule



op die wiki-pagina dat:
Y = 10*1 matrix (ofwel vector) met je y-waarden
X = 10*3 matrix, met voor elke regel i: X[i,0]=1, X[i,1]=xi, X[i,2]=xi^2

De best passende fit is dan:



In je voorbeeld kom ik dan uit op
a = B2 = -1475.7575757576
b = B1 = 104.890909090909
c = B0 = 21.682450757576

met de top bij x = 0.03553798767967

Als je de 3 linker waarden minder betrouwbaar vindt, kan je die weggooien,
vervolgens de overgebleven 7 punten spiegelen in de lijn x = 0.03553798767967,
en dan opnieuw de lineaire regressie toepassen op de 14 punten die je nu hebt.

Zoek je zoiets?

[Marcoe]

Ja dit is wat ik zoek, zal wel eerst weer aan de studie moeten om het rekenen met een matrix onder de knie te krijgen.

[arie]

De theorie lijkt me nog te overzien:
het gaat in ons geval om de getransponeerde (T), de inverse (-1) en het matrixproduct.

Voor de praktijk:
Er zijn veel rekenprogramma's die deze berekeningen kunnen maken.
Ook in Excel lijkt het al te lukken, zie bv
http://web.iyte.edu.tr/food/Course_Web_ ... 0Excel.pdf

[Marcoe]

Hi Arie,

Thanks voor je hulp, weer een stapje verder.

De regressiecoefficienten a, b en c kun je in excel berekenen met de matrix functie: lijnsch

Code: Selecteer alles

={LIJNSCH(B2:B11,A2:A11^{1,2})}

Daaruit kan ik dan de top berekenen met -b/2a.

Het punt is dat niet alle waardes geldig zijn, met name de laatste 7 punten komen mogelijk te vervallen. Dit is afhankelijk van de beste fit in de punten voor de top. De moeilijkheid zit juist in de eerste drie punten, omdat daar de grootste onnauwkeurigheid zit.

[Marcoe]

Ik heb nu door hoe je met een matrix kan rekenen. Alleen begrijp ik nog niet hoe je met een vaste top de beste fit kan bepalen (bijvoorbeeld het laatstse punt uit de tabel als top)

[arie]

Bedenk eerst of je zeker weet dat je dat echt wilt. Het punt dat je als top uitkiest is immers ook een meetwaarde (met een (onbekende) meetfout erin).

Indien je de top van je parabool vastlegt, wordt de hele berekening wel veel eenvoudiger.
Neem bijvoorbeeld de eerste 7 metingen uit je voorbeeld:

Code: Selecteer alles

x        y   
0.0025   21.95   
0.0075   22.32   
0.0125   22.67   
0.0175   23.24   
0.0225   23.42   
0.0275   23.43   
0.0325   23.46   

waarbij je het laatste punt als top T = (xt,yt) kiest, dus xt=0.0325 en yt=23.46.

Transleer alle punten nu zo dat de top in de oorsprong O = (0,0) komt te liggen:
trek van alle punten de coördinaten van de top af, je krijgt dan voor elk van de oude punten (xi, yi) een nieuw punt (Xi,Yi) waarbij:
Xi = xi - xt
Yi = yi - yt

Je zoekt nu de best passende parabool
f(X) = a*X^2 + b*X + c

Omdat deze parabool door de oorsprong gaat moet f(0) = 0 zijn, dus c = 0.
Omdat de top ligt bij X=0, moet ook -b/2a gelijk zijn aan 0, zodat ook b = 0.

De best passende parabool heeft nu dus de vorm:
f(X) = a*X^2
ofwel



Nu zoeken we de a, zodanig dat











minimaal is. Dan moet de afgeleide van E(a) naar a nul zijn:



ofwel



dus geheel zonder matrices te berekenen.

In ons voorbeeld kom ik uit op a = -1717.274725... :

Code: Selecteer alles

xi      yi      Xi       Yi     Xi^2 * Yi   Xi^4
0.0025  21.95  -0.0300  -1.51  -0.00135900  0.0000008100
0.0075  22.32  -0.0250  -1.14  -0.00071250  0.0000003906
0.0125  22.67  -0.0200  -0.79  -0.00031600  0.0000001600
0.0175  23.24  -0.0150  -0.22  -0.00004950  0.0000000506
0.0225  23.42  -0.0100  -0.04  -0.00000400  0.0000000100
0.0275  23.43  -0.0050  -0.03  -0.00000075  0.0000000006
0.0325  23.46   0.0000   0.00   0.00000000  0.0000000000
                               -----------  ------------
                               -0.00244175  0.0000014219

waardoor a = -0.00244175 / 0.0000014219 = -1717.274725

Gebruik tenslotte weer
X = x - xt
Y = y - yt
om de parabool terug te transleren naar de oorspronkelijke positie van je puntenverzameling.

Dit levert in ons geval:

[Marcoe]

Chapeau, ik was even vergeten hoe leuk wiskunde kan zijn.

Je hebt gelijk, een vaste top is toch niet zo'n goed idee. Ik heb nu genoeg data om vast te stellen dat de meetfout te groot is om de beste fit te bepalen met een vaste top.

Arie bedankt!

grt,
Marco

[Marcoe]

Weer een stapje verder.

Nu ik alle meetreeksen verwerkt heb kom ik erachter dat de methode van het te bepalen aantal bruikbare meetparen mbv de maximale determinatiecoefficient niet goed genoeg is. Als voorwaarde heb ik er al ingezet dat het laatste bruikbare meetpaar links van top moet liggen. Als nl bij een van de meetreeksen een mogelijke uitbijter zit hoe haal ik deze uitbijter dan eruit, of is het mogelijk de waarde te corrigeren.
Met behulp van een residuplot kan ik wel de spreiding van de meetpunten t.o.v. het geschatte model bepalen, maar ik kan nog geen uitspraak doen over welk punt nu een uitbijter is of niet.

[arie]

Hier staat wat informatie over uitbijters: http://nl.wikipedia.org/wiki/Uitbijter.
Doorgaans verwijder je uitbijters uit je dataset omdat corrigeren niet mogelijk is (zie het voorbeeld met de kamertemperatuur op die wiki pagina: de waarde 200 graden is duidelijk een uitbijter, maar daaruit kunnen we niet afleiden wat de werkelijke waarde is).

Bedenk dat we hier bovendien werken met een kwadratisch model: .
Als je werkt met een foutmaat epsilon in de vorm
kan je al je residuen direct vergelijken.
Als je echter hebt

met een fout delta klein t.o.v. x, dan is

en is de afwijking in y_i door delta ook afhankelijk van x en a.
Je zou je residuen daar dan nog op moeten corrigeren.

Wellicht dat een statisticus of iemand met meer ervaring op dit gebied hier meer over kan vertellen.