Aperçu

Participation à la session de lecture "Introduction à la modélisation statistique pour l'analyse des données" (dite Midorimoto) organisée par l'entreprise, C'était un livre très facile à comprendre, mais malheureusement pour les utilisateurs de Mac, l'exemple d'implémentation est WinBUGS Parce que c'était L'approche d'estimation bayésienne du modèle linéaire généralisé du chapitre 9 a été implémentée dans Python + STAN.

Ce que j'ai fait

Suivez à peu près les étapes ci-dessous.

Générer des données factices à partir de la distribution de probabilité en fonction de paramètres spécifiques
Définissez le modèle de prédiction
À partir des données factices et du modèle de prédiction, estimer les paramètres (distribution postérieure) qui ont généré les données avec MCMC et faire correspondre les réponses.

Plus précisément, la taille corporelle d'une certaine plante (en prenant des valeurs discrètes par incréments de 0,1 de 3,0 à 7,0) est utilisée comme variable explicative. Estimer la distribution de probabilité du nombre de graines (entier 0 ou plus) qui suit la distribution de Poisson.

Outils et bibliothèques utilisés

Échantillonneur MCMC: STAN (PyStan)
Calcul matriciel, opération de tableau: NumPy
Dessin graphique: matplotlib

Méthode d'estimation

MCMC La méthode MCMC est une abréviation de la méthode Marcov Chain Monte Carlo. En japonais, on l'appelle la méthode de Monte Carlo en chaîne de Markov. Chaîne de Markov La méthode de Monte Carlo est une chaîne de Markov qui utilise la méthode de Monte Carlo. ... Je suis devenu une totologie, alors je vais vous expliquer un peu plus.

Chaîne de Markov

La propriété selon laquelle un certain état ne dépend que de l'état immédiatement précédent est appelée ** propriété Markov **. ** Chaîne de Markov ** est un modèle probabiliste dans lequel des états qui ne dépendent que de l'état immédiatement précédent apparaissent dans une chaîne. Du point de vue d'un ingénieur, c'est plus facile à comprendre si vous le considérez comme un ** type d'automate **.

Méthode de Monte Carlo

La méthode de Monte Carlo est un nom général pour les calculs numériques et les simulations utilisant des nombres aléatoires.

En d'autres termes, la ** méthode de Monte Carlo par chaîne de Markov ** est une méthode pour générer une chaîne de Markov à l'aide de nombres aléatoires. Ici, en particulier, il fait référence à un algorithme qui génère une distribution de probabilité (distribution postérieure des paramètres) en utilisant les propriétés de la chaîne de Markov.

Estimation bayésienne x MCMC

En combinant le cadre d'estimation bayésien (distribution a posteriori ∝ vraisemblance x distribution a priori) et MCMC qui échantillonne la distribution de probabilité proportionnelle à la vraisemblance, Même si le modèle ne peut pas être résolu analytiquement, la distribution postérieure peut être estimée tant qu'elle peut être exprimée par une formule mathématique.

la mise en oeuvre

Une explication détaillée du modèle linéaire généralisé et du MCMC semble être très longue, nous allons donc commencer à l'implémenter. "Introduction à la modélisation statistique pour l'analyse des données - modèle linéaire généralisé, modèle hiérarchique de Bayes, MCMC (science des probabilités et de l'information)" est facile à comprendre.

Génération de données d'entraînement

La taille du corps est de 3,0 ~ 7,0, Moyenne μ = 1,5 + 0,1 * pour chaque individu Elle est générée à partir de la distribution de Poisson de la taille corporelle.

def generate(n):
	for i in range(n):
		x = round(random.random() * 4 + 3, 1) # 3.0 ~ 7.Nombre aléatoire jusqu'à 0
		mu = math.exp(1.5 + 0.1*x)
		print (x, np.random.poisson(mu))

"Taille corporelle" et "nombre de graines".

STAN mcmc.stan

data {
  int<lower=1> N;
  vector[N] x;
  int<lower=0> y[N];
}
parameters {
  real beta1;
  real beta2;
}
model {
  for (i in 1:N) {
    y[i] ~ poisson(exp(beta1 + beta2 * x[i])); //Distribution de Poisson x fonction de lien logarithmique
  }
  beta1 ~ normal(0, 1000); //Moyenne 0,Distribution normale de la variance 1000 ≒ aucune information distribution antérieure
  beta2 ~ normal(0, 1000); //Moyenne 0,Distribution normale de la variance 1000 ≒ aucune information distribution antérieure
}

Comment lire

data Ce sont les données à transmettre au programme stan. Déclarez au format ** {type de données} {nom de la variable}; **. Transmettez les données de Python en spécifiant le nom de la variable écrit ici.

parameters Il s'agit de la variable utilisée dans le modèle décrit dans stan. Cette fois, la section β1 et le coefficient β2 de la fonction de lien logarithmique de la distribution de Poisson sont fixés à l'intérieur du STAN. Elle est générée à partir d'une distribution a priori non informative (une distribution normale avec une très grande variance approximative).

model C'est un modèle de prédiction. Les opérateurs ** '~' ** signifient que le terme de gauche suit la distribution de probabilité du terme de droite. Ici, le nombre de graines ** y ** suit une distribution de Poisson avec ** exp (β1 + β2x) ** comme fonction de lien (c'est-à-dire fonction de lien logarithmique).

Python

import numpy as np
import pystan
import matplotlib.pyplot as plt

