XProfan

WindowTitle "Polynom-Regression mit der Methode der kleinsten Quadrate"
'---------------------------------------------------
' VisualBasic-Original (C) Prof.em. Rolf HIRTHE, https://www.rhirte.de/vb/home.htm
' Übertragen nach XProfan-11.2a 2014-10 by P.Specht, Wien; Ohne jede Gewähr!
'---------------------------------------------------
WindowStyle 24:Window 0,0-%maxx,%maxy-40:randomize:set("decimals",16):font 0
print "\n\n Die Gauss'sche Methode der minimierten mittleren quadratischen Abweichung"
print " wird hier dazu verwendet, eine Polynomkurve zwischen Messpunkte einzupassen. "
print " Ziel ist, durch Anpassung des Polynomgrades trotz Messfehlern die dem "
print " Punktverlauf best-angepasste Funktion als 'Gesetzmäßigkeit' zu finden.\n"
Declare a![50,51],b![50,51],x![50],y![50],Xv![25],Yv![25],grad&,paar&
'---------------------------------------------------
' Einige Messpunkte P(xv,yv) als Testbeispiel:
'---------------------------------------------------
paar&=5
Xv![1] = 1   : Yv![1]= -0.3
Xv![2] = 1.9 : Yv![2] = 0.3
Xv![3] = 2.9 : Yv![3] = 1.3
Xv![4] = 4   : Yv![4] = 3.3
Xv![5] = 5.1 : Yv![5] = 4.3
'---------------------------------------------------
grad& = 2' des Polynoms
'---------------------------------------------------
font 2
Print "\n\n Anpassung eines Polynoms "+str$(grad&)+". Grades"
print " an ";paar&;" Messpunkte (Einprogrammmiertes Beispiel)\n\n"
'---------------------------------------------------
print AUSGLEICH(paar&,grad&,Xv![],Yv![])
'---------------------------------------------------
'Hier potuto eine grafische Ausgabe den Kurvenverlauf zeigen
waitinput
END
'===================================================

Proc AUSGLEICH :parameters paar&,Pgrad&,Xv![],Yv![]

    declare i&,j&,k&,Ggrad&,Mat![25,25],Fak!,A![25],sig!,txt$

    If paar& < (Pgrad& + 1)'Vorbereitende Eingaben

        txt$=" FEHLER: Zahl der Wertepaare nicht ausreichend!"
        goto "Exit_Sub"

    EndIf

    Ggrad&=Pgrad&+1

    Whileloop Ggrad&:i&=&Loop

        WhileLoop Ggrad&+1:j&=&Loop

            Mat![i&,j&]=0

        endwhile

    endwhile

    '---------------------------------------------------
    Mat![1,1]=paar&' Aufbau der Matrix

    WhileLoop paar&:i&=&Loop

        Mat![1,Ggrad&+1]=Mat![1,Ggrad&+1]+Yv![i&]
        Fak!=1

        Whileloop 2,Ggrad&:j&=&Loop

            Fak! = Fak! * Xv![i&]
            Mat![1,j&]=Mat![1,j&]+Fak!
            Mat![j&,Ggrad&+1] = Mat![j&,Ggrad&+1]+Fak!*Yv![i&]

        Endwhile

        Whileloop 2,Ggrad&:j&=&Loop

            Fak! = Fak! * Xv![i&]
            Mat![j&,Ggrad&]=Mat![j&,Ggrad&]+Fak!

        endwhile

    endwhile

    whileloop 2,Ggrad&:i&=&Loop

        WhileLoop Ggrad&-1:j&=&Loop

            Mat![i&,j&]=Mat![i&-1,j&+1]

        endwhile

    endwhile

    '---------------------------------------------------
    ' Aufruf der Routine zur Lösung des lin. GS
    GAUSS Ggrad&,Mat![],a![]
    '---------------------------------------------------
    txt$=""'Ergebnisanzeige

    WhileLoop Ggrad&:i&=&Loop

        txt$=txt$+" a"+str$(int(i&-1))+" = "+str$(a![i&])+"\n"

    endwhile

    txt$=txt$+"\n"
    sig!=0
    '---------------------------------------------------
    ' Berechnung und Anzeige zur Qualität der Anpassung
    '---------------------------------------------------

    whileloop paar&:i&=&Loop

        Fak! = a![Ggrad&]

        WhileLoop Ggrad&-1,1,-1:j&=&Loop

            Fak! = a![j&] + Xv![i&] * Fak!

        endwhile

        sig!=sig!+sqr(Yv![i&]-Fak!)
        txt$=txt$+" y("+str$(i&)+")="+str$(Yv![i&])+" -> "+format$("%g",Fak!)+"\n"

    endwhile

    txt$=txt$+"\n Durchschnittliche Abweichung [dsig] = " + str$(Sqrt(sig!/(paar&-2)))+"\n"
    '---------------------------------------------------
    Exit_Sub:
    return txt$

