J'ai des données qui sont un mélange 50:50 d'une distribution normale et d'une valeur constante:

numdata = 10000
data = np.random.normal(0.0,1.0,numdata).astype(np.float32)
data[int(numdata/2):] = 0.0
plt.hist(data,30,density=True)

sample data

Ma tâche est d'adapter une densité de mélange à ces données. J'utilise un tfd.Mixture avec tfd.Normal et tfd.Deterministic Le rapport connu (dans le cas des échantillons de données) de Normal à Déterministe est de 0,5 Mon MCMC renvoie à la place un ratio de 0,83 en faveur de Normal.

Existe-t-il une meilleure façon d'ajuster cette distribution avec le bon ratio?

Voici un exemple de code complet:

import os
os.environ['CUDA_VISIBLE_DEVICES'] = '-1'
import tensorflow as tf
import tensorflow_probability as tfp
import matplotlib.pyplot as plt
tfd = tfp.distributions
tfb = tfp.bijectors

import numpy as np
from time import time

numdata = 10000
data = np.random.normal(0.0,1.0,numdata).astype(np.float32)
data[int(numdata/2):] = 0.0
_=plt.hist(data,30,density=True)

root = tfd.JointDistributionCoroutine.Root
def dist_fn(rv_p,rv_mu):
    rv_cat = tfd.Categorical(probs=tf.stack([rv_p, 1.-rv_p],-1))
    rv_norm  = tfd.Normal(rv_mu,1.0)
    rv_zero =  tfd.Deterministic(tf.zeros_like(rv_mu))
    
    rv_mix = tfd.Independent(
                tfd.Mixture(cat=rv_cat,
                            components=[rv_norm,rv_zero]),
                reinterpreted_batch_ndims=1)
    return rv_mix


def model_fn():
    rv_p    = yield root(tfd.Sample(tfd.Uniform(0.0,1.0),1))
    rv_mu   = yield root(tfd.Sample(tfd.Uniform(-1.,1. ),1))
    
    rv_mix  = yield dist_fn(rv_p,rv_mu)
    
jd = tfd.JointDistributionCoroutine(model_fn)
unnormalized_posterior_log_prob = lambda *args: jd.log_prob(args + (data,))

n_chains = 1

p_init = [0.3]
p_init = tf.cast(p_init,dtype=tf.float32)

mu_init = 0.1
mu_init = tf.stack([mu_init]*n_chains,axis=0)

initial_chain_state = [
    p_init,
    mu_init,
]

bijectors = [
    tfb.Sigmoid(),  # p
    tfb.Identity(),  # mu
]

step_size = 0.01

num_results = 50000
num_burnin_steps = 50000


kernel=tfp.mcmc.TransformedTransitionKernel(
    inner_kernel=tfp.mcmc.HamiltonianMonteCarlo(
    target_log_prob_fn=unnormalized_posterior_log_prob,
    num_leapfrog_steps=2,
    step_size=step_size,
    state_gradients_are_stopped=True),
    bijector=bijectors)

kernel = tfp.mcmc.SimpleStepSizeAdaptation(
    inner_kernel=kernel, num_adaptation_steps=int(num_burnin_steps * 0.8))

#XLA optim
@tf.function(autograph=False, experimental_compile=True)
def graph_sample_chain(*args, **kwargs):
  return tfp.mcmc.sample_chain(*args, **kwargs)


st = time()
trace,stats = graph_sample_chain(
      num_results=num_results,
      num_burnin_steps=num_burnin_steps,
      current_state=initial_chain_state,
      kernel=kernel)
et = time()
print(et-st)


ptrace, mutrace = trace
plt.subplot(121)
_=plt.hist(ptrace.numpy(),100,density=True)
plt.subplot(122)
_=plt.hist(mutrace.numpy(),100,density=True)
print(np.mean(ptrace),np.mean(mutrace))

Les distributions résultantes de p et mu sont comme ceci: Histogrammes de p et mu

Évidemment, il devrait avoir une moyenne de p à 0,5. Je soupçonne qu'il y a peut-être quelque chose qui ne va pas avec model_fn (). J'ai essayé d'évaluer log_prob du modèle à différentes valeurs p et en effet le "optimal" est autour de 0,83 Je ne comprends tout simplement pas pourquoi et comment le corriger afin de reconstruire le mélange d'origine.

[EDIT] Un code de démonstration "plus simple" avec pymc3. Toujours le même comportement, le résultat est de 0,83 au lieu de 0,5

import pymc3 as pm
import numpy as np
import arviz as az
import matplotlib.pyplot as plt


numdata = 1000
data1 = np.random.normal(0.0,1.0,numdata).astype(np.float32)
data2 = np.zeros(numdata).astype(np.float32)
data = np.concatenate((data1,data2))


_=plt.hist(data,30,density=True)

with pm.Model() as model:
    norm = pm.Normal.dist(0.0,1.0)
    zero = pm.Constant.dist(0.0)
    
    components = [norm,zero]
    w = pm.Dirichlet('p', a=np.array([1, 1]))  # two mixture component weights.
    like = pm.Mixture('data', w=w, comp_dists=components, observed=data)
    
    posterior = pm.sample()
    
    idata = az.from_pymc3(posterior)
    az.plot_posterior(posterior)
1
rv_normal 28 févr. 2021 à 12:24

1 réponse

Meilleure réponse

Incommensurabilité de la densité et de la masse de probabilité

