Lineare Approximation von vielen (x,y)-Wertpaaren

[Froggerprogger]

Der folgende Code nutzt eine Liste von n beliebigen (x,y)-Wertpaaren und berechnet daraus die Parameter a und b einer Geraden y = a*x + b, so dass die durchschnittliche quadratische Abweichung der Punkte von dieser minimal ist. Das nennt sich auch lineare Approximation oder lineare Ausgleichsrechnung. Weitere Infos hier:
http://de.wikipedia.org/wiki/Methode_de ... n_Quadrate

Die berechnete Gerade und die Punkte werden abschließend in ein Fenster gemalt.

Code: Alles auswählen

EnableExplicit

Structure Pos2D
  x.f
  y.f
EndStructure

NewList A.Pos2D()

; fill the list with some points that nearly fit the line y = a' * x + b'
Define i.l
For i=1 To 10
  AddElement(A())
  A()\x = Random(50000) / 100.0 - 250.0
  A()\y = -0.5*A()\x - 50 + (Random(4000)/100-20.0); a'=-0.5, b'= -50 plus random distub
Next

; now calculate a and b of the approximating line
Define n.f = 0 ; counts number of elements in list (same result as ListSize(A()))
Define avgX.d = 0 ; the average x-value
Define avgY.d = 0 ; the average y-value
Define sumDividend.d = 0 ; the cumulated dividend
Define sumDivisor.d = 0 ; the cumulated divisor
Define a.f, b.f ; the calculated a and b

ForEach A()
  sumDividend + (A()\x * A()\y)
  sumDivisor + (A()\x * A()\x)
  avgX + A()\x
  avgY + A()\y
  n+1.0
Next
avgX = avgX / n
avgY = avgY / n

; now we have all intermediate results for calculating a and b
a.f = (sumDividend - (n * avgX * avgY)) / (sumDivisor - (n * avgX * avgX))
b.f = avgY - a*avgX

; we're finished.
;End

; let's draw the output:
InitSprite()
OpenWindow(0, 0, 0, 600, 600, "Approximation by a line (a=" + StrF(a) + " / b=" + StrF(b) + ")")
OpenWindowedScreen(WindowID(0), 0, 0, 600, 600, 1, 0, 0)
Repeat
  StartDrawing(ScreenOutput())  
    ; draw x- and y-axis  
    Line(300, 0, 0, 600, $999999) 
    Line(0, 300, 600, 0, $999999)
    ; draw points
    ForEach A()
      Circle(A()\x + 300, A()\y + 300, 2, $0000FF)
    Next
    ; draw approximating line
    Define yFrom.f = a * -300 + b;
    Define yTo.f = a * 300 + b;
    Line(0, yFrom+300, 600, yTo - yFrom, $00FF00)
  StopDrawing()
  FlipBuffers()
  
  Define event = WaitWindowEvent()
Until event = #PB_Event_CloseWindow