[Statistiques] Expliquons en détail l'exécution de la régression logistique en stan (avec jeu de données Titanic)

C'est un article qui utilise Stan pour renvoyer logistiquement les données du Titanic et évaluer davantage les performances de la classification.

Le langage de programmation probabiliste "Stan" utilisé dans cet article utilise la méthode Hamiltonian Monte Carlo (méthode HMC) et NUTS pour estimer les paramètres de distribution. Strictement parlant, le principe de génération de nombres aléatoires est différent, mais une méthode légèrement plus simple est la méthode de Monte Carlo en chaîne de Markov, la méthode de Metropolis Hastings (méthode MH). J'ai écrit @ kenmatsu4 à propos de ce principe de fonctionnement

Il y a deux points, veuillez donc vous y référer si vous le souhaitez. Je pense que la méthode MH et la méthode HMC ne sont pas si différentes si l'intention est de donner une image de ce qu'elles font. (Fondamentalement, la différence d'efficacité d'échantillonnage aléatoire)

environnement

OSX Yosemite 10.10.5
Python 2.7
Anaconda 3.18.9
PyStan 2.8.0

Le code complet est publié sur GitHub.

Introduction de PyStan

J'avais peur que les choses n'allaient pas bien, mais c'était simple et cela fonctionnait.

conda install pystan

Préparation de la bibliothèque

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

Préparation des données

[Apprentissage automatique] Résumé et exécution de l'évaluation / index du modèle (avec l'ensemble de données Titanic) Utilise les données Titanic qui ont été converties en valeurs numériques. Faire. Nous l'avons préparé sur GitHub, veuillez donc l'utiliser également.

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

Divisez les données en la variable explicative «data» et la variable objective «target».

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

#Données d'entraînement(80%),données de test(20%)Diviser en
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,)]

Autrement dit, les données d'entraînement $ X $ sont

C'est une matrice avec le nombre de données $ N = 712 $ et le nombre de caractéristiques $ M = 6 $. Le vecteur horizontal $ x_i $ est les données pour une personne. Les six quantités de caractéristiques sont pclass 　　sex　　 age　　 sibsp　　 parch　　 fare est.

Les données de l'enseignant sont les informations de survie du Titanic. Survie = 1 et la personne qui n'est jamais revenue = 0 $ N $ dimension vector

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

est.

À propos de la régression logistique

Maintenant, quand il existe de telles données, le [modèle de régression logistique (lien wikipedeia)](https://ja.wikipedia.org/wiki/Logistic regression) est exprimé comme suit.

\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 $ peut être considéré comme la probabilité de survie des $ i $ èmes données d'enregistrement. Réorganisons un peu la formule pour voir pourquoi elle peut être considérée comme une probabilité de survie. Tout d'abord, mettez les deux côtés dans $ \ exp (\ cdot) $

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

Résoudre cela pour $ p_i $

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

Voici la fonction sigmoïde standard $ \sigma(x) = {1 \over 1+ e^{-x}} $

Utiliser. De forme similaire à la formule précédente, en remplaçant $ x $ par $ (\ beta + \ beta_ {1} x_ {1, i} + \ cdots + \ beta_ {M} x_ {M, i}) $ J'ai fait

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

Vous pouvez le traiter comme.

Cela ressemble à ceci dans un graphique.

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)

Comme vous pouvez le voir dans le graphique, la fonction sigmoïde $ \ sigma (x) $ prend une valeur dans la plage (0, 1), donc cette valeur peut être interprétée comme une probabilité et utilisée.

De ce qui précède, la partie qui ressemble à une régression linéaire normale,

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

En utilisant la fonction sigmoïde bien basée sur

\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)

On peut dire que la régression logistique doit prendre la forme et interpréter la valeur comme une probabilité.

Stan est introduit en essayant d'estimer les paramètres $ \ beta, \ beta_1, \ cdots, \ beta_M $ de cette expression en utilisant cette fois Stan.

Courir sur Stan