Le problème ici est que la probabilité de provenir de chaque modèle implique une densité de probabilité pour le gaussien et une masse pour le discret, qui ne sont pas proportionnelles. Plus précisément, le calcul pour comparer l'origine d'une observation nulle impliquera des probabilités

P[x=0|Normal[0,1]] = 1/sqrt(2*pi) = 0.3989422804014327
P[x=0|   Zero    ] = 1

Qui les comparera (pondérés par p) comme s'ils avaient la même unité. Cependant, le premier est une densité et donc infinitésimale par rapport au second. Si l'on ignore cette incommensurabilité, on agit effectivement comme si le gaussien avait 40% de chances de générer des zéros, alors qu'en réalité il presque jamais ne génère exactement un zéro.

Solution de contournement: distribution pseudo-discrète

Nous devons convertir les unités d'une manière ou d'une autre. Une manière simple de faire ceci est d'approximer la distribution discrète avec une distribution continue, de sorte que les probabilités qu'elle génère soient en unités de densité. Par exemple, l'utilisation d'une distribution gaussienne ou de Laplace de haute précision (étroite) centrée sur la valeur discrète donne un postérieur sur p centré autour de 0,5:

with pm.Model() as model:
    norm = pm.Normal.dist(0.0, 1.0)
    pseudo_zero = pm.Laplace.dist(0.0, 1e-16)
    
    components = [norm, pseudo_zero]
    w = pm.Dirichlet('p', a=np.array([1, 1]))  # two mixture component weights.
    like = pm.Mixture('data', w=w, comp_dists=components, observed=data)
    
    posterior = pm.sample()
    
    idata = az.from_pymc3(posterior)
    az.plot_posterior(posterior)

enter image description here


Pourquoi p=0.83?

Le postérieur que nous observons en mélangeant le discret et le continu n'est pas arbitraire. Voici quelques façons de l'obtenir. Pour ce qui suit, nous utiliserons simplement un p pour désigner la probabilité de provenir de la Gaussienne.

Estimation MAP

En ignorant l'incommensurabilité, nous pouvons dériver l'estimation MAP pour p comme suit. Désignons les observations combinées par D = { D_1 | D_2 }, où D_1 est le sous-ensemble de la Gaussienne, etc., et n est le nombre d'observations de chaque sous-ensemble. Ensuite, nous pouvons écrire la probabilité

P[p|D] ~ P[D|p]P[p]

Puisque le Dirichlet est uniforme, nous pouvons ignorer P[p] et étendre nos données

P[D|p] = P[D_1|p]P[D_2|p]
       = (Normal[D_1|0,1]*(p^n))(Normal[0|0,1]*p + 1*(1-p))^n
       = Normal[D_1|0,1]*(p^n)(0.3989*p + 1 - p)^n
       = Normal[D_1|0,1]*(p - 0.6011*(p^2))^n

Prenant le dérivé w.r.t. p et mise à zéro, nous avons

0 = n*(1-1.2021*p)(p-0.6011*p^2)^(n-1)

Qui prend un zéro (non trivial) en p = 1/1.2021 = 0.8318669.

Expérience de la pensée d'échantillonnage

Une autre façon d'aborder cela serait par une expérience d'échantillonnage. Supposons que nous utilisions le schéma suivant pour échantillonner p.

  1. Commencez par un p donné.
  2. Pour chaque observation, dessinez un échantillon de Bernoulli en utilisant la vraisemblance des deux modèles, pondérée par la valeur p précédente.
  3. Calculez un nouveau p comme la moyenne de tous ces tirages de Bernoulli.
  4. Aller à l'étape 1.

Essentiellement, un échantillonneur de Gibbs pour p, mais robuste à des affectations de modèles d'observation impossibles.

Pour la première itération, commençons par p=0.5. Pour toutes les observations véritablement gaussiennes, elles auront une vraisemblance nulle pour le modèle discret, donc, au minimum, la moitié de tous nos tirages de Bernoulli sera 1 (pour le gaussien). Pour toutes les observations issues du modèle discret, il s'agira d'une comparaison de la probabilité d'observer un zéro dans chaque modèle. C'est 1 pour le modèle discret et 0,3989422804014327 pour le gaussien. Normaliser cela signifie que nous avons Bernoulli dessine avec une probabilité de

p*0.3989/((1-p)*1 + p*0.3989)
# 0.2851742248343187

En faveur du gaussien. Maintenant, nous pouvons mettre à jour p, et ici nous allons simplement travailler avec les valeurs attendues, à savoir:

p = 0.5*1 + 0.5*0.2851742248343187
# 0.6425871124171594

Que se passe-t-il si nous commençons à itérer cela?

# likelihood for zero from normal
lnorm = np.exp(pm.Normal.dist(0,1).logp(0).eval())

# history
p_n = np.zeros(101)

# initial value
p_n[0] = 0.5

for i in range(100):
    # update
    p_n[1+i] = 0.5 + 0.5*p_n[i]*lnorm/((1-p_n[i])+p_n[i]*lnorm)

plt.plot(p_n);
p_n[100]
# 0.8318668635076404

enter image description here

Encore une fois, les valeurs attendues convergent vers notre moyenne postérieure de p = 0.83.

Par conséquent, mis à part le fait que les PDF et les PMF ont des unités différentes pour leurs codomaines, il semble que TFP et PyMC3 se comportent correctement.

2
merv 7 mars 2021 à 03:45