Estimation la plus probable de diverses distributions avec Pyro

Estimation la plus probable de divers paramètres de distribution de probabilité continue unidimensionnelle avec Pyro.

Bien que certaines des estimations les plus probables fournissent des solutions analytiques (comme la distribution normale), Certains sont itératifs (comme la distribution gamma). De plus, la méthode diffère selon la distribution et l'implémentation est différente. Faites-le ici en utilisant pyro, Définissez simplement la distribution de probabilité dans le modèle Après cela, l'estimation la plus probable est effectuée par la méthode itérative avec le même code.

Seule la distribution normale sera expliquée en détail. D'autres codes peuvent être trouvés dans l'essentiel ci-dessous. https://gist.github.com/tttamaki/b061f64ad1c0f640acb2bccb88b5087e

Tips

--Utilisez ʻAutoDelta` comme guide (estimation ponctuelle) --Le moniteur pendant l'apprentissage est accroché

Préparation

Préparation. Pour ordinateur portable


import matplotlib.pyplot as plt
%matplotlib inline

from collections import defaultdict
import os

import numpy as np


import torch
import torch.distributions.constraints as constraints

import pyro
from pyro.optim import Adam
import pyro.distributions as dist
from pyro.infer import SVI, Trace_ELBO
from pyro.infer.autoguide import AutoDelta, AutoNormal
from pyro.infer.autoguide.initialization import init_to_feasible

# import pyro.poutine as poutine
from pyro import poutine

pyro.enable_validation(True)
pyro.set_rng_seed(101)


from tqdm.notebook import tqdm
import daft

Distribution normale unidimensionnelle

Vrais paramètres

Créer aléatoirement de vrais m et std


true_m = np.random.rand() * 10
true_std = np.abs(np.random.rand() * 2)
print('m =', true_m, 'std =', true_std)

x_range = np.arange(true_m - 4*true_std, true_m + 4*true_std, 0.01)
x_max = x_range.max()
x_min = x_range.min()
print('x_range: from {0:.3f} to {1:.3f}'.format(x_min, x_max))

python


m = 2.3235366181476067 std = 0.16712286732668735
x_range: from 1.655 to 2.985

Création de données

Échantillonnage de données à partir d'une vraie distribution


# sampling
target_dist = dist.Normal(true_m, true_std)
data = torch.tensor([target_dist() for i in range(100)])
#pdf: fonction de distribution
true_y = [target_dist.log_prob(torch.tensor([x])).exp() for x in x_range]


fig = plt.figure(figsize=(10,5))
plt.plot(x_range, true_y, c='k', label=r'Normal(m,$\sigma$)')
plt.hist(data, range=(x_min, x_max), bins=100, density=True, alpha=0.2, label=r'observed data $\{x_i\}$')

plt.title(r'm={0:.3f} $\sigma$={1:.3f}'.format(true_m, true_std))

plt.ylim(0,)
plt.xlabel('x')
plt.legend()
plt.savefig('obs-data')
plt.show()

normal_obs-data.png

À propos, la solution analytique de l'estimation la plus probable est la moyenne et la variance dans le cas de la distribution normale.

python


ml_m = data.numpy().mean()
ml_std = data.numpy().std()
print('m={0:.3f}, sigma={1:.3f}'.format(ml_m, ml_std))

python


m=2.318, sigma=0.157

Modèle graphique

Le modèle «model» a la structure suivante.

normal-graph.png

Dessinez un modèle graphique avec daft


pgm = daft.PGM()

pgm.add_node("xn", r"$x_n$", 2, 1, observed=True)
pgm.add_node("m", "m", 1.5, 2, fixed=True)
pgm.add_node("std", r"$\sigma$", 2.5, 2, fixed=True)
pgm.add_edge("m", "xn")
pgm.add_edge("std", "xn")
pgm.add_plate([1, 0.5, 2, 1], label=r"$n = 1, \ldots, N$", shift=-0.2)

pgm.show()  # render and show

Modèles et guides

Pour le modèle, définissez m et std comme paramètres dans pyro.param.

modèle


def model(data=None):
        
    #La valeur initiale estimée est 0.Mettre à 0
    m = pyro.param("m", torch.tensor(0.0))

    #La valeur initiale estimée est 1.Mettre à 0
    # std >Puisqu'il vaut 0, ajoutez une contrainte positive
    std = pyro.param("std", torch.tensor(1.0), constraint=constraints.positive)

    #data représente toutes les données (c.-à-d. tableau de liste vectorielle)
    with pyro.plate('observe_data'):  #Étant donné que chaque obs
        pyro.sample('obs', dist.Normal(m, std), obs=data)  # vector-Utiliser la plaque

Utilisez l'estimation des points AutoDelta pour vous guider

python


guide = AutoDelta(
    poutine.block(model, hide=['obs']),  #Cachez les obs du modèle et utilisez les autres (m et std)
    init_loc_fn=init_to_feasible  #La valeur initiale est une valeur raisonnable
)

Préparation estimée

python


adam_params = {"lr": 0.005, "betas": (0.95, 0.999)}
optimizer = Adam(adam_params)  #Adam pour le moment


svi = SVI(model=model,
          guide=guide,
          optim=optimizer,
          loss=Trace_ELBO()
         )

Initialisez les paramètres et définissez les hooks. Pour enregistrer la valeur au milieu du paramètre.

python


pyro.clear_param_store()  # pyro.get_param_store()Sera vide

svi.loss(model, guide, data)  #Je dois faire ça une fois piro.get_param_store()Est une maman vide

trace_dic = defaultdict(list)  #Sauvegardez la valeur (cette liste sera appelée trace ici)
for name, value in pyro.get_param_store().named_parameters():
    print('tracing', name, type(value), value)
    #hook est appelé pour chaque calcul de gradient. Ignorer la valeur du dégradé grad, utiliser le nom pour param(name)Obtenez la valeur avec et ajoutez-la à la liste de dic
    value.register_hook(lambda grad, name=name: trace_dic[name].append(pyro.param(name).item()))

python


tracing m <class 'torch.Tensor'> tensor(0., requires_grad=True)
tracing std <class 'torch.Tensor'> tensor(0., requires_grad=True)

Estimation

Estimation itérative.

python


num_steps = 5000

with tqdm(range(num_steps)) as pbar:
    for i,p in enumerate(pbar):

        loss = svi.step(data)  #Calcul du gradient, mise à jour

        trace_dic['loss'].append(loss)

        if i > 100 and np.isclose(trace_dic['loss'][-100], loss):
            break  #Convergence si la valeur ne change pas d'il y a 100 fois
            
        if i % 10 == 0:  #Si vous l'affichez à chaque fois, il est trop tôt pour le voir, alors affichez-le une fois toutes les 10 fois
            pbar.set_postfix(loss=loss)

État de convergence

python


fig = plt.figure(figsize=(10,5))
plt.plot(trace_dic['loss'])
plt.title("ELBO")
plt.xlabel("step")
plt.ylabel("loss");
savefig("loss")
plt.show()


fig = plt.figure(figsize=(15,5))

for i,param in enumerate(['m', 'std']):
    plt.subplot(1,2,i+1)
    plt.plot(trace_dic[param])
    plt.ylabel(param)

plt.tight_layout()
savefig("params")
plt.show()

État de diminution de la perte normal_loss.png

Convergence des paramètres m et std normal_params.png

résultat

Comparaison de la valeur estimée la plus probable et de la valeur réelle

python


trace = poutine.trace(model).get_trace()

print('m', trace.nodes['m']['value'].item(), 'true m', true_m)
print('std', trace.nodes['std']['value'].item(), 'true std', true_std)

C'est proche de la vraie valeur.

python


m 5.083353042602539 true m 5.163986277024462
std 1.1655045747756958 true std 1.1413351737362796

Comparaison de la distribution des estimations les plus probables et de la distribution réelle

python


fig = plt.figure(figsize=(10,5))


#Obtenir des paramètres avec trace
trace = poutine.trace(model).get_trace()
m = trace.nodes['m']['value'].item()
std = trace.nodes['std']['value'].item()

# generate a pdf
estimated_dist = dist.Normal(m, std)
y = [estimated_dist.log_prob(torch.tensor([x])).exp() for x in x_range]


# plot
plt.plot(x_range, y, c='k', label='ML estimated pdf')
plt.plot(x_range, true_y, c='r', label=r'true pdf')
plt.hist(data, range=(x_min, x_max), bins=100, alpha=0.2, density=True, l![normal_sampled-pdf.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44835/0d36d454-fce6-9b62-a42c-a8aa546733f6.png)
abel=r'observed data $\{x_i\}$')


plt.xlabel('x')
plt.legend()
savefig("sampled-pdf")
plt.show()

normal_sampled-pdf.png

Échantillon de la distribution estimée et comparaison avec les données originales

python


fig = plt.figure(figsize=(10,5))

obs = []
for _ in range(5000):
    trace = poutine.trace(model).get_trace()
    obs.append(trace.nodes['obs']['value'].item())

plt.xlabel('posterior samples')

plt.hist(data, range=(x_min, x_max), bins=100, alpha=0.2, density=True, label=r'observed data $\{x_i\}$')
plt.hist(obs, range=(x_min, x_max), bins=100, alpha=0.2, density=True, color='r', label=r'sampled data$')

plt.legend()
savefig("pdf-obs")
plt.show()

normal_pdf-obs.png

Autres distributions

Ce qui suit ne montre que les données d'origine (bleu), la vraie distribution (rouge) et la distribution estimée avec le paramètre d'estimation le plus probable (noir).

Le code est https://gist.github.com/tttamaki/b061f64ad1c0f640acb2bccb88b5087e C'est dedans.

Distribution bêta

beta_sampled-pdf.png

distribution t

student-t_sampled-pdf.png

Distribution de Laplace

laplace_sampled-pdf.png

Distribution de Cauchy

cauchy_sampled-pdf.png

Distribution exponentielle

Exponential_sampled-pdf.png

Distribution gamma

gamma_sampled-pdf.png

Distribution Gambel

gumbel_sampled-pdf.png

Recommended Posts

Estimation la plus probable de diverses distributions avec Pyro
Estimation la plus probable de la moyenne et de la variance avec TensorFlow
Avantages et inconvénients de la méthode d'estimation la plus probable
Implémentation d'estimation la plus probable du modèle de sujet en python
1. Statistiques apprises avec Python 1-3. Calcul de diverses statistiques (statistiques)
Exemple de code python pour la distribution exponentielle et l'estimation la plus probable (MLE)
Essayons à nouveau Estimation de la plupart des probabilités et ajustement du modèle (distribution de probabilité) ① Distribution de probabilité discrète
Essayons à nouveau La plupart des estimations de probabilité et ajustement du modèle (distribution de probabilité) ② Distribution de probabilité continue
Jusqu'à ce que la méthode d'estimation la plus probable trouve le vrai paramètre