[PYTHON] [Statistik] Lassen Sie uns die Ausführung der logistischen Regression in Stan im Detail erklären (mit Titanic-Datensatz).

In diesem Artikel wird Stan verwendet, um Titanic-Daten logistisch zurückzugeben und die Leistung der Klassifizierung weiter zu bewerten.

Die in diesem Artikel verwendete probabilistische Programmiersprache "Stan" verwendet die Hamilton-Monte-Carlo-Methode (HMC-Methode) und NUTS, um die Verteilungsparameter zu schätzen. Genau genommen ist das Prinzip der Erzeugung von Zufallszahlen unterschiedlich, aber eine etwas einfachere Methode ist die Markov-Ketten-Monte-Carlo-Methode, die Metropolis-Hastings-Methode (MH-Methode). Ich habe @ kenmatsu4 über dieses Funktionsprinzip geschrieben

Es gibt zwei Punkte. Bitte beziehen Sie sich auf diese, wenn Sie möchten. Ich denke, dass die MH-Methode und die HMC-Methode nicht so unterschiedlich sind, wenn die Absicht besteht, ein Bild davon zu geben, was sie tun. (Grundsätzlich der Unterschied in der Effizienz der Zufallsstichprobe)

Umgebung

OSX Yosemite 10.10.5
Python 2.7
Anaconda 3.18.9
PyStan 2.8.0

Der vollständige Code wird auf [GitHub] veröffentlicht (https://github.com/matsuken92/Qiita_Contents/blob/master/PyStan-Titanic/PyStan-Logistic-titanic-rev03.ipynb).

Einführung von PyStan

Ich hatte Angst, dass die Dinge nicht gut liefen, aber das war einfach und funktionierte.

conda install pystan

Bibliotheksvorbereitung

import numpy as np
import numpy.random as rd
import pandas as pd

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import animation as ani
import matplotlib.cm as cm
import seaborn as sns
sns.set(style="whitegrid", palette="muted", color_codes=True)

from tabulate import tabulate
from time import time

import pystan
from pystan import StanModel

from sklearn import cross_validation
from scipy.stats import gaussian_kde

Datenaufbereitung

[Maschinelles Lernen] Zusammenfassung und Ausführung der Modellbewertung / des Index (mit Titanic-Datensatz) Verwendet die Titanic-Daten, die in numerische Werte konvertiert wurden. Machen. Wir haben es auf GitHub vorbereitet, also benutze es bitte auch.

titanic = pd.read_table("https://raw.githubusercontent.com/matsuken92/Qiita_Contents/master/PyStan-Titanic/data/titanic_converted.csv",
              sep=",", header=0)
titanic.head(10)

ID	survived	pclass	sex	age	sibsp	parch	fare	embarked	class	who	adult_male	deck	embark_town	alone
0	0	3	1	22	1	0	7.25	1	3	1	1	0	3	1
1	1	1	0	38	1	0	71.2833	2	1	2	0	3	1	0
2	1	3	0	26	0	0	7.925	1	3	2	0	0	3	0
3	1	1	0	35	1	0	53.1	1	1	2	0	3	3	0
4	0	3	1	35	0	0	8.05	1	3	1	1	0	3	1
5	0	3	1	28	0	0	8.4583	3	3	1	1	0	2	1
6	0	1	1	54	0	0	51.8625	1	1	1	1	5	3	1
7	0	3	1	2	3	1	21.075	1	3	0	0	0	3	0
8	1	3	0	27	0	2	11.1333	1	3	2	0	0	3	0
9	1	2	0	14	1	0	30.0708	2	2	0	0	0	1	0

Teilen Sie die Daten in die erklärende Variable "Daten" und die Zielvariable "Ziel".

target = titanic.ix[:, 0]
data = titanic.ix[:, [1,2,3,4,5,6]]

#Trainingsdaten(80%),Testdaten(20%)Teilen in
X_train, X_test, y_train, y_test = cross_validation.train_test_split(data, target, test_size=0.2, random_state=0)
print [d.shape for d in [X_train, X_test, y_train, y_test]]

`out`


[(712, 6), (179, 6), (712,), (179,)]

Das heißt, die Trainingsdaten $ X $ sind

Es ist eine Matrix mit der Anzahl der Daten $ N = 712 $ und der Anzahl der Merkmale $ M = 6 $. Der horizontale Vektor $ x_i $ ist die Daten für eine Person. Die sechs Merkmalsgrößen sind pclass 　　sex　　 age　　 sibsp　　 parch　　 fare ist.

Die Lehrerdaten sind die Überlebensinformationen der Titanic. Überleben = 1 und die Person, die nie zurückgegeben hat = 0 $ N $ Dimensionsvektor

y= (y_1, \cdots, y_N)^{{\rm T}}

ist.

Über logistische Regression

Wenn es solche Daten gibt, wird das [logistische Regressionsmodell (Wikipedeia-Link)](https://ja.wikipedia.org/wiki/Logistic Regression) wie folgt ausgedrückt.

\operatorname{logit} (p_{i}) = \ln \left({\frac {p_{i}}{1-p_{i}}}\right)
= \beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} \\
(i=1, 2, \cdots, N)

$ p_i $ kann als Überlebenswahrscheinlichkeit der $ i $ -ten Datensatzdaten angesehen werden. Lassen Sie uns die Formel ein wenig umstellen, um zu sehen, warum sie als Überlebenswahrscheinlichkeit angesehen werden kann. Legen Sie zuerst beide Seiten in $ \ exp (\ cdot) $

\left({\frac {p_{i}}{1-p_{i}}}\right) = \exp ( \beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} )

