introduction

À la suite de la collecte de diverses données d'observation, il peut y avoir un désir d'appliquer la distribution de probabilité derrière la valeur d'observation au modèle et de l'estimer. C'est pratique pour calculer diverses statistiques et calculs de simulation. Il existe une distribution normale, une distribution exponentielle, une distribution binomiale, une distribution bêta, etc. comme modèles de distribution de probabilité, mais cette fois-ci, en implémentant la méthode MCMC (méthode Markov Chain Monte Carlo) "simplement" à partir du principe de distribution gamma. , Je voudrais essayer combien il peut être adapté. En tant que bon outil pour exécuter MCMC, WinBUGS et PyMC3 Il existe //docs.pymc.io/) etc., vous n'avez donc pas à le faire vous-même **. Cependant, il y a beaucoup de choses que vous pouvez voir en le fabriquant vous-même, et je pense qu'il y a un mérite à promouvoir la compréhension.

Qu'est-ce que la distribution gamma?

Qu'est-ce que Gamma Distribution? C'est une distribution caractérisée par deux paramètres: le paramètre de forme $ k $ et le paramètre d'échelle $ q $.

P(x) = \frac{1}{\Gamma(k)q^k}x^{k-1}e^{-x/q},\ x \geq 0 \ (k>0,\ q>0)

C'est une distribution qui prend des nombres réels non négatifs, et c'est un excellent gars qui peut également exprimer une distribution exponentielle, une distribution du chi carré ou une distribution d'Alan. Surtout, j'aime le fait qu'il puisse exprimer une distribution déformée avec un large ourlet. Il existe de nombreuses données d'observation dans le monde qui sont représentées par des valeurs non négatives telles que le nombre, la longueur, le poids, l'intervalle de temps, la population, le montant d'argent, etc., et cela semble utile pour exprimer la distribution de ces variables. En fait, il semble être utilisé en ingénierie et en économie, comme la distribution des temps d'attente et la distribution des revenus.

Qu'est-ce que MCMC?

La méthode MCMC (méthode Markov Chain Monte Carlo) est une méthode qui génère efficacement un échantillon d'une certaine distribution de probabilité et estime la nature de la distribution de probabilité à partir de l'échantillon. La chaîne de Markov est responsable de la génération de l'échantillon selon la distribution de probabilité, et Monte Carlo est responsable de l'estimation à partir de l'échantillon. Pour plus de détails, la Vidéo du commentaire du Dr Yukito Iba de l'Institut de mathématiques statistiques est très facile à comprendre. Dans cet article, nous utiliserons MCMC pour estimer la distribution gamma, nous allons donc définir les termes suivants:

$ \ Theta $: paramètre de dimension $ N $ que vous souhaitez estimer
$ P (x | \ theta) $: Modèle de distribution de probabilités calculé selon le paramètre $ \ theta $ (je le connais car c'est un modèle)
$ X_1, X_2, ..., X_n $: valeur d'échantillon observée. Cependant, l'indépendance est supposée.
$ D_n = (X_1, X_2, ..., X_n) $: Données observées (ensemble de valeurs)
$ P (\ theta | D_n) $: Distribution de probabilité du paramètre $ \ theta $ estimée à partir des données $ D_n $

Ce que je veux faire ici est de collecter beaucoup de données $ D_n $, et si la distribution de probabilité $ P (\ theta | D_n) $ du paramètre $ \ theta $ peut être bien calculée, le (du paramètre $ \ theta $ du modèle supposé Cela signifie que vous pouvez calculer la valeur attendue (basée sur la probabilité ex post facto). La distribution de probabilité $ P (\ theta | D_n) $ est exprimée par l'équation suivante du théorème de Bayes (ou de la condition de normalisation de la distribution de probabilité).

P(\theta|D_n) = P(D_n | \theta)P(\theta) / \int_{\theta} d\theta P(D_n | \theta)P(\theta)

Ici, l'intégration doit être prise pour tout $ \ theta $ avec une probabilité non nulle de produire les données $ D_n $. L'un des moyens les plus simples d'estimer cette distribution est le calcul de la grille. En d'autres termes, si $ \ theta $ est uniformément réparti sur une grille appropriée et représenté par la valeur discrète $ \ theta_i (i = 1,2, ..., I) $,

P(\theta_i|D_n) \simeq P(D_n | \theta_i)P(\theta_i) / \sum_{j \in I}  P(D_n | \theta_j)P(\theta_j)

L'intégrale du dénominateur peut être approximée par la somme, comme dans. Cela ressemble à ceci lorsque je l'écris dans un diagramme.

Cependant, pour améliorer la précision d'estimation de $ \ theta $, le nombre de valeurs discrètes $ I $ doit être suffisamment grand. Si vous essayez simplement de représenter un paramètre de dimension $ N $ avec des mailles $ M $, vous avez besoin de valeurs discrètes $ I = M ^ N $, qui exploseront à mesure que la dimension augmentera.

Ainsi, ** MCMC change la façon de penser et ne fixe pas la valeur discrète $ \ theta_i $, mais la déplace régulièrement **. Déplacer $ k (k = 1,2, ...) $, $ k $ th $ \ theta_i (i = 1, .., I) $ vers $ \ theta_i (k) $, $ \ theta_i ( Soit $ P_k (\ theta) $ la distribution de probabilité que k) $ suit. Aussi, laissez $ \ pi (x \ rightarrow x ') $ passer de $ \ theta_i (k) $ à $ \ theta_i (k + 1) $. Ensuite, il peut être illustré comme suit. P_k_k+1_map1.jpg

J'ai conçu cette façon de bien déplacer $ \ pi (x \ rightarrow x ') $, et après l'avoir déplacée plusieurs fois,

P_k(\theta) \rightarrow P(\theta | D_n) \  ({\rm as \ } k \rightarrow \infty)

Si c'est vrai, vous pouvez atteindre votre objectif. Ensuite, regardons l'algorithme Metropolis comme l'un des moyens de le déplacer.

Algorithme de Metropolis

L'algorithme Metropolis est un algorithme d'échantillonnage à partir d'une distribution de probabilité. Dans la figure précédente, nous augmenterons le nombre de caractères.

$ P _ * (\ theta) : Distribution selon le modèle que vous souhaitez estimer ( P (\ theta | D_n) $ dans cet article)
$ P_ {o_k} (\ theta) $: Distribution selon autre que $ P_ * (\ theta) $ dans $ P_k (\ theta) $

En d'autres termes, disons que $ \ theta_i $ est généré en synthétisant les distributions qui suivent et ne suivent pas le modèle. À ce stade, supposons que le mouvement par mappage de $ \ pi (x \ rightarrow x ') $ est limité comme indiqué dans la figure ci-dessous.

P_k_k+1_map2.jpg

En d'autres termes, la restriction est que $ \ theta_i $ qui a suivi ** $ P_ * (\ theta) $ se déplacera toujours vers $ P_ * (\ theta) $ ** ((1) dans la figure). De plus, supposons que $ \ theta_i $ qui a suivi ** $ P_ {o_k} (\ theta) $ peut toujours être déplacé vers $ P_ * (\ theta) $ ** (② sur la figure). Alors, si ② existe modérément, à partir de la limitation de ①, on peut voir que lorsque $ k $ augmente, il devient $ P_k (\ theta) \ rightarrow P_ * (\ theta) $. Strictement parlant, (1) s'appelle stabilité, (2) s'appelle contractivité, et il semble que la convergence de la distribution de probabilité à partir de ces deux conditions s'appelle ergodisme. (Pour plus de détails, voir Vidéo du commentaire du Dr Yukito Iba).

L'un des mappages $ \ pi (x \ rightarrow x ') $ qui satisfont ces propriétés est l'algorithme Metropolis, qui est représenté par le flux suivant.

