Das deutsche Python-Forum

Hallo Zusammen

Ich habe mich von pydeal's Idee inspirieren lassen und möchte mal mein kleines Wochenendergebnis hinsichtlich des Codes sowie hinsichtlich des Ergebnisses zur Diskussion stellen. Mich jedenfalls haut ein Standardfehler von 16.25 % nicht gerade vom Hocker.

Code: Alles auswählen

from matplotlib import pyplot as plt

def determinante(matrix = []):
    summe = 0
    zeilen = len(matrix)
    for liste in matrix:
        if len(liste) != zeilen:
            print("Keine quadratische Matrix!")
            return
    if zeilen == 2:
        return matrix[0][0]*matrix[1][1]-matrix[1][0]*matrix[0][1]
    elif zeilen < 2:
        print("Benötigt werden mindestens 2 Zeilen und 2 Spalten.")
        return
    else:
        index = 0
        vorzeichen = 1
        for zahl in matrix[0]:
            if vorzeichen:
                summe += zahl*determinante(bildeSubmatrix(matrix,index))
                vorzeichen = 0
            else:
                summe -= zahl*determinante(bildeSubmatrix(matrix,index))
                vorzeichen = 1
            index += 1
    return summe

def bildeSubmatrix(matrix = [], index = 0):
    submatrix = []
    for i in range(1, len(matrix)):
        zeile = []
        for j in range(len(matrix[i])):
            if j != index:
                zeile.append(matrix[i][j])
        submatrix.append(zeile)
    return submatrix
    
def cramerscheRegel(koeffizientenmatrix = [], loesungsvektor = []):
    zeilen = len(koeffizientenmatrix)
    for liste in koeffizientenmatrix:
        if len(liste) != zeilen:
            print("Keine quadratische Matrix!")
            return
    nenner = determinante(koeffizientenmatrix)
    if nenner == 0.0:
        print("Gleichungssystem nicht eindeutig loesbar!")
        return
    loesung = []
    for i in range(len(koeffizientenmatrix[0])):
        tmpMatrix = [item*1 for item in koeffizientenmatrix]
        j = 0
        for zeile in tmpMatrix:
            zeile[i] = loesungsvektor[j]
            j += 1
        loesung.append(determinante(tmpMatrix)/nenner)
    return loesung

if __name__ == "__main__":
    
    Dax = [0.0956,0.0265,0.2548,0.2906,-0.1469,
           0.1606,0.2385,-0.4037,0.2229,0.2198,
           0.2707,0.0734,0.3708,-0.4394,-0.1979]
    Gold = [-0.1211,0.0012,-0.2733,0.0826,0.0893,
            0.2924,0.2504,0.0432,0.3192,0.2320,
            0.1777,0.0465,0.1989,0.2557,0.0141]
    Jahre = []
    jahr = 2015
    while jahr > 2001:
        Jahre.append(jahr)
        jahr -= 1
    Jahre.append(jahr)

    xWerteProzente = [item*100 for item in Dax]
    yWerteProzente = [item*100 for item in Gold]

    summex = sum(xWerteProzente)
    summey = sum(yWerteProzente)
    summexy = 0
    summex2 = 0
    summex3 = 0
    summex4 = 0
    summex2y = 0
    for i in range(len(xWerteProzente)):
        summexy += xWerteProzente[i]*yWerteProzente[i]
        summex2 += xWerteProzente[i]**2
        summex3 += xWerteProzente[i]**3
        summex4 += xWerteProzente[i]**4
        summex2y += (xWerteProzente[i]**2)*yWerteProzente[i]

    koeffMatrix = [[len(xWerteProzente),summex,summex2],
                   [summex,summex2,summex3],
                   [summex2,summex3,summex4]]
    lVektor = [summey,summexy,summex2y]

    funktionswerte = cramerscheRegel(koeffMatrix,lVektor)

    xWerteFunktion = [item for item in xWerteProzente]
    yWerteFunktion = []
    xWerteFunktion.sort()

    for item in xWerteFunktion:
        y = funktionswerte[0]+(funktionswerte[1]*item)+(funktionswerte[2]*(item**2))
        yWerteFunktion.append(y)

    quadrResiduum = 0
    for i in range(len(xWerteProzente)):
        quadrResiduum += ((yWerteProzente[i])
                          -(funktionswerte[0]
                            +(funktionswerte[1]*xWerteProzente[i])
                            +(funktionswerte[2]*(xWerteProzente[i]**2))))**2

    standarderror = (quadrResiduum/(len(yWerteProzente)-2))**0.5

    yWerteFunktionStdErrUp = []
    yWerteFunktionStdErrDwn = []
    for item in yWerteFunktion:
        yWerteFunktionStdErrUp.append(item+standarderror)
        yWerteFunktionStdErrDwn.append(item-standarderror)

    stringTitelOben = "Regression Dax/Gold [2001-2015] annual performance\n"
    stringTitelUnten = "".join(["(",str(round(funktionswerte[0],2)),
                               "+",str(round(funktionswerte[1],2)),
                               "x+",str(round(funktionswerte[2],2)),
                               "x^2) +/- standard error (",
                               str(round(standarderror,2)),")"])

    plt.grid(color='grey', linestyle='-', linewidth=0.25)
    plt.plot(xWerteFunktion,yWerteFunktion,"r",
             xWerteFunktion,yWerteFunktionStdErrUp,"r--",
             xWerteFunktion,yWerteFunktionStdErrDwn,"r--",
             xWerteProzente,yWerteProzente,"go")
    
    text = []
    for i in range(len(Jahre)):
        plt.text(xWerteProzente[i],yWerteProzente[i],str(Jahre[i]))
    
    plt.title(stringTitelOben+stringTitelUnten)
    plt.xlabel("Dax")
    plt.ylabel("Gold")
    plt.show()

Ich freue mich auf Eure Beiträge.

@Bob Billsworth: Fast die Hälfte des Codes läuft im Namensraum des Moduls und legt dort Variablen an. Das Hauptprogramm steht üblicherweise in einer Funktion die `main()` heisst.

Viele Namensschreibweisen entsprechen nicht dem Style Guide for Python Code.

Veränderbare Defaultargumente sind ein „code smell“ weil die Ausdrücke dafür nur *einmal* ausgewertet werden, wenn die Funktion definiert wird, und nicht bei jedem Aufruf.

Allerdings machen die Defaultwerte hier auch allesamt überhaupt keinen Sinn. Die Funktionen können nicht sinnvoll ohne Matrixargumente aufgerufen werden. Noch unsinniger wird es wenn der Defaultwert *falsch* ist. Was soll das?

``for i in range(len(sequence)):`` um dann `i` als Index in das Sequenzobjekt zu verwenden ist in Python ein „anti pattern“. Man kann *direkt* über die Elemente von Sequenzen iterieren. Und falls man den Index *zusätzlich* zum Element benötigt, gibt es die `enumerate()`-Funktion. Zum ”parallelen” iterieren über mehr als eine Sequenz gibt es die `zip()`-Funktion (und `itertools.izip()` in Python 2). Auch dafür braucht man keine Indexvariable. Du hast da ja fast alle Möglichkeiten durchprobiert die man in Python *nicht* schreiben würde — Iterieren über Index, Index manuell in der ``for``-Schleife hochgezählt, und sogar eine ``while``-Schleife in der eine Zählvariable manuell aktualisiert wird.

`bilde_submatrix()` kann man mit „slicing“ und einer „list comprehension“ wesentlich kompakter schreiben, so dass das weniger nach einer anderen Programmiersprache aussieht. Das ist dann ein Einzeiler.

Eine `print()`-Ausgabe und ein ``return`` ohne expliziten Rückgabewert sind keine sinvolle Art mit Fehlern in Funktionen umzugehen. Das führt zu einer Ausnahme an anderer Stelle weil der Aufrufer nichts mit dem `None`-Wert anfangen kann. Sinnvoller wäre es selbst eine Ausnahme auszulösen die neben dem Text auch einen passenden Typ hat.

`zeilen` in `determinante()` ist ein mehrdeutiger Name. Die meisten würden hinter dem Namen tatsächlich die Zeilen erwarten und nur deren Anzahl. `liste` ist ein sehr/zu generischer Name wenn damit doch eine Zeile der Matrix gemeint ist.

`summe` wird nur in einem ``else``-Zweig verwendet, aber schon ganz am Anfang der Funktion definiert.

`vorzeichen` ist eigentlich ein Wahrheitswert und keine Zahl und darum sollte man das auch so belegen. Alternieren kann man das dann mit ``not``.

Der Code für den Test ob eine Matrix quadratisch ist, kommt zweimal im Programm vor. Das sollte man in eine Funktion herausziehen. Weder Code noch Daten sollte man wiederholen, das verletzt das DRY-Prinzip („Don't Repeat Yourself“).

Eine Liste mit eins zu multiplizieren um eine (flache) Kopie zu erstellen ist ein netter Trick für „code golf“, hat in regulären Programmen aber nichts zu suchen.

`Jahre` wird fast am Anfang des Hauptprogramms definiert, aber erst fast am Ende verwendet. Das liesse sich auch deutlich einfacher erstellen — es gibt `range()`, `reversed()` und `list()`. Damit wäre das ein Einzeiler. Letztendlich will man hier ja aber Jahre von 2015 an rückwärts zählen, was mit `itertools.count()` geht.

Für eine flache Kopie einer Liste benötigt man keine „list comprehension“ (LC), da reicht ein Aufruf von `list()`. Falls man eine (flache) sortierte Kopie benötigt, gibt es die `sorted()`-Funktion.

Python kennt Zeichenkettenformatierung. Werte und Zeichenketten mit `str()` und `round()` und ``+`` oder ``''.join(…)`` zusammenzubasteln ist nicht wirklich gut lesbar und eher Basic als Python.

Das Hauptprogramm steht auf Modulebene und hinterlässt viele Namen direkt im Modul.

Ich lande dann ungefähr bei so etwas:

Code: Alles auswählen

#!/usr/bin/env python
# coding: utf8
from __future__ import absolute_import, division, print_function
from itertools import count, cycle, izip
from matplotlib import pyplot as plt


def test_auf_quadratisch(matrix):
    zeilenanzahl = len(matrix)
    if any(len(zeile) != zeilenanzahl for zeile in matrix):
        raise ValueError('Keine quadratische Matrix!')


def bilde_submatrix(matrix, index):
    return [zeile[:index] + zeile[index + 1:] for zeile in matrix[1:]]


def determinante(matrix):
    test_auf_quadratisch(matrix)
    zeilenanzahl = len(matrix)
    if zeilenanzahl < 2:
        raise ValueError('Benötigt werden mindestens 2 Zeilen und 2 Spalten.')
    
    if zeilenanzahl == 2:
        return matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1]
    else:
        return sum(
            x * determinante(bilde_submatrix(matrix, i)) * factor
            for i, (x, factor) in enumerate(izip(matrix[0], cycle([1, -1])))
        )


def cramersche_regel(koeffizientenmatrix, loesungsvektor):
    test_auf_quadratisch(koeffizientenmatrix)
    nenner = determinante(koeffizientenmatrix)
    if nenner == 0.0:
        raise ValueError('Gleichungssystem nicht eindeutig loesbar!')

    loesung = list()
    for i in xrange(len(koeffizientenmatrix[0])):
        temporaere_matrix = map(list, koeffizientenmatrix)
        for j, zeile in enumerate(temporaere_matrix):
            zeile[i] = loesungsvektor[j]
        loesung.append(determinante(temporaere_matrix) / nenner)
    return loesung


def main():
    # 
    # Jahr für den ersten Wert in den Wertelisten.  Bei den Folgewerten
    # ist das Jahr absteigend.
    # 
    endjahr = 2015
    dax = [
        0.0956, 0.0265, 0.2548, 0.2906, -0.1469, 0.1606, 0.2385, -0.4037,
        0.2229, 0.2198, 0.2707, 0.0734, 0.3708, -0.4394, -0.1979,
    ]
    gold = [
        -0.1211, 0.0012, -0.2733, 0.0826, 0.0893, 0.2924, 0.2504, 0.0432,
        0.3192, 0.2320, 0.1777, 0.0465, 0.1989, 0.2557, 0.0141,
    ]

    dax_prozente = [wert * 100 for wert in dax]
    gold_prozente = [wert * 100 for wert in gold]

    summe_x = sum(dax_prozente)
    summe_y = sum(gold_prozente)
    summe_xy = 0
    summe_x2 = 0
    summe_x3 = 0
    summe_x4 = 0
    summe_x2y = 0
    for dax_prozent, gold_prozent in izip(dax_prozente, gold_prozente):
        summe_xy += dax_prozent * gold_prozent
        summe_x2 += dax_prozent**2
        summe_x3 += dax_prozent**3
        summe_x4 += dax_prozent**4
        summe_x2y += dax_prozent**2 * gold_prozent

    koeffizientenmatrix = [
        [len(dax_prozente), summe_x, summe_x2],
        [summe_x, summe_x2, summe_x3],
        [summe_x2, summe_x3, summe_x4],
    ]
    loesungsvektor = [summe_y, summe_xy, summe_x2y]

    funktionswerte = cramersche_regel(koeffizientenmatrix, loesungsvektor)

    sortierte_dax_prozente = sorted(dax_prozente)
    y_werte = [
        (
            funktionswerte[0]
            + (funktionswerte[1] * wert)
            + (funktionswerte[2] * wert**2)
        )
        for wert in sortierte_dax_prozente
    ]

    quadr_residuum = 0
    for dax_prozent, gold_prozent in izip(dax_prozente, gold_prozente):
        quadr_residuum += (
            gold_prozent - (
                funktionswerte[0]
                + (funktionswerte[1] * dax_prozent)
                + (funktionswerte[2] * dax_prozent**2)
            )
        )**2
    standard_error = (quadr_residuum / (len(gold_prozente) - 2))**0.5

    y_werte_std_err_up = list()
    y_werte_std_err_down = list()
    for wert in y_werte:
        y_werte_std_err_up.append(wert + standard_error)
        y_werte_std_err_down.append(wert - standard_error)

    titel = (
        u'Regression Dax/Gold [{0}-{1}] annual performance\n'
        u'(${2[0]:.2f}+{2[1]:.2f}x+{2[2]:.2f}x^2$) ± standard error ({3:.2f})'
            .format(
                endjahr - len(dax_prozente) + 1,
                endjahr,
                funktionswerte,
                standard_error,
            )
    )

    plt.grid(color='grey', linestyle='-', linewidth=0.25)
    plt.plot(
        sortierte_dax_prozente, y_werte, 'r',
        sortierte_dax_prozente, y_werte_std_err_up, 'r--',
        sortierte_dax_prozente, y_werte_std_err_down, 'r--',
        dax_prozente, gold_prozente, 'go'
    )

    for dax_prozent, gold_prozent, jahr in izip(
        dax_prozente, gold_prozente, count(endjahr, -1)
    ):
        plt.text(dax_prozent, gold_prozent, str(jahr))
   
    plt.title(titel)
    plt.xlabel('Dax')
    plt.ylabel('Gold')
    plt.show()


if __name__ == '__main__':
    main()

Das Hauptprogramm ist auch noch sehr lang, das kann man sicher noch sinnvoll auf Funktionen aufteilen.

Da man mit `matplotlib` sowieso schon Numpy als Abhängigkeit hat, sollte man das vielleicht auch nutzen, statt das Rad neu zu erfinden. Da fallen dann gleich zwei komplette Funktionen weg, denn für die Determinante gibt es bereits `numpy.linalg.det()`. Und man kann sich ein paar Schleifen sparen:

Code: Alles auswählen

#!/usr/bin/env python
# coding: utf8
from __future__ import absolute_import, division, print_function
from itertools import count, cycle, izip
from matplotlib import pyplot as plt
import numpy as np


def cramersche_regel(koeffizientenmatrix, loesungsvektor):
    nenner = np.linalg.det(koeffizientenmatrix)
    if nenner == 0.0:
        raise ValueError('Gleichungssystem nicht eindeutig loesbar!')

    loesung = np.zeros(len(koeffizientenmatrix))
    for i in xrange(len(koeffizientenmatrix)):
        temporaere_matrix = koeffizientenmatrix.copy()
        temporaere_matrix[i, :] = loesungsvektor
        loesung[i] = np.linalg.det(temporaere_matrix)
    return loesung / nenner


def main():
    # 
    # Jahr für den ersten Wert in den Wertelisten.  Bei den Folgewerten
    # ist das Jahr absteigend.
    # 
    endjahr = 2015
    dax = np.array(
        [
            0.0956, 0.0265, 0.2548, 0.2906, -0.1469, 0.1606, 0.2385, -0.4037,
            0.2229, 0.2198, 0.2707, 0.0734, 0.3708, -0.4394, -0.1979,
        ]
    )
    gold = np.array(
        [
            -0.1211, 0.0012, -0.2733, 0.0826, 0.0893, 0.2924, 0.2504, 0.0432,
            0.3192, 0.2320, 0.1777, 0.0465, 0.1989, 0.2557, 0.0141,
        ]
    )

    dax_prozente = dax * 100
    gold_prozente = gold * 100

    summe_x = dax_prozente.sum()
    summe_y = gold_prozente.sum()
    summe_xy = (dax_prozente * gold_prozente).sum()
    summe_x2 = (dax_prozente**2).sum()
    summe_x3 = (dax_prozente**3).sum()
    summe_x4 = (dax_prozente**4).sum()
    summe_x2y = (dax_prozente**2 * gold_prozente).sum()

    koeffizientenmatrix = np.array(
        [
            [len(dax_prozente), summe_x, summe_x2],
            [summe_x, summe_x2, summe_x3],
            [summe_x2, summe_x3, summe_x4],
        ]
    )
    loesungsvektor = np.array([summe_y, summe_xy, summe_x2y])

    funktionswerte = cramersche_regel(koeffizientenmatrix, loesungsvektor)

    sortierte_dax_prozente = dax_prozente[dax_prozente.argsort()]
    y_werte = (
        funktionswerte[0]
        + (funktionswerte[1] * sortierte_dax_prozente)
        + (funktionswerte[2] * sortierte_dax_prozente**2)
    )

    quadr_residuum = (
        (
            gold_prozente
            - (
                funktionswerte[0]
                + (funktionswerte[1] * dax_prozente)
                + (funktionswerte[2] * dax_prozente**2)
            )
        )**2
    ).sum()
    standard_error = (quadr_residuum / (len(gold_prozente) - 2))**0.5

    y_werte_std_err_up = y_werte + standard_error
    y_werte_std_err_down = y_werte - standard_error

    titel = (
        u'Regression Dax/Gold [{0}-{1}] annual performance\n'
        u'(${2[0]:.2f}+{2[1]:.2f}x+{2[2]:.2f}x^2$) ± standard error ({3:.2f})'
            .format(
                endjahr - len(dax_prozente) + 1,
                endjahr,
                funktionswerte,
                standard_error,
            )
    )

    plt.grid(color='grey', linestyle='-', linewidth=0.25)
    plt.plot(
        sortierte_dax_prozente, y_werte, 'r',
        sortierte_dax_prozente, y_werte_std_err_up, 'r--',
        sortierte_dax_prozente, y_werte_std_err_down, 'r--',
        dax_prozente, gold_prozente, 'go'
    )

    for dax_prozent, gold_prozent, jahr in izip(
        dax_prozente, gold_prozente, count(endjahr, -1)
    ):
        plt.text(dax_prozent, gold_prozent, str(jahr))
   
    plt.title(titel)
    plt.xlabel('Dax')
    plt.ylabel('Gold')
    plt.show()


if __name__ == '__main__':
    main()

Hey BlackJack

Vielen Dank für Deinen ausführlichen Beitrag. Ich werde ihn verinnerlichen.

Das deutsche Python-Forum

Nichtlineare Regression/Cramer'sche Regel

Nichtlineare Regression/Cramer'sche Regel

Re: Nichtlineare Regression/Cramer'sche Regel

Re: Nichtlineare Regression/Cramer'sche Regel