Löse dies für $ p_i $

p_i = {1 \over 1 + \exp (\ -(\beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} )\ ) }

Hier die Standard-Sigmoid-Funktion $ \sigma(x) = {1 \over 1+ e^{-x}} $

Benutzen. Ähnlich wie in der vorherigen Formel, wobei $ x $ durch $ (\ beta + \ beta_ {1} x_ {1, i} + \ cdots + \ beta_ {M} x_ {M, i}) $ ersetzt wird Ich tat

p_i = \sigma(\beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i})

Sie können es als behandeln.

In einem Diagramm sieht es so aus.

def sigmoid(x):
    return 1 / (1+np.exp(-x))

x_range = np.linspace(-10, 10, 501)

plt.figure(figsize=(9,5))
plt.ylim(-0.1, 1.1)
plt.xlim(-10, 10)

plt.plot([-10,10],[0,0], "k", lw=1)
plt.plot([0,0],[-1,1.5], "k", lw=1)
plt.plot(x_range, sigmoid(x_range), zorder=100)

Wie Sie in der Grafik sehen können, nimmt die Sigmoidfunktion $ \ sigma (x) $ einen Wert im Bereich (0, 1) an, sodass dieser Wert als Wahrscheinlichkeit interpretiert und verwendet werden kann.

Von oben, der Teil, der wie eine normale lineare Regression aussieht,

\beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i}

Durch die Verwendung der Sigmoid-Funktion gut basierend auf

\operatorname{logit} (p_{i}) = \ln \left({\frac {p_{i}}{1-p_{i}}}\right)
= \beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} \\
(i=1, 2, \cdots, N)

Man kann sagen, dass logistische Regression die Form hat und den Wert als Wahrscheinlichkeit interpretiert.

Stan wird eingeführt, indem versucht wird, die Parameter $ \ beta, \ beta_1, \ cdots, \ beta_M $ dieses Ausdrucks dieses Mal mit Stan zu schätzen.

Laufen auf Stan

Erstellen Sie ein Wörterbuchobjekt mit Trainingsdaten für die Übergabe an Stan.

dat = {'N': X_train.shape[0], 'M': X_train.shape[1], 'X': X_train, 'y': y_train}