Générez $ x '$ à partir de la distribution de probabilité symétrique $ Q (x, x') (= Q (x ', x)) $ avec le point courant $ x $.
Génère un nombre aléatoire uniforme $ 0 \ leq r <1 $.
Si $ r <{P_ * (x ')} / {P_ * (x)} $, déplacez-vous vers $ x \ rightarrow x' . Sinon, ne bougez pas ( x \ rightarrow x $).
Répétez les étapes 1 à 3.

Tant que la distribution de probabilité de 1 est symétrique, par exemple, une distribution normale peut être utilisée.

Q(x,x') = \frac{1}{\sqrt{2\pi \sigma^2}}exp\left(-\frac{(x-x')^2}{2\sigma^2}\right)

D'ailleurs, l'algorithme étendu au cas non symétrique ($ Q (x, x ') \ neq Q (x', x) $) est appelé l'algorithme de Metropolis-Hastings. En regardant l'étape 3,

$ P_ * (x ') \ geq P_ * (x) $ se déplacera toujours de x à x'. En d'autres termes, nous passerons à celui avec la probabilité la plus élevée. (Équivalent à dans la figure ci-dessus)
$ P_ * (x ') <P_ * (x) $ se déplace vers x'ou x, selon le rapport de probabilité $ P_ * (x') $ à $ P_ * (x) $. (Équivalent à dans la figure ci-dessus)

Essayez de calculer avec Python

Maintenant, estimons la distribution gamma lors de l'implémentation de MCMC en Python.

1. Préparation

Définissez la bibliothèque.

import numpy as np
import matplotlib.pyplot as plt
from math import gamma

Définit une fonction gamma, une génération d'histogramme et une fonction qui affiche la forme de la distribution de probabilité.

def gamma_pdf(x, k, q):
    p = 1./(gamma(k) * pow(q, k)) * pow(x, k-1) * np.exp( -x/q )
    return p

def genPdfHist(f, xmin, xmax, N):
    x = np.linspace(xmin,xmax, N)
    p = [f(v) for v in x]
    return x,p

def showPdf(f, xmin, xmax, N, yscale=1.0):
    x, p = genPdfHist(f, xmin, xmax, N)
    y = [e*yscale for e in p]
    plt.plot(x,y)
    plt.show()

2. Préparez les données d'observation

Ensuite, créons des données d'observation en utilisant la distribution gamma. Définissez une fonction à échantillonner pour créer des données d'observation.

def PlainSampler(pdf, xmin, xmax, N):
    n = 0
    xs = []
    while n < N:
        x = xmin + np.random.rand()*(xmax-xmin)
        p = pdf(x)
        if np.random.rand() < p:
            # accept
            xs.append(x)
            n = n + 1
    return xs

Utilisez la fonction ci-dessus pour générer des données d'observation. (Signification des données d'observation expérimentales) Tout d'abord, générons le modèle correct avec les paramètres $ (k ^ *, q ^ *) = (2,4) $.

k_true = 2.0
q_true = 4.0
gp1 = lambda x : gamma_pdf(x, k_true, q_true)

smp1 = PlainSampler(gp1, 0, 50, 1000)
h1 = plt.hist(smp1, bins=50)
showPdf(gp1, 0, 50, 100, h1[0].sum()/(h1[1][1]-h1[1][0]))
plt.show()

La ligne de la fonction showPdf est multipliée par le facteur d'ajustement pour aligner l'axe Y de l'histogramme avec la distribution de probabilité. Les données smp1 sont la valeur observée et son histogramme ressemble à ceci. Bien entendu, comme il est généré avec des nombres aléatoires, la forme de l'histogramme change à chaque fois qu'il est exécuté. pdf_sample_k2.0_q4.0.png

3. Préparez une fonction pour MCMC

Nous définirons les fonctions utilisées pour calculer le MCMC. Premièrement, la distribution de probabilité postérieureP(D_n |\theta)Est une fonction qui calcule. Probabilité simultanée en sélectionnant en outre 50 valeurs d'échantillon à partir de l'échantillon d'origine pour éviter que la valeur de probabilité calculée ne devienne trop petite et pour garantir l'indépendance des donnéesP(D_n |\theta)=P(X_1 | \theta)P(X_2 | \theta) \cdots P(X_n | \theta)Est en cours de calcul.

def pdf_gamma_sample(x, sample):
    k = x[0]
    q = x[1]
    #Empêcher P de devenir trop petit
    sample1 = np.random.choice(sample, 50)
    p = [gamma_pdf(s, k, q) for s in sample1]
    return np.prod(p)

Une fonction qui donne une distribution a priori (distribution initiale) du paramètre $ \ theta = (k, q) $. Puisqu'il n'y a aucune information préalable, générons-la à partir d'un nombre aléatoire uniforme. Cependant, je ne donnerai que la valeur maximale.

def genInit( cdim, xdim, xmax):
    x = np.random.rand(cdim, xdim)*xmax
    return x

Vient ensuite la partie principale. Une fonction qui exécute MCMC avec un mapping $ \ pi (x \ rightarrow x ') $ par l'algorithme Metropolis.

def MHmap(pdf, x, xdim, sig):
    #Échantillonnage de Gibbs (changer une seule variable)
    xnew = x + np.eye(1, xdim, np.random.randint(xdim))[0,:] * np.random.normal(loc=0, scale=sig)
    #Limiter la plage de valeurs de x
    xnew = np.maximum(xnew, 0.01)
    #Sélectionné selon les critères Metropolis
    r = np.random.rand()
    if r < pdf(xnew)/pdf(x):
        x = xnew
    return x

def MCMC_gamma_01(sample, x0, sig, N):
    #
    pdf1 = lambda x: pdf_gamma_sample(x, sample)
    #
    x = x0
    xs = []
    for i in range(N):
        #Cartographie des points selon les normes Metropolis
        x = [MHmap(pdf1, e, 2, sig) for e in x]
        #Rééchantillonnage aléatoire
        s = np.random.randint(0,len(x),(len(x)))
        x = [x[j] for j in s]
        #
        xs.append(x)
    return np.array(xs)

C'est le point lors du codage.

Utiliser l'échantillonnage Gibbs (peut-être pas nécessaire)
La fonction de transition $ Q (x, x ') $ utilise une distribution normale avec un écart type $ \ sigma $
Limiter la plage de la variable x (ne pas sauter au domaine négatif) Après la cartographie avec + π, rééchantillonner davantage (en théorie, cela ne devrait pas être nécessaire, mais je veux réduire la possibilité que les points qui volent trop loin ne reviennent pas facilement).

** Je ne suis pas sûr que cette implémentation soit optimale **, mais je peux imaginer que les bibliothèques réelles (PyMC3, WinBUGS, etc.) auraient été conçues de différentes manières. C'est là que le fossé entre la théorie et la mise en œuvre devient visible.

4. Essayez MCMC

Maintenant, nous allons MCMC avec la fonction créée en 3 pour estimer le paramètre $ \ theta = (k, q) $ du modèle de distribution gamma à partir des données d'observation smp1 créées en 2. Définissez la taille du groupe de points sur 500, la dimension du paramètre sur 2 dimensions et la valeur maximale sur 10. Soit également l'écart type de la fonction de transition $ Q (x, x ') $ 3. Effectuez le calcul MCMC 100 fois.

p0 = genInit( 500, 2, 10)
sig = 3
xs = MCMC_gamma_01(smp1, p0, sig, 100)

Après quelques minutes, le résultat sortira, alors faisons un graphique.

fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.hist(xs[-1,:,0], bins=50)
ax2.hist(xs[-1,:,1], bins=50)
ax1.set_xlabel("k")
ax2.set_xlabel("q")
k_est = xs[-1,:,0].mean()
q_est = xs[-1,:,1].mean()
print("k*:{}, q*:{}".format(k_est, q_est))
plt.show()

Le résultat est présenté ci-dessous. Ce n'est pas un histogramme très propre, mais lorsqu'il est calculé avec la valeur moyenne, il est $ (\ bar {k}, \ bar {q}) = (2,21, 4,46) $. MCMC_sample_k2.0_q4.0.png

Une comparaison de la distribution générée d'origine et de la distribution utilisant $ (\ bar {k}, \ bar {q}) = (2.21, 4.46) $ est montrée dans la figure ci-dessous. C'est un peu décalé, mais il semble que la distribution à longue queue puisse être ajustée dans une certaine mesure. MCMC_match_k2.0_q4.0.png

5. Autres résultats

Si vous essayez de changer le paramètre $ \ theta = (k, q) $ de différentes manières, il semble convenir. L'impression est qu'elle convergera après l'avoir répétée environ 100 fois.

Considération

De ce qui précède, les tendances suivantes ont été observées à partir des résultats de l'estimation de la distribution gamma en mettant simplement en œuvre MCMC.

Il est possible (dans une certaine mesure) d'estimer les paramètres de la distribution gamma originale à partir de la moyenne de l'échantillon en utilisant les données d'observation générées à partir de la distribution gamma.
Cependant, ** il y avait une tendance à dépendre de la distribution initiale comme la distribution antérieure donnée à MCMC et la fonction de transition $ Q (x, x ') $ (dispersion) **.
Cela semble utile car il peut être utilisé pour une fonction à forte non-linéarité appelée distribution gamma, et il peut être utilisé en principe même si la dimension du paramètre est élevée.

En outre ...

La partie qui calcule le rapport de probabilité à partir du modèle à partir des données $ D_n $ a été écrite très simplement, mais il semble que l'on puisse concevoir diverses mesures telles que des mesures lorsque la probabilité devient très faible et s'accélère.
L'estimation de la distribution de probabilité derrière l'ensemble d'échantillons dispersés en tant que groupe de points semble être quelque peu similaire aux caractéristiques des particules et des ondes de la lumière, et j'estime qu'elles sont liées.
Que se passe-t-il lorsque la distribution de probabilité change avec le temps? Est également intéressant.

Lien de référence

Je me suis référé à la page suivante.

[PYTHON] Essayez d'estimer les paramètres de la distribution gamma tout en implémentant simplement MCMC