Créez un objet Dictionary contenant des données d'entraînement à transmettre à 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)));
        }
        """

Pour le code stan, je me suis référé à l'article ici.

Le code de Stan est divisé en plusieurs blocs, mais cette fois

data
parameters
model

Il est divisé en trois parties.

Regardons chacun d'eux.

bloc de données

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

Il s'agit d'un rôle de type espace réservé où les valeurs sont regroupées du côté Python. Du côté Python

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

Puisque j'attribuais une valeur comme celle-ci, elle est utilisée du côté Stan. cette fois,

N: nombre de données
M: nombre de fonctionnalités
X: données d'entraînement de la matrice NxM
y: données de l'enseignant à N dimensions. Il y a une limite de 1 limite supérieure et 0 limite inférieure.

Chaque variable est utilisée comme.

bloc de paramètres

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

Le bloc de paramètres est un bloc qui déclare les variables cibles que vous souhaitez estimer, telles que les paramètres de la distribution de probabilité. Cette fois, je veux estimer les paramètres de régression (section $ \ beta $ et coefficient de régression $ \ beta_i \ (i = 1,2, \ cdots, M) $), donc je les déclare ici.

bloc modèle

C'est la clé. Déclaration d'un modèle statistique.

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

Il existe plusieurs fonctions Stan, je vais donc les présenter en premier.

real dot_product(vector x, vector y) Stan Manual v2.9 p358 Calculez le produit intérieur. Ici, cela correspond à la partie $ \ beta_ {1} x_ {1, i} + \ cdots + \ beta_ {M} x_ {M, i} $ de la formule, et chaque enregistrement et paramètre sont multipliés et ajoutés.

real inv_logit(real y) Stan Manual v2.9 p338 C'est une fonction qui effectue le calcul exprimé par la formule suivante. En d'autres termes, c'est une fonction sigmoïde standard. $ \operatorname{inv \\_ logit}(y) = {1 \over 1 + \exp(−y)} $

y ~ bernoulli(theta) Stan Manual v2.9 p384 Cette expression représentée par ~ est appelée Instructions d'échantillonnage, et c'est en fait l'expression suivante increment_log_prob(bernoulli_log(y, theta)); Est appelé et exécuté. Puisque la probabilité a posteriori par estimation bayésienne est calculée, «y» est entré comme données observées et «thêta» est entré comme paramètre à estimer. theta est la partie de ʻinv_logit (beta0 + dot_product (X [i], beta))` qui joint linéairement la variable explicative et le paramètre $ \ beta_i $ expliqués précédemment et les définit comme l'argument de la fonction sigmoïde standard. Droite. En résumé, cette ligne est la partie qui calcule la vraisemblance logarithmique et la fonction de probabilité de la distribution de bernoulli

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

Et le logarithmique est

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

De plus, puisqu'il s'agit d'une probabilité simultanée du nombre de données $ N $, si vous prenez le logarithme, il sera ajouté.

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

est. Cette partie additionnelle est la partie traitée par ʻincrement_log_prob, et elle est ajoutée par empilement en tournant avec l'instruction for` (donc incrément).

De plus, la distribution antérieure n'est pas explicitement définie pour les paramètres «beta0, beta», mais dans ce cas, Stan définit la distribution a priori non informative (distribution uniforme de l'intervalle (−∞, ∞)). ..

À partir du code Stan ci-dessus, la valeur proportionnelle à la distribution postérieure (probabilité x distribution antérieure) est prête. Ceci est utilisé pour générer des nombres aléatoires par HMC et NUTS.

Référence: Stan Manual v2.9 　　　26.3. Sampling Statements 　　　6.2. Priors for Coefficients and Scales

Compiler le modèle

Stan traduit ce code en code C ++, le compile et l'exécute pour accélérer les choses. Il est inefficace de compiler le code à chaque fois, donc une fois que le code est décidé, il est pratique de le compiler et de l'avoir en tant qu'objet Python.

# 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

Une fois le modèle construit, il est temps d'échantillonner. C'est là que cela prend le plus de temps. Les paramètres relatifs au nombre d'échantillonnages sont les suivants.

Nombre d'échantillons: 3000
Nombre d'échauffements [^ 1](nombre de burn-ins): 200
Nombre de chaînes: 200

(3000-200) * 2 = Échantillonnage jusqu'à ce que le nombre d'échantillons atteigne 5600.

Comme cadre d'exécution

Nombre de courses parallèles: réglez -1 pour exécuter au nombre maximum de courses parallèles possible sur la machine.
Utilisez NUTS pour l'algorithme d'échantillonnage.

[^ 1]: Selon BDA3 (Notes on Bayesian Data Analysis 3e édition p282), le burn-in peut être trompeur, il semble donc recommandé de l'appeler échauffement.

n_itr = 3000
n_warmup = 200
chains = 2

#Effectuer un échantillonnage
%time fit = stm.sampling(data=dat, iter=n_itr, chains=chains, n_jobs=-1, warmup=n_warmup, algorithm="NUTS", verbose=False)

#Extraire la colonne d'échantillon
la    = fit.extract(permuted=True)  # return a dictionary of arrays
#Le nom du paramètre
names = fit.model_pars 
#Nombre de paramètres
n_param = np.sum([1 if len(x) == 0 else x[0] for x in fit.par_dims])

Il a fallu environ 5 minutes pour exécuter l'échantillonnage.

`out`


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

Affiche un résumé des résultats d'échantillonnage.

print fit

Vous pouvez voir la moyenne, l'erreur standard, l'écart type, chaque centile, etc. pour chaque paramètre. Des informations sur la probabilité sont également affichées en bas. Puisque la valeur «Rhat» qui juge la convergence du résultat de l'échantillonnage est 1,0, on peut voir que la distribution est suffisamment convergée. (En général, si Rhat est égal ou inférieur à 1,1, on peut juger qu'il a suffisamment convergé [^ 2])

[^ 2]: Référence: Citation mémo sur la valeur spécifique de Rhat dans le diagnostic de convergence de 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).

La valeur attendue de la distribution de chaque paramètre est appelée estimation EAP (Expected A Posteriori). L'estimation du PAE peut être obtenue en faisant simplement la moyenne des résultats de l'échantillon. C'est facile: sourire:

#Extraction de la liste d'estimations EAP pour chaque paramètre
mean_list = np.array(fit.summary()['summary'])[:,0]

Ensuite, trouvez l'estimation MAP (Maximum A Posteriori). Ceci a été calculé en estimant la densité du noyau et en introduisant la valeur qui prend la valeur maximale. (Si quelqu'un connaît un moyen plus simple, merci de me le faire savoir m (_ _) m)

#Génération de listes d'estimations MAP pour chaque paramètre
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]))

Enfin, le résultat d'échantillonnage de chaque paramètre est visualisé dans un graphique. Le côté gauche est le résultat échantillonné en densité du noyau, et le côté droit est appelé trace plot, qui trace les valeurs des nombres aléatoires dans l'ordre d'échantillonnage.

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()

De même, j'ai écrit un graphique concernant la probabilité.

# 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()

Enfin, regardons les performances des résultats de l'estimation. Au début, les données étaient divisées en 8: 2 par la méthode d'exclusion, voyons donc les performances lors de la réaffectation et les performances de généralisation. Comparons en utilisant les estimations EAP et MAP comme paramètres à utiliser.

#Remplacez la valeur estimée du paramètre et la probabilité p_Une fonction pour calculer 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)

#Jugez 0 ou 1 avec le seuil défini. Le seuil par défaut est 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]

#Taux de réponse correct lors de la réaffectation
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))

#Voir les performances de généralisation par taux de réponse correct des données de test
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

Au moment de la réaffectation, le taux de réponse correcte est plus élevé pour l'estimation MAP, mais la méthode d'estimation du PAE semble être meilleure en termes de performance de généralisation. Je suis curieux de savoir si cela se produit ou s'il existe une telle tendance.

référence

http://aaaazzzz036.hatenablog.com/entry/2013/11/06/225417 → J'ai fait référence au code Stan.

Manuel Stan (Guide de l’utilisateur et manuel de référence du langage de modélisation, v2.9.0) https://github.com/stan-dev/stan/releases/download/v2.9.0/stan-reference-2.9.0.pdf