code = """
        data {
            int<lower=0> N;
            int<lower=0> M;
            matrix[N, M] X;
            int<lower=0, upper=1> y[N];
        } parameters {
            real beta0;
            vector[M] beta; 
        } model {
            for (i in 1:N)
                y[i] ~ bernoulli(inv_logit (beta0 + dot_product(X[i] , beta)));
        }
        """

Für den Stan-Code habe ich auf den Artikel [hier] verwiesen (http://aaaazzzz036.hatenablog.com/entry/2013/11/06/225417).

Stans Code ist in mehrere Blöcke unterteilt, diesmal jedoch

data
parameters
model

Es ist in drei Teile gegliedert.

Schauen wir uns jeden an.

Datenblock

data {
    int<lower=0> N;
    int<lower=0> M;
    matrix[N, M] X;
    int<lower=0, upper=1> y[N];
}

Dies ist eine sogenannte platzhalterähnliche Rolle, bei der Werte von der Python-Seite gepackt werden. Auf der Python-Seite

dat = {'N': X_train.shape[0], 'M': X_train.shape[1], 'X': X_train, 'y': y_train}

Da ich einen solchen Wert zugewiesen habe, wird er auf der Stan-Seite verwendet. diesmal,

N: Anzahl der Daten
M: Anzahl der Funktionen
X: NxM-Matrix-Trainingsdaten
y: N-dimensionale Lehrerdaten. Es gibt eine Grenze von 1 Obergrenze und 0 Untergrenze.

Jede Variable wird wie verwendet.

Parameterblock

parameters {
    real beta0;
    vector[M] beta; 
}

Der Parameterblock ist ein Block, der die Zielvariablen deklariert, die Sie schätzen möchten, z. B. die Parameter der Wahrscheinlichkeitsverteilung. Dieses Mal möchte ich die Regressionsparameter (Abschnitt $ \ beta $ und Regressionskoeffizient $ \ beta_i \ (i = 1,2, \ cdots, M) $) schätzen, also deklariere ich sie hier.

Modellblock

Das ist der Schlüssel. Erklärung eines statistischen Modells.

model {
    for (i in 1:N)
        y[i] ~ bernoulli(inv_logit (beta0 + dot_product(X[i] , beta)));
}

Es gibt mehrere Stan-Funktionen, daher werde ich sie zuerst vorstellen.

real dot_product(vector x, vector y) Stan Manual v2.9 p358 Berechnen Sie das innere Produkt. Hier entspricht es dem Teil $ \ beta_ {1} x_ {1, i} + \ cdots + \ beta_ {M} x_ {M, i} $ der Formel, und jeder Datensatz und Parameter wird multipliziert und addiert.

real inv_logit(real y) Stan Manual v2.9 p338 Diese Funktion führt die Berechnung durch, die durch die folgende Formel ausgedrückt wird. Mit anderen Worten, es ist eine Standard-Sigmoid-Funktion. $ \operatorname{inv \\_ logit}(y) = {1 \over 1 + \exp(−y)} $

y ~ bernoulli(theta) Stan Manual v2.9 p384 Dieser Ausdruck, der durch "~" dargestellt wird, heißt "Sampling Statements" und ist tatsächlich der folgende Ausdruck: increment_log_prob(bernoulli_log(y, theta)); Wird aufgerufen und ausgeführt. Da die hintere Wahrscheinlichkeit durch Bayes'sche Schätzung berechnet wird, wird "y" als beobachtete Daten und "Theta" als zu schätzender Parameter eingegeben. theta ist der Teil von inv_logit (beta0 + dot_product (X [i], beta)), der die zuvor erläuterte erklärende Variable und den zuvor erläuterten Parameter $ \ beta_i $ linear verbindet und sie als Argument für die Standard-Sigmoid-Funktion festlegt. Richtig. Zusammenfassend ist diese eine Zeile der Teil, der die logarithmische Wahrscheinlichkeit und die Wahrscheinlichkeitsfunktion der Bernoulli-Verteilung berechnet

B(y_i|\theta) = \theta^{y_i}(1-\theta)^{1-y_i}

Und der logarithmische ist

\ln(B(y_i|\theta)) = y_i \ln (\theta) + (1-y_i) \ln(1-\theta)

Da die Anzahl der Daten $ N $ beträgt, handelt es sich um eine gleichzeitige Wahrscheinlichkeit. Wenn Sie also den Logarithmus verwenden, wird er hinzugefügt.

\sum_{i=1}^N\ln(B(y_i|\theta)) = \sum_{i=1}^N\ \\{ y_i \ln (\theta) + (1-y_i) \ln(1-\theta) \\}

ist. Dieser Additionsteil ist der Teil, der von "increment_log_prob" verarbeitet wird, und er wird durch Stapeln hinzugefügt, während mit der "for" -Anweisung gedreht wird (also inkrementieren).

Außerdem wird die vorherige Verteilung nicht explizit für die Parameter "beta0, beta" festgelegt, aber in diesem Fall legt Stan die vorherige Verteilung ohne Informationen fest (gleichmäßige Verteilung des Bereichs (−∞, ∞)). ..

Aus dem obigen Stan-Code ist der Wert proportional zur posterioren Verteilung (Wahrscheinlichkeit x vorherige Verteilung) bereit. Dies wird verwendet, um Zufallszahlen von HMC und NUTS zu generieren.

Referenz: Stan Manual v2.9 　　　26.3. Sampling Statements 　　　6.2. Priors for Coefficients and Scales

Modell kompilieren

Stan übersetzt diesen Code in C ++ - Code, kompiliert ihn und führt ihn aus, um die Dinge zu beschleunigen. Es ist ineffizient, den Code jedes Mal zu kompilieren. Sobald der Code festgelegt ist, ist es praktisch, ihn zu kompilieren und als Python-Objekt zu verwenden.

# Build Stan model
%time stm = StanModel(model_code=code)

`out`


CPU times: user 1.36 s, sys: 113 ms, total: 1.48 s
Wall time: 15.4 s

Sobald das Modell erstellt ist, ist es Zeit zum Probieren. Hier dauert es am meisten. Die Einstellungen für die Anzahl der Proben sind wie folgt.

Anzahl der Proben: 3000
Anzahl der Aufwärmübungen [^ 1](Anzahl der Einbrennungen): 200
Anzahl der Ketten: 200

(3000-200) * 2 = Probenahme, bis die Anzahl der Proben 5600 erreicht.

Als Einstellung der Ausführung

Anzahl paralleler Läufe: Auf -1 setzen, um die maximal mögliche Anzahl paralleler Läufe auf der Maschine auszuführen.
Verwenden Sie NUTS für den Stichprobenalgorithmus.

[^ 1]: Laut BDA3 (Notes on Bayesian Data Analysis, 3. Ausgabe, S. 282) kann das Einbrennen irreführend sein. Es wird daher empfohlen, es als Aufwärmen zu bezeichnen.

n_itr = 3000
n_warmup = 200
chains = 2

#Probenahme durchführen
%time fit = stm.sampling(data=dat, iter=n_itr, chains=chains, n_jobs=-1, warmup=n_warmup, algorithm="NUTS", verbose=False)

#Probensäule extrahieren
la    = fit.extract(permuted=True)  # return a dictionary of arrays
#Parametername
names = fit.model_pars 
#Anzahl der Parameter
n_param = np.sum([1 if len(x) == 0 else x[0] for x in fit.par_dims])

Die Durchführung der Probenahme dauerte ca. 5 Minuten.

`out`


CPU times: user 156 ms, sys: 59.1 ms, total: 215 ms
Wall time: 5min 24s

Zeigt eine Zusammenfassung der Stichprobenergebnisse an.

print fit

Sie können den Mittelwert, den Standardfehler, die Standardabweichung, jedes Perzentil usw. für jeden Parameter anzeigen. Informationen zur Wahrscheinlichkeit werden auch unten angezeigt. Da der Wert "Rhat", der die Konvergenz des Stichprobenergebnisses beurteilt, 1,0 ist, ist ersichtlich, dass die Verteilung ausreichend konvergiert. (Wenn Rhat 1,1 oder weniger ist, kann im Allgemeinen beurteilt werden, dass es ausreichend konvergiert hat [^ 2])

[^ 2]: Referenz: Zitiernotiz über den spezifischen Wert von Rhat in der Konvergenzdiagnose von MCMC

`out`


Inference for Stan model: anon_model_b12e3f2368679a0c562b9f74618b2f82.
2 chains, each with iter=3000; warmup=200; thin=1; 
post-warmup draws per chain=2800, total post-warmup draws=5600.

          mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
beta0     5.02  8.0e-3    0.6   3.89   4.61   5.01   5.42   6.23 5600.0    1.0
beta[0]  -1.04  2.1e-3   0.15  -1.35  -1.15  -1.04  -0.93  -0.75 5600.0    1.0
beta[1]  -2.77  3.0e-3   0.23  -3.22  -2.92  -2.76  -2.61  -2.33 5600.0    1.0
beta[2]  -0.04  1.2e-4 8.9e-3  -0.06  -0.05  -0.04  -0.04  -0.03 5600.0    1.0
beta[3]  -0.41  1.6e-3   0.12  -0.66  -0.49  -0.41  -0.33  -0.18 5600.0    1.0
beta[4]  -0.08  1.8e-3   0.14  -0.35  -0.17  -0.08   0.02   0.18 5600.0    1.0
beta[5] 2.5e-3  3.4e-5 2.5e-3-2.3e-3 7.0e-4 2.4e-3 4.1e-3 7.7e-3 5600.0    1.0
lp__    -322.4    0.02   1.86 -326.8 -323.4 -322.0 -321.0 -319.7 5600.0    1.0

Samples were drawn using NUTS(diag_e) at Sun Dec 13 02:53:45 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Der erwartete Wert der Verteilung jedes Parameters wird als EAP-Schätzung (Expected A Posteriori) bezeichnet. Die EAP-Schätzung kann durch einfaches Mitteln der Probenergebnisse erhalten werden. Es ist einfach: grinsen:

#Extraktion der EAP-Schätzliste für jeden Parameter
mean_list = np.array(fit.summary()['summary'])[:,0]

Als nächstes finden Sie die MAP-Schätzung (Maximum A Posteriori). Dies wurde mit ein wenig Einfallsreichtum berechnet, indem die Kerneldichte geschätzt und der Wert eingegeben wurde, der den Maximalwert annimmt. (Wenn jemand einen einfacheren Weg kennt, lass es mich wissen, m (_ _) m)

#Listen Sie die Generierung von MAP-Schätzungen für jeden Parameter auf
ddd = la['beta0']
def find_MAP(data):
    kde = gaussian_kde(data)
    x_range = np.linspace(min(data), max(data), 501)
    eval_kde = kde.evaluate(x_range)
    #plt.plot(x_range, eval_kde)
    return x_range[np.argmax(eval_kde)]

MAP_list = []
MAP_list.append(find_MAP(ddd))
for i in range(n_param-1):
    MAP_list.append(find_MAP(la['beta'][:,i]))

Schließlich wird das Stichprobenergebnis jedes Parameters in einem Diagramm dargestellt. Die linke Seite ist das abgetastete Ergebnis der Kerneldichte, und die rechte Seite wird als Trace-Plot bezeichnet, bei dem die Werte der Zufallszahlen in der Reihenfolge der Abtastung aufgetragen werden.

f, axes = plt.subplots(n_param, 2, figsize=(15, 4*n_param))
cnt = 0
for name in names:
    dat = la[name]
    if dat.ndim == 2:
        for j in range(dat.shape[1]):
            d = dat[:,j]
            sns.distplot(d, hist=False, rug=True, ax=axes[cnt, 0])
            sns.tsplot(d,   alpha=0.8, lw=.5, ax=axes[cnt, 1])
            cnt += 1
    else:
        # Intercept
        sns.distplot(dat, hist=False, rug=True, ax=axes[cnt, 0])
        sns.tsplot(dat,   alpha=0.8, lw=.5, ax=axes[cnt, 1])
        cnt += 1

name_list = []
for name in names:
    if la[name].ndim == 2:
        for i in range(dat.shape[1]):
            name_list.append("{}{}".format(name,i+1))
    else:
        name_list.append(name)

for i in range(2):
    for j, t in enumerate(name_list):
        axes[j, i].set_title(t)
plt.show()

Ebenso habe ich eine Grafik zur Wahrscheinlichkeit geschrieben.

# Likelihood
f, axes = plt.subplots(1, 2, figsize=(15, 4))

sns.distplot(la['lp__'], hist=False, rug=True, ax=axes[0])
sns.tsplot(la['lp__'],   alpha=0.8, lw=.5, ax=axes[1])
axes[0].set_title("Histgram of likelihood.")
axes[1].set_title("Traceplot of likelihood.")
plt.show()

Lassen Sie uns abschließend die Leistung der Schätzergebnisse betrachten. Zuerst wurden die Daten durch die Holdout-Methode in 8: 2 unterteilt. Schauen wir uns also die Leistung bei der Neuzuweisung und die Generalisierungsleistung an. Vergleichen wir mit EAP- und MAP-Schätzungen als zu verwendenden Parametern.

#Ersetzen Sie den geschätzten Parameterwert und die Wahrscheinlichkeit p_Eine Funktion zur Berechnung von i.
def logis(x, beta):
    assert len(beta) == 7
    assert len(beta) == 7
    if type(x) != 'np.array':
        x = np.array(x)
    tmp = [1]
    tmp.extend(x)
    x = tmp
    return (1+np.exp(-np.dot(x, beta)))**(-1)

#Beurteilen Sie 0 oder 1 mit der eingestellten Schwelle. Der Standardschwellenwert ist 0.5
def check_accuracy(data, target, param, threshold = 0.5):
    ans_list = []
    for i in range(len(data)):
        idx = data.index[i]
        res = logis(data.ix[idx], param)
        ans = 1 if res > threshold else 0
        ans_list.append(ans == target.ix[idx])

    return np.mean(ans_list)


param = mean_list[0:7]

#Richtige Antwortrate bei Neuzuweisung
print u"[train][EAP] Accuracy:{0:.5f}".format(check_accuracy(X_train, y_train, param))
print u"[train][MAP] Accuracy:{0:.5f}".format(check_accuracy(X_train, y_train, MAP_list))

#Siehe Generalisierungsleistung durch korrekte Antwortrate der Testdaten
print "[test][EAP] Accuracy:{0:.5f}".format(check_accuracy(X_test, y_test, param))
print "[test][MAP] Accuracy:{0:.5f}".format(check_accuracy(X_test, y_test, MAP_list))

`out`


[train][EAP] Accuracy:0.79073
[train][MAP] Accuracy:0.80618
[test][EAP] Accuracy:0.81006
[test][MAP] Accuracy:0.79888

Zum Zeitpunkt der Neuzuweisung ist die korrekte Antwortrate für die MAP-Schätzung höher, aber die EAP-Schätzmethode scheint hinsichtlich der Generalisierungsleistung besser zu sein. Ich bin gespannt, ob dies passiert oder ob es eine solche Tendenz gibt.

Referenz

http://aaaazzzz036.hatenablog.com/entry/2013/11/06/225417 → Ich habe mich auf den Stan-Code bezogen.

Stan Manual (Modeling Language Benutzerhandbuch und Referenzhandbuch, v2.9.0) https://github.com/stan-dev/stan/releases/download/v2.9.0/stan-reference-2.9.0.pdf