EndProc

Proc Gauss :parameters n&,A![],x![]

    declare i&,j&,k&,jmax&,kmax&,merk&[]
    declare s!,max!,skal![]
    setsize merk&[],n&
    setsize skal![],n&
    '---------------------------------------------------
    ' Reihenfolge sichern

    whileloop n&:i&=&Loop

        merk&[i&] = i&

    endwhile

    '---------------------------------------------------
    ' Normalisierung

    whileloop n&:i&=&Loop

        s! = 0

        whileloop n&:j&=&Loop

            s! = s! + Abs(A![i&,j&])

        endwhile

        skal![i&]=1/s!

    endwhile

    '---------------------------------------------------
    ' Vorwärtselimination

    whileloop n&-1:k&=&Loop

        max! = skal![k&]*Abs(A![k&,k&])
        kmax& = k&'Spalte mit max
        jmax& = k&'Zeile mit max
        'Pivotzelle suchen:

        whileloop k&,n&:j&=&Loop

            whileloop k&,n&:i&=&Loop

                If (skal![j&]*Abs(A![j&,i&]))> max!

                    jmax& = j&
                    kmax& = i&
                    max! = skal![j&]*Abs(A![j&,i&])

                EndIf

            endwhile

        endwhile

        '---------------------------------------------------

        If jmax& <> k&' Zeilentausch, wenn nötig

            whileloop k&,n&+1:j&=&Loop

                s! = A![k&,j&]
                A![k&,j&] = A![jmax&,j&]
                A![jmax&,j&] = s!

            endwhile

            s! = skal![k&]
            skal![k&] = skal![jmax&]
            skal![jmax&] = s!

        EndIf

        '---------------------------------------------------

        If kmax& <> k&'Spaltentausch, wenn nötig

            whileloop n&:i&=&Loop

                s! = A![i&,k&]
                A![i&,k&] = A![i&,kmax&]
                A![i&,kmax&] = s!

            endwhile

            j& = merk&[k&]
            merk&[k&] = merk&[kmax&]
            merk&[kmax&] = j&

        EndIf

        '---------------------------------------------------
        ' eigentliche Elimination

        whileloop k&+1,n&:i&=&Loop

            s! = A![i&,k&]/A![k&,k&]
            A![i&,k&]=0

            whileloop k&+1,n&+1:j&=&Loop

                A![i&,j&] = A![i&,j&] - s! * A![k&,j&]

            endwhile

        endwhile

    endwhile

    '---------------------------------------------------
    ' Auflösung rückwärts
    x![merk&[n&]] = A![n&,n&+1]/A![n&,n&]

    whileloop n&-1,1,-1:i&=&Loop

        s! = A![i&,n&+1]

        whileloop i&+1,n&:j&=&Loop

            s! = s! - A![i&,j&]*x![merk&[j&]]

        endwhile

        x![merk&[i&]]=s!/A![i&,i&]

    endwhile

Endproc