data = np.loadtxt('./data.txt', delimiter=' ')

#Génération de données à passer à l'interface Stan
x = data[:,0] #notation numpy, découpez la première colonne de données
y = data[:,1].astype(np.int) #Puisqu'il s'agit de la variable objective de la distribution de Poisson, elle est convertie en une valeur entière.
N = data.shape[0] #Le nombre de données

stan_data = {'N': N, 'x': x, 'y': y}

fit = pystan.stan(file='./mcmc.stan',\
	data=stan_data, iter=5000, chains=3, thin=10)
    # iter =Numéro de chaque échantillon
    # chain =Répéter n ensembles d'échantillons spécifiés par iter
    # thin =Numéro d'éclaircissage de l'échantillon

Résultat d'exécution

Si cela se passe bien, un journal comme celui-ci sortira. Puisque STAN lui-même est implémenté en C ++, la compilation est en cours d'exécution.

INFO:pystan:COMPILING THE C++ CODE FOR MODEL ~ NOW.
Chain 1, Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 0, Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2, Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1, Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 2, Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 0, Iteration: 1000 / 10000 [ 10%]  (Warmup)
...
Chain 0, Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2, Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1, Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 0, Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 2, Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1, Iteration: 10000 / 10000 [100%]  (Sampling)# 
...
#  Elapsed Time: 9.51488 seconds (Warm-up)
#                10.7133 seconds (Sampling)
#                20.2282 seconds (Total)
#

Graphique

Statistiques

summary

Inference for Stan model: anon_model_381b30a63720cfb3906aa9ce3e051d13.
3 chains, each with iter=10000; warmup=5000; thin=10; 
post-warmup draws per chain=500, total post-warmup draws=1500.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
beta1   1.41  2.4e-3   0.05   1.31   1.38   1.41   1.45   1.52    481    1.0
beta2   0.12  4.6e-4   0.01    0.1   0.11   0.12   0.13   0.14    478    1.0
lp__  7821.4    0.04   0.97 7818.8 7821.1 7821.7 7822.1 7822.4    496    1.0

Samples were drawn using NUTS(diag_e) at Tue Feb  9 23:31:02 2016.
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).

Lorsque la méthode fit.summary () est exécutée, la sortie sera la suivante. Premièrement, en regardant β1, cela montre que la moyenne des échantillons est de 1,41 et qu'il y a 95% de chances qu'elle soit comprise entre 1,31 et 1,52. (Cela semble être appelé une section de crédit en termes bayésiens.)

Puisqu'il n'y a qu'une seule montagne dans la distribution cette fois, je vais regarder la valeur représentative en moyenne, β1 (1,5) est 1,41 et β2 (0,1) est 0,12, qui ne sont pas significativement différents, et les valeurs numériques peuvent être estimées.

Rhat Lors de l'appel du code de Stan depuis Python, j'ai répété l'échantillonnage 3 fois avec le paramètre * chain *. Dans l'estimation de la distribution postérieure des paramètres par MCMC, les valeurs initiales des paramètres sont déterminées de manière appropriée. Selon le modèle, la distribution de probabilité estimée peut différer en fonction de la valeur initiale. Répétez l'échantillonnage plusieurs fois et vérifiez si la distribution de probabilité converge vers ** Rhat **, qui quantifie la variation entre les ensembles. Il semble que ce soit OK si c'est 1.1 ou moins, mais cette fois, on peut juger qu'il n'y a pas de problème car les deux bêta1 et bêta2 sont inférieurs à cela.

Livre de référence

Comme l'accent est mis sur l'introduction de l'implémentation, la signification de chaque argument tel que * thin * lors de l'appel de STAN, et Pourquoi la première moitié de l'échantillonnage est-elle utilisée pour le * réchauffement *? En outre, MCMC est un terme général pour les méthodes en premier lieu, et quel est l'algorithme spécifique? Si vous souhaitez en savoir plus, veuillez lire ce qui suit.

[Introduction à la modélisation statistique pour l'analyse des données - modèle linéaire généralisé, modèle hiérarchique de Bayes, MCMC (science des probabilités et de l'information)](http://www.amazon.co.jp/%E3%83%87%E3 % 83% BC% E3% 82% BF% E8% A7% A3% E6% 9E% 90% E3% 81% AE% E3% 81% 9F% E3% 82% 81% E3% 81% AE% E7% B5 % B1% E8% A8% 88% E3% 83% A2% E3% 83% 87% E3% 83% AA% E3% 83% B3% E3% 82% B0% E5% 85% A5% E9% 96% 80 % E2% 80% 95% E2% 80% 95% E4% B8% 80% E8% 88% AC% E5% 8C% 96% E7% B7% 9A% E5% BD% A2% E3% 83% A2% E3 % 83% 87% E3% 83% AB% E3% 83% BB% E9% 9A% 8E% E5% B1% A4% E3% 83% 99% E3% 82% A4% E3% 82% BA% E3% 83 % A2% E3% 83% 87% E3% 83% AB% E3% 83% BBMCMC-% E7% A2% BA% E7% 8E% 87% E3% 81% A8% E6% 83% 85% E5% A0% B1% E3% 81% AE% E7% A7% 91% E5% AD% A6-% E4% B9% 85% E4% BF% 9D-% E6% 8B% 93% E5% BC% A5 / dp / 400006973X)

[PYTHON] La première méthode de Monte Carlo en chaîne de Markov par PyStan