Propabilistische Programmierung

… hier wird noch gearbeitet …

Parameterberechnung mittels Wahrscheinlichkeiten

In dem Beitrag zur klassischen Ausgleichsrechnung (Seite verlinken) haben wir uns mit dem Finden von optimalen Parametern anhand einer klassischen mathematischen Methoden befasst. Dazu wurde das Ohm’sche Gesetz
\[ U = R \cdot I + U_0 \]
als Problem mit einem und später zwei Freiheitsgraden betrachtet, und aus gegebenen Testdaten die gesuchten Werte bestimmt. In diesem Artikel werden wir das dort beschriebene Modell mit zwei Freiheitsgraden verwenden, um statistische Methoden zur Parameterfindung zu verstehen.

Das Grundprinzip der propabilistischen Ansätze beruht auf dem Theorem von Bayes.
\[ P\left( \Theta | d \right) = \frac{P\left( d | \Theta \right) P\left( d \right)}{P\left( \Theta \right)} \]
Wenn man \( P\left( \Theta | d \right) \) als die Wahrscheinlichkeit versteht, dass \( \Theta \) aus \( d \) folgt und \( P\left( d | \Theta \right) \) als die, dass \( d\) aus \( \Theta \) folgt, dann stellt dieses Theorem eine verblüffend leichte Weise dar, die Logik „aus \( d \) folgt \( \Theta \)“ mittels Wahrscheinlichkeitsexperimenten zu invertieren.

Seien \( d \) die Messdaten der Messdatenreihe und \( \Theta \) die gesuchten Parameter, dann beschreibt \( P(d) \) die Wahrscheinlichkeit der Messdaten und \( P(\Theta) \) die der Parameter. Nach dem Theorem von Bayes lässt sich aus der Wahrscheinlichkeit mit der ein Parameter \( \Theta \) die Daten \( d \) bewirkt die Umkehrung direkt berechnen:
Die Wahrscheinlichkeit mit der die Daten \( d \) den Parameter \( \Theta \) bewirken.
In dieser Formel steht kein Differential. Sie stellt die direkte Berechnung der Wahrscheinlichkeit dar.

Die Herleitung des Theorems und die Umsetzung in nutzbare Ansätze liegt außerhalb unseres Fokus. Allein die Wahrnehmung der letzten Jahre, dass Datenanalyse eine stark gewachsene Rolle spielt, und das Etablieren von Methoden in gängigen Programmiersprachen genügen, um neugierig verstehen zu wollen, was diese Entwicklung für das Bestimmen von Modellparametern bedeutet.

Einer der Berechnungskerne, die entwickelt zur statistischen Datenanalyse entwicklet wurden, ist MC-Stan, das als PyStan in Python zur Verfügung steht. PyStan ist viel mächtiger als wir hier beleuchten können. Das folgende Beispiel ist aus einem einfachen Demo-Code von PyStan abgeleitet.

Stan beinhaltet eine eigene Spezifikationssprache, in der das zu untersuchende Problem beschrieben wird. Aus der Spezifiaktion wird wegen der effizienteren Laufzeit C-Code generiert, der als dynamische Library über eine Python-Schnittstelle angesprochen wird.

Die Spezifikation für unser Modell sieht wie folgt aus:


data {
    int N;      // Anzahl der gemessenen Punkte
    real U[N];           // Messdaten der Spannung
    real I[N];           // Messdaten des Stromes
    real sigma; // und eine Streubreute
}
parameters {
    // Hier die Parameter zum Anpassen
    real R;
    real U0;
}
transformed parameters {
    real U_berechnet[N];
    // Hier steht unser Modell
    for (j in 1:N)
      U_berechnet[j] = R * I[j] + U0; // das Ohm'sche Gesetz
}
model {
    // prior: erwartete Verteilung der gesuchten Parameter
    R ~ std_normal();   // R ist normalverteilt
    U0 ~ std_normal();  // U0 ist normalverteilt
    // likelihood: die Verknuepfung der berechneten Werte mit den gemessenen
    U ~ normal(U_berechnet, sigma); 
}

Man erkennt eine C-ähnliche Syntax. Die vier Blöcke „data“, „parameters“, „transformed parameters“ und „model“ enthalten Anweisungen, die mit „,“ getrennt sind.

  • „data“ spezifizier die Schnittstelle zu Variablen im Python-Programm. Eine Größe „N“ enthält die Anzahl der Messpunkte. Dieser folgen zwei Arrays, jeweils eines für die Messdaten des Stromes und der Spannung. Auch eine Streubreite „sigma“ mit auf positive Zahlen eingeschränktem Wertebereich wird hier deklariert. Später wird man sehen, wie deren Vorgabe aus dem Python-Programm erfolgt.
  • „parameters“ deklariert die zu suchenden Parameter „R“ und „U0“.
  • „transformed parameters“ enthält unser Modell. Dazu wird das Array „U_berechnet“ angelegt und mittels einer Vorschleife entsprechend dem Modell befüllt.
  • „model“ fixiert die Stochastischen Anforderungen an die Parameter: „R“ und „U0“ sollen normalverteilt sein. Außerdem wird das Array mit den Messdaten der Spannung auch über eine Normalverteilung mit dem mittels Modell berechneten Array verknüpft. Die hier verwendete Streubreite „sigma“ wird explizit aus Python heraus vorgegeben.

Diese Spezifikation findet sich im folgenden Python-Programm als String nach der üblichen Eröffnungssequenz mit Import der nötigen Module.

In der main-Abfrage finden sich die Messdaten und deren Aufbereitung in Numpy-Arrays, die die importierte pystan-Schnittstelle benötigt. In der Funktion „doMyStaff“ werden die Python-Seite und die Stan-Spezifikation miteinander verknüpft. Dann wird der Compiler aufgerufen und die Berechnung gestartet. In der pystan-Doku finden sich umfangreiche Erläuterungen und Beispiele dazu. 


#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.use("Agg") # force Matplotlib backend to Agg
import pystan

#
# Dies ist das Propabilistische Programm 'MC-Stan'
# 
line_code = """
data {
    int N;      // Anzahl der gemessenen Punkte
    real U[N];           // Messdaten der Spannung
    real I[N];           // Messdaten des Stromes
    real sigma; // und eine Streubreute
}
parameters {
    // Hier die Parameter zum Anpassen
    real R;
    real U0;
}
transformed parameters {
    real U_berechnet[N];
    // Hier steht unser Modell
    for (j in 1:N)
      U_berechnet[j] = R * I[j] + U0; // das Ohm'sche Gesetz
}
model {
    // prior: erwartete Verteilung der gesuchten Parameter
    R ~ std_normal();   // R ist normalverteilt
    U0 ~ std_normal();  // U0 ist normalverteilt
    // likelihood: die Verknuepfung der berechneten Werte mit den gemessenen
    U ~ normal(U_berechnet, sigma); 
}
"""



def doMyStaff(I,U):
    # Daten zwichen Python und MC-Stan verknüpfen
    linear_data = {'N': len(I),  # Anzahl der Punkte
                   'I': I,       # I Messdaten als numpy array
                   'U': U,       # U Messdaten als numpy array
                   'sigma': 1.00} # Startwert für sigma
    Nsamples = 1000000 # Die Lösung wird aus 10e6 Versuchen generiert
    chains = 1         # und das in einem Prozess (nicht mehrere)
    sm = pystan.StanModel(model_code=line_code) # Modell kompilieren
    #
    # und nun kommt die Anpassung
    fit = sm.sampling(data=linear_data, iter=Nsamples, chains=chains); 
    la = fit.extract(permuted=True)  # Das Ergebnis lesen

    #
    for k in la.keys():
        print("Name: {} Länge: {}".format(k,len(la[k])))
        print(la[k])
    postsamples = np.vstack((la['R'], la['U0'])).T
    # print(postsamples)
    print('Number of posterior samples is {}'.format(postsamples.shape[0]))
    # Das Ergebnis der Parameter wird in einem Corner-Plot dargestellt.
    try:
        import corner 
    except ImportError:
        sys.exit(1)
    fig = corner.corner(postsamples, labels=[r"$R$ in kOhm",r"$U0$ in V"])
    fig.savefig('stR_Plot_2.png')


if __name__ == '__main__':
    #
    # Unsre Daten
    #
    u=[-5.03,-2.99,-1.00, 0.99, 3.02, 4.98]
    i=[-2.03,-1.24,-0.41, 0.40, 1.21, 2.01]
    I=np.array(i)*0.001
    U=np.array(u)
    plt.plot(I*1000.0,U, 'o', color='tab:brown')
    plt.xlabel("Strom in mA")
    plt.ylabel("Spannung in V")
    #plt.show()
    plt.savefig("U_I.png", dpi=300)
    doMyStaff(I*1000,U)

Wie funktioniert dieses Programm?
Aus der Spezifikation, den Daten und der vorgegebenen Streubreite „sigma“ würfelt es solange Parameter „R“ und „U0“ bis 1000000 erfolgreiche Zuordnungen vorliegen. Also es würfelt solange, bis 1000000 mal die aus dem Modell berechneten Spannungen innerhalt der vorgegebenen Streubreite der gemessenen Spannungen liegen. Über diese erfolgreichen Daten-Parameter-Kombinationen erfolgt dann eine statistische Auswertung.

Falls das Modul „corner“ installiert ist entsteht mit „sigma=1.0″folgendes Bild.

sigma = 1.0 V
sigma = 0.02 V

Das Ergebnis mit kleinerer Streubreite „sigma“ weist in den Messdaten einen Offset aus. Deshalb könnte man das Modell auch bezüglich des Stromes um einen Offset erweitern.
\[ U = R \cdot \left( I +I_0 \right) + U_0 \]
In einer Vorlesung mit einigen aufmerksamen Studenten gäbe es jetzt Protest: Der Stromoffset müsse negativ verrechnet werden und das System habe einen Freiheitsgrad zuviel. Trotzdem: Wie reagiert der Propabilistische Ansatz?

Zunächst wieder die Auswertung mit einer Streubreite von 1.0V. Die drei Parameter und ihre wechselweise Kopplung zeigen Auffälligkeiten: Die Streubreiten von \( U_0 \) und \( I_0 \) sind erheblich größer und die Korrelation der Punkte zeigt sich in der mittleren Grafik in der untersten Zeile.

Bei einer vorgegeben Streubreite von 0,02 V ändern sich diese Auffälligkeiten nicht. Der Korrelationseffekt zeigt sich erheblich schärfer.

Was hätten und die aufmerksamen Zuhörer mitgeteilt? Das Modell kann umgeschrieben werden zu
\[ U = R \cdot I + \left(U_0 + R \cdot I_0 \right) \]
Nun sieht man die Lineare Abhängigkeit
\[ \hat{U_0} = U_0 + R \cdot I_0 \]
Wenn dieses \( \hat{U_0}\) konstant gehalten wird, erhält man immer gleichwertige Ergebnisse.

Hinweise für Ihre eigenen Experimente

Die Installation von pystan unter Linux erfolgte in meinen Versuchen problemlos. Unter Windows ergeben sich wegen der Abhängigkeit von einem Gnu-C-Compiler gcc Randbedingungen. Hier müssen Sie sehr genau auf die Installationsanleitung von MC-Stan achten.

Der in vielen propabilistischen Beispielen dargestellte Arbeitsstil ist interaktiv. Die Programmzeilen werden einzeln in interaktive Shells eingegeben. Damit werden viele Programmzeilen global angelegt. Dies kann beim Aktivieren mehrerer Prozesse, Chains, unter Windows zu unangenehmen Wechselwirkungen führen. Der obige Sourcecode kann hinsichtlich eines robusteren Verhaltens als Vorlage dienen.