Cette page a été traduite par l'API Cloud Translation.
Switch to English

Ajustement du modèle de mélange de processus Dirichlet à l'aide de la dynamique de Langevin à gradient stochastique préconditionné

Dans ce cahier, nous montrerons comment regrouper un grand nombre d'échantillons et déduire le nombre de groupes simultanément en ajustant un mélange de processus de Dirichlet de distribution gaussienne. Nous utilisons la dynamique de Langevin à gradient stochastique préconditionné (pSGLD) pour l'inférence.

Table des matières

  1. Échantillons

  2. Modèle

  3. Optimisation

  4. Visualisez le résultat

    4.1. Résultat groupé

    4.2. Visualisez l'incertitude

    4.3. Moyenne et échelle du composant de mélange sélectionné

    4.4. Poids du mélange de chaque composant du mélange

    4.5. Convergence de $ \ alpha $

    4.6. Nombre inféré de clusters au cours des itérations

    4.7. Ajuster le modèle avec RMSProp

  5. Conclusion


1. Échantillons

Tout d'abord, nous avons créé un jeu de données sur les jouets. Nous générons 50 000 échantillons aléatoires à partir de trois distributions gaussiennes bivariées.

 from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import time
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import tensorflow_probability as tfp
 
 plt.style.use('ggplot')
tfd = tfp.distributions
 
 def session_options(enable_gpu_ram_resizing=True):
  """Convenience function which sets common `tf.Session` options."""
  config = tf.ConfigProto()
  config.log_device_placement = True
  if enable_gpu_ram_resizing:
    # `allow_growth=True` makes it possible to connect multiple colabs to your
    # GPU. Otherwise the colab malloc's all GPU ram.
    config.gpu_options.allow_growth = True
  return config

def reset_sess(config=None):
  """Convenience function to create the TF graph and session, or reset them."""
  if config is None:
    config = session_options()
  tf.reset_default_graph()
  global sess
  try:
    sess.close()
  except:
    pass
  sess = tf.InteractiveSession(config=config)
 
 # For reproducibility
rng = np.random.RandomState(seed=45)
tf.set_random_seed(76)

# Precision
dtype = np.float64

# Number of training samples
num_samples = 50000

# Ground truth loc values which we will infer later on. The scale is 1.
true_loc = np.array([[-4, -4],
                     [0, 0],
                     [4, 4]], dtype)

true_components_num, dims = true_loc.shape

# Generate training samples from ground truth loc
true_hidden_component = rng.randint(0, true_components_num, num_samples)
observations = (true_loc[true_hidden_component]

                + rng.randn(num_samples, dims).astype(dtype))
 
 # Visualize samples
plt.scatter(observations[:, 0], observations[:, 1], 1)
plt.axis([-10, 10, -10, 10])
plt.show()
 

png

2. Modèle

Ici, nous définissons un mélange de processus de Dirichlet de distribution gaussienne avec Symmetric Dirichlet Prior. Tout au long du cahier, les quantités vectorielles sont écrites en gras. Sur $ i \ in {1, \ ldots, N} $ samples, le modèle avec un mélange de $ j \ in {1, \ ldots, K} $ distributions gaussiennes est formulé comme suit:

$$ \ begin {align *} p (\ boldsymbol {x} _1, \ cdots, \ boldsymbol {x} _N) & = \ prod_ {i = 1} ^ N \ text {GMM} (x_i), \\ & \, \ quad \ text {avec} \; \ text {GMM} (x_i) = \ sum_ {j = 1} ^ K \ pi_j \ text {Normal} (x_i \, | \, \ text {loc} = \ boldsymbol {\ mu_ {j}}, \, \ text {scale} = \ boldsymbol {\ sigma_ {j}}) \\ \ end {align *} $$

où:

$$ \ begin {align *} x_i & \ sim \ text {Normal} (\ text {loc} = \ boldsymbol {\ mu} _ {z_i}, \, \ text {scale} = \ boldsymbol {\ sigma} _ { z_i}) \\ z_i & = \ text {Catégorique} (\ text {prob} = \ boldsymbol {\ pi}), \\ & \, \ quad \ text {avec} \; \ boldsymbol {\ pi} = \ {\ pi_1, \ cdots, \ pi_K \} \\ \ boldsymbol {\ pi} & \ sim \ text {Dirichlet} (\ text {concentration} = \ {\ frac {\ alpha} {K}, \ cdots, \ frac {\ alpha} {K} \}) \\ \ alpha & \ sim \ text {InverseGamma} (\ text {concentration} = 1, \, \ text {rate} = 1) \\ \ boldsymbol {\ mu_j} & \ sim \ text {Normal} (\ text {loc} = \ boldsymbol {0}, \, \ text {scale} = \ boldsymbol {1}) \\ \ boldsymbol {\ sigma_j} & \ sim \ text {InverseGamma} (\ text {concentration} = \ boldsymbol {1}, \, \ text {rate} = \ boldsymbol {1}) \\ \ end {align *} $$

Notre objectif est d'assigner chaque $ x_i $ au $ j $ ème cluster via $ z_i $ qui représente l'index inféré d'un cluster.

Pour un modèle de mélange Dirichlet idéal, $ K $ est défini sur $ \ infty $. Cependant, on sait que l'on peut approcher un modèle de mélange de Dirichlet avec un $ K $ suffisamment grand. Notez que bien que nous ayons fixé arbitrairement une valeur initiale de $ K $, un nombre optimal de clusters est également déduit par optimisation, contrairement à un simple modèle de mélange gaussien.

Dans ce cahier, nous utilisons une distribution gaussienne bivariée comme composant de mélange et définissons $ K $ à 30.

 reset_sess()

# Upperbound on K
max_cluster_num = 30

# Define trainable variables.
mix_probs = tf.nn.softmax(
    tf.Variable(
        name='mix_probs',
        initial_value=np.ones([max_cluster_num], dtype) / max_cluster_num))

loc = tf.Variable(
    name='loc',
    initial_value=np.random.uniform(
        low=-9, #set around minimum value of sample value
        high=9, #set around maximum value of sample value
        size=[max_cluster_num, dims]))

precision = tf.nn.softplus(tf.Variable(
    name='precision',
    initial_value=
    np.ones([max_cluster_num, dims], dtype=dtype)))

alpha = tf.nn.softplus(tf.Variable(
    name='alpha',
    initial_value=
    np.ones([1], dtype=dtype)))

training_vals = [mix_probs, alpha, loc, precision]


# Prior distributions of the training variables

#Use symmetric Dirichlet prior as finite approximation of Dirichlet process.
rv_symmetric_dirichlet_process = tfd.Dirichlet(
    concentration=np.ones(max_cluster_num, dtype) * alpha / max_cluster_num,
    name='rv_sdp')

rv_loc = tfd.Independent(
    tfd.Normal(
        loc=tf.zeros([max_cluster_num, dims], dtype=dtype),
        scale=tf.ones([max_cluster_num, dims], dtype=dtype)),
    reinterpreted_batch_ndims=1,
    name='rv_loc')


rv_precision = tfd.Independent(
    tfd.InverseGamma(
        concentration=np.ones([max_cluster_num, dims], dtype),
        rate=np.ones([max_cluster_num, dims], dtype)),
    reinterpreted_batch_ndims=1,
    name='rv_precision')

rv_alpha = tfd.InverseGamma(
    concentration=np.ones([1], dtype=dtype),
    rate=np.ones([1]),
    name='rv_alpha')

# Define mixture model
rv_observations = tfd.MixtureSameFamily(
    mixture_distribution=tfd.Categorical(probs=mix_probs),
    components_distribution=tfd.MultivariateNormalDiag(
        loc=loc,
        scale_diag=precision))
 

3. Optimisation

Nous optimisons le modèle avec Preconditioned Stochastic Gradient Langevin Dynamics (pSGLD), ce qui nous permet d'optimiser un modèle sur un grand nombre d'échantillons par descente de gradient en mini-lots.

Pour mettre à jour les paramètres $ \ boldsymbol {\ theta} \ equiv {\ boldsymbol {\ pi}, \, \ alpha, \, \ boldsymbol {\ mu_j}, \, \ boldsymbol {\ sigma_j}} $ in $ t \, $ ème itération avec une taille de mini-lot $ M $, la mise à jour est échantillonnée comme:

$$ \ begin {align *} \ Delta \ boldsymbol {\ theta} _ {t} & \ sim \ frac {\ epsilon _ {t}} {2} \ bigl [G \ left (\ boldsymbol {\ theta} _ {t} \ right) \ bigl (\ nabla _ {\ boldsymbol {\ theta}} \ log p \ left (\ boldsymbol {\ theta} _ {t} \ right) + \ frac {N} {M} \ sum _ {k = 1} ^ {M} \ nabla _ \ boldsymbol {\ theta} \ log \ text {GMM} (x_ {t_k}) \ bigr) + \ sum_ \ boldsymbol {\ theta} \ nabla_ \ theta G \ left (\ boldsymbol {\ theta} _ {t} \ right) \ bigr] \\ & + G ^ {\ frac {1} {2}} \ left (\ boldsymbol {\ theta} _ {t} \ right) \ text {Normal} \ left (\ text {loc} = \ boldsymbol {0}, \, \ text {scale} = \ epsilon _ {t} \ boldsymbol {1} ​​\ right) \\ \ end {align *} $$

Dans l'équation ci-dessus, $ \ epsilon _ {t} $ est le taux d'apprentissage à $ t \, $ ème itération et $ \ log p (\ theta_t) $ est une somme des distributions antérieures log de $ \ theta $. $ G (\ boldsymbol {\ theta} _ {t}) $ est un préconditionneur qui ajuste l'échelle du gradient de chaque paramètre.

 # Learning rates and decay
starter_learning_rate = 1e-6
end_learning_rate = 1e-10
decay_steps = 1e4

# Number of training steps
training_steps = 10000

# Mini-batch size
batch_size = 20

# Sample size for parameter posteriors
sample_size = 100
 

Nous utiliserons la probabilité logarithmique conjointe de la vraisemblance $ \ text {GMM} (x_ {t_k}) $ et les probabilités antérieures $ p (\ theta_t) $ comme fonction de perte pour pSGLD.

Notez que comme spécifié dans l' API de pSGLD , nous devons diviser la somme des probabilités antérieures par la taille de l'échantillon $ N $.

 # Placeholder for mini-batch
observations_tensor = tf.compat.v1.placeholder(dtype, shape=[batch_size, dims])

# Define joint log probabilities
# Notice that each prior probability should be divided by num_samples and
# likelihood is divided by batch_size for pSGLD optimization.
log_prob_parts = [
    rv_loc.log_prob(loc) / num_samples,
    rv_precision.log_prob(precision) / num_samples,
    rv_alpha.log_prob(alpha) / num_samples,
    rv_symmetric_dirichlet_process.log_prob(mix_probs)[..., tf.newaxis]
    / num_samples,
    rv_observations.log_prob(observations_tensor) / batch_size
]
joint_log_prob = tf.reduce_sum(tf.concat(log_prob_parts, axis=-1), axis=-1)
 
 # Make mini-batch generator
dx = tf.compat.v1.data.Dataset.from_tensor_slices(observations)\
  .shuffle(500).repeat().batch(batch_size)
iterator = tf.compat.v1.data.make_one_shot_iterator(dx)
next_batch = iterator.get_next()

# Define learning rate scheduling
global_step = tf.Variable(0, trainable=False)
learning_rate = tf.train.polynomial_decay(
    starter_learning_rate,
    global_step, decay_steps,
    end_learning_rate, power=1.)

# Set up the optimizer. Don't forget to set data_size=num_samples.
optimizer_kernel = tfp.optimizer.StochasticGradientLangevinDynamics(
    learning_rate=learning_rate,
    preconditioner_decay_rate=0.99,
    burnin=1500,
    data_size=num_samples)

train_op = optimizer_kernel.minimize(-joint_log_prob)

# Arrays to store samples
mean_mix_probs_mtx = np.zeros([training_steps, max_cluster_num])
mean_alpha_mtx = np.zeros([training_steps, 1])
mean_loc_mtx = np.zeros([training_steps, max_cluster_num, dims])
mean_precision_mtx = np.zeros([training_steps, max_cluster_num, dims])

init = tf.global_variables_initializer()
sess.run(init)

start = time.time()
for it in range(training_steps):
  [
      mean_mix_probs_mtx[it, :],
      mean_alpha_mtx[it, 0],
      mean_loc_mtx[it, :, :],
      mean_precision_mtx[it, :, :],
      _
  ] = sess.run([
      *training_vals,
      train_op
  ], feed_dict={
      observations_tensor: sess.run(next_batch)})

elapsed_time_psgld = time.time() - start
print("Elapsed time: {} seconds".format(elapsed_time_psgld))

# Take mean over the last sample_size iterations
mean_mix_probs_ = mean_mix_probs_mtx[-sample_size:, :].mean(axis=0)
mean_alpha_ = mean_alpha_mtx[-sample_size:, :].mean(axis=0)
mean_loc_ = mean_loc_mtx[-sample_size:, :].mean(axis=0)
mean_precision_ = mean_precision_mtx[-sample_size:, :].mean(axis=0)
 
Elapsed time: 309.8013095855713 seconds

4. Visualisez le résultat

4.1. Résultat groupé

Tout d'abord, nous visualisons le résultat du clustering.

Pour affecter chaque échantillon $ x_i $ à un cluster $ j $, nous calculons le postérieur de $ z_i $ comme:

$$ \ begin {align *} j = \ underet {z_i} {\ arg \ max} \, p (z_i \, | \, x_i, \, \ boldsymbol {\ theta}) \ end {align *} $$
 loc_for_posterior = tf.compat.v1.placeholder(
    dtype, [None, max_cluster_num, dims], name='loc_for_posterior')
precision_for_posterior = tf.compat.v1.placeholder(
    dtype, [None, max_cluster_num, dims], name='precision_for_posterior')
mix_probs_for_posterior = tf.compat.v1.placeholder(
    dtype, [None, max_cluster_num], name='mix_probs_for_posterior')

# Posterior of z (unnormalized)
unnomarlized_posterior = tfd.MultivariateNormalDiag(
    loc=loc_for_posterior, scale_diag=precision_for_posterior)\
   .log_prob(tf.expand_dims(tf.expand_dims(observations, axis=1), axis=1))\

   + tf.log(mix_probs_for_posterior[tf.newaxis, ...])

# Posterior of z (normarizad over latent states)
posterior = unnomarlized_posterior\

  - tf.reduce_logsumexp(unnomarlized_posterior, axis=-1)[..., tf.newaxis]

cluster_asgmt = sess.run(tf.argmax(
    tf.reduce_mean(posterior, axis=1), axis=1), feed_dict={
        loc_for_posterior: mean_loc_mtx[-sample_size:, :],
        precision_for_posterior: mean_precision_mtx[-sample_size:, :],
        mix_probs_for_posterior: mean_mix_probs_mtx[-sample_size:, :]})

idxs, count = np.unique(cluster_asgmt, return_counts=True)

print('Number of inferred clusters = {}\n'.format(len(count)))
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})

print('Number of elements in each cluster = {}\n'.format(count))

def convert_int_elements_to_consecutive_numbers_in(array):
  unique_int_elements = np.unique(array)
  for consecutive_number, unique_int_element in enumerate(unique_int_elements):
    array[array == unique_int_element] = consecutive_number
  return array

cmap = plt.get_cmap('tab10')
plt.scatter(
    observations[:, 0], observations[:, 1],
    1,
    c=cmap(convert_int_elements_to_consecutive_numbers_in(cluster_asgmt)))
plt.axis([-10, 10, -10, 10])
plt.show()
 
Number of inferred clusters = 3

Number of elements in each cluster = [16911 16645 16444]


png

Nous pouvons voir qu'un nombre presque égal d'échantillons sont assignés aux grappes appropriées et le modèle a également réussi à déduire le nombre correct de grappes.

4.2. Visualisez l'incertitude

Ici, nous examinons l'incertitude du résultat du regroupement en le visualisant pour chaque échantillon.

Nous calculons l'incertitude en utilisant l'entropie:

$$ \ begin {align *} \ text {Incertitude} _ \ text {entropie} = - \ frac {1} {K} \ sum ^ {K} _ {z_i = 1} \ sum ^ {O} _ {l = 1} p (z_i \, | \, x_i, \, \ boldsymbol {\ theta} _l) \ log p (z_i \, | \, x_i, \, \ boldsymbol {\ theta} _l) \ end {align * } $$

Dans pSGLD, nous traitons la valeur d'un paramètre d'apprentissage à chaque itération comme un échantillon de sa distribution postérieure. Ainsi, nous calculons l'entropie sur des valeurs à partir d'itérations $ O $ pour chaque paramètre. La valeur d'entropie finale est calculée en faisant la moyenne des entropies de toutes les affectations de cluster.

 # Calculate entropy
posterior_in_exponential = tf.exp(posterior)
uncertainty_in_entropy = tf.reduce_mean(-tf.reduce_sum(
    posterior_in_exponential

    * posterior,
    axis=1), axis=1)

uncertainty_in_entropy_ = sess.run(uncertainty_in_entropy, feed_dict={
    loc_for_posterior: mean_loc_mtx[-sample_size:, :],
    precision_for_posterior: mean_precision_mtx[-sample_size:, :],
    mix_probs_for_posterior: mean_mix_probs_mtx[-sample_size:, :]
})
 
 plt.title('Entropy')
sc = plt.scatter(observations[:, 0],
                 observations[:, 1],
                 1,
                 c=uncertainty_in_entropy_,
                 cmap=plt.cm.viridis_r)
cbar = plt.colorbar(sc,
                    fraction=0.046,
                    pad=0.04,
                    ticks=[uncertainty_in_entropy_.min(),
                           uncertainty_in_entropy_.max()])
cbar.ax.set_yticklabels(['low', 'high'])
cbar.set_label('Uncertainty', rotation=270)
plt.show()
 

png

Dans le graphique ci-dessus, moins de luminance représente plus d'incertitude. Nous pouvons voir que les échantillons près des limites des grappes ont une incertitude particulièrement élevée. C'est intuitivement vrai, que ces échantillons sont difficiles à regrouper.

4.3. Moyenne et échelle du composant de mélange sélectionné

Ensuite, nous examinons les clusters sélectionnés '$ \ mu_j $ et $ \ sigma_j $.

 for idx, numbe_of_samples in zip(idxs, count):
  print(
      'Component id = {}, Number of elements = {}'
      .format(idx, numbe_of_samples))
  print(
      'Mean loc = {}, Mean scale = {}\n'
      .format(mean_loc_[idx, :], mean_precision_[idx, :]))
 
Component id = 0, Number of elements = 16911
Mean loc = [-4.030 -4.113], Mean scale = [ 0.994  0.972]

Component id = 4, Number of elements = 16645
Mean loc = [ 3.999  4.069], Mean scale = [ 1.038  1.046]

Component id = 5, Number of elements = 16444
Mean loc = [-0.005 -0.023], Mean scale = [ 0.967  1.025]


Encore une fois, les $ \ boldsymbol {\ mu_j} $ et $ \ boldsymbol {\ sigma_j} $ se rapprochent de la vérité fondamentale.

4.4 Poids du mélange de chaque composant du mélange

Nous examinons également les poids de mélange inférés.

 plt.ylabel('Mean posterior of mixture weight')
plt.xlabel('Component')
plt.bar(range(0, max_cluster_num), mean_mix_probs_)
plt.show()
 

png

Nous voyons que quelques (trois) composants de mélange ont des poids significatifs et le reste des poids ont des valeurs proches de zéro. Cela montre également que le modèle a réussi à déduire le nombre correct de composants du mélange qui constitue la distribution des échantillons.

4.5. Convergence de $ \ alpha $

Nous regardons la convergence du paramètre de concentration $ \ alpha $ de la distribution de Dirichlet.

 print('Value of inferred alpha = {0:.3f}\n'.format(mean_alpha_[0]))
plt.ylabel('Sample value of alpha')
plt.xlabel('Iteration')
plt.plot(mean_alpha_mtx)
plt.show()
 
Value of inferred alpha = 0.679


png

Compte tenu du fait que $ \ alpha $ plus petit entraîne un nombre moins attendu de clusters dans un modèle de mélange de Dirichlet, le modèle semble apprendre le nombre optimal de clusters sur les itérations.

4.6. Nombre inféré de clusters au cours des itérations

Nous visualisons comment le nombre inféré de clusters change au fil des itérations.

Pour ce faire, nous déduisons le nombre de clusters sur les itérations.

 step = sample_size
num_of_iterations = 50
estimated_num_of_clusters = []
interval = (training_steps - step) // (num_of_iterations - 1)
iterations = np.asarray(range(step, training_steps+1, interval))
for iteration in iterations:
  start_position = iteration-step
  end_position = iteration

  result = sess.run(tf.argmax(
      tf.reduce_mean(posterior, axis=1), axis=1), feed_dict={
          loc_for_posterior:
              mean_loc_mtx[start_position:end_position, :],
          precision_for_posterior:
              mean_precision_mtx[start_position:end_position, :],
          mix_probs_for_posterior:
              mean_mix_probs_mtx[start_position:end_position, :]})

  idxs, count = np.unique(result, return_counts=True)
  estimated_num_of_clusters.append(len(count))
 
 plt.ylabel('Number of inferred clusters')
plt.xlabel('Iteration')
plt.yticks(np.arange(1, max(estimated_num_of_clusters) + 1, 1))
plt.plot(iterations - 1, estimated_num_of_clusters)
plt.show()
 

png

Au fil des itérations, le nombre de grappes se rapproche de trois. Avec le résultat de la convergence de $ \ alpha $ vers une valeur plus petite sur les itérations, nous pouvons voir que le modèle apprend avec succès les paramètres pour déduire un nombre optimal de clusters.

Fait intéressant, nous pouvons voir que l'inférence a déjà convergé vers le nombre correct de clusters dans les premières itérations, contrairement à $ \ alpha $ qui a convergé dans des itérations beaucoup plus tardives.

4.7. Ajuster le modèle avec RMSProp

Dans cette section, pour voir l'efficacité du schéma d'échantillonnage Monte Carlo de pSGLD, nous utilisons RMSProp pour ajuster le modèle. Nous choisissons RMSProp pour la comparaison car il est livré sans le schéma d'échantillonnage et pSGLD est basé sur RMSProp.

 # Learning rates and decay
starter_learning_rate_rmsprop = 1e-2
end_learning_rate_rmsprop = 1e-4
decay_steps_rmsprop = 1e4

# Number of training steps
training_steps_rmsprop = 50000

# Mini-batch size
batch_size_rmsprop = 20
 
 # Define trainable variables.
mix_probs_rmsprop = tf.nn.softmax(
    tf.Variable(
        name='mix_probs_rmsprop',
        initial_value=np.ones([max_cluster_num], dtype) / max_cluster_num))

loc_rmsprop = tf.Variable(
    name='loc_rmsprop',
    initial_value=np.zeros([max_cluster_num, dims], dtype)

    + np.random.uniform(
        low=-9, #set around minimum value of sample value
        high=9, #set around maximum value of sample value
        size=[max_cluster_num, dims]))

precision_rmsprop = tf.nn.softplus(tf.Variable(
    name='precision_rmsprop',
    initial_value=
    np.ones([max_cluster_num, dims], dtype=dtype)))

alpha_rmsprop = tf.nn.softplus(tf.Variable(
    name='alpha_rmsprop',
    initial_value=
    np.ones([1], dtype=dtype)))

training_vals_rmsprop =\
    [mix_probs_rmsprop, alpha_rmsprop, loc_rmsprop, precision_rmsprop]

# Prior distributions of the training variables

#Use symmetric Dirichlet prior as finite approximation of Dirichlet process.
rv_symmetric_dirichlet_process_rmsprop = tfd.Dirichlet(
    concentration=np.ones(max_cluster_num, dtype)

    * alpha_rmsprop / max_cluster_num,
    name='rv_sdp_rmsprop')

rv_loc_rmsprop = tfd.Independent(
    tfd.Normal(
        loc=tf.zeros([max_cluster_num, dims], dtype=dtype),
        scale=tf.ones([max_cluster_num, dims], dtype=dtype)),
    reinterpreted_batch_ndims=1,
    name='rv_loc_rmsprop')


rv_precision_rmsprop = tfd.Independent(
    tfd.InverseGamma(
        concentration=np.ones([max_cluster_num, dims], dtype),
        rate=np.ones([max_cluster_num, dims], dtype)),
    reinterpreted_batch_ndims=1,
    name='rv_precision_rmsprop')

rv_alpha_rmsprop = tfd.InverseGamma(
    concentration=np.ones([1], dtype=dtype),
    rate=np.ones([1]),
    name='rv_alpha_rmsprop')

# Define mixture model
rv_observations_rmsprop = tfd.MixtureSameFamily(
    mixture_distribution=tfd.Categorical(probs=mix_probs_rmsprop),
    components_distribution=tfd.MultivariateNormalDiag(
        loc=loc_rmsprop,
        scale_diag=precision_rmsprop))
 
 og_prob_parts_rmsprop = [
    rv_loc_rmsprop.log_prob(loc_rmsprop),
    rv_precision_rmsprop.log_prob(precision_rmsprop),
    rv_alpha_rmsprop.log_prob(alpha_rmsprop),
    rv_symmetric_dirichlet_process_rmsprop
        .log_prob(mix_probs_rmsprop)[..., tf.newaxis],
    rv_observations_rmsprop.log_prob(observations_tensor)

    * num_samples / batch_size
]
joint_log_prob_rmsprop = tf.reduce_sum(
    tf.concat(log_prob_parts_rmsprop, axis=-1), axis=-1)
 
 # Define learning rate scheduling
global_step_rmsprop = tf.Variable(0, trainable=False)
learning_rate = tf.train.polynomial_decay(
    starter_learning_rate_rmsprop,
    global_step_rmsprop, decay_steps_rmsprop,
    end_learning_rate_rmsprop, power=1.)

# Set up the optimizer. Don't forget to set data_size=num_samples.
optimizer_kernel_rmsprop = tf.train.RMSPropOptimizer(
    learning_rate=learning_rate,
    decay=0.99)

train_op_rmsprop = optimizer_kernel_rmsprop.minimize(-joint_log_prob_rmsprop)

init_rmsprop = tf.global_variables_initializer()
sess.run(init_rmsprop)

start = time.time()
for it in range(training_steps_rmsprop):
  [
      _
  ] = sess.run([
      train_op_rmsprop
  ], feed_dict={
      observations_tensor: sess.run(next_batch)})

elapsed_time_rmsprop = time.time() - start
print("RMSProp elapsed_time: {} seconds ({} iterations)"
      .format(elapsed_time_rmsprop, training_steps_rmsprop))
print("pSGLD elapsed_time: {} seconds ({} iterations)"
      .format(elapsed_time_psgld, training_steps))

mix_probs_rmsprop_, alpha_rmsprop_, loc_rmsprop_, precision_rmsprop_ =\
  sess.run(training_vals_rmsprop)
 
RMSProp elapsed_time: 53.7574200630188 seconds (50000 iterations)
pSGLD elapsed_time: 309.8013095855713 seconds (10000 iterations)

Comparé à pSGLD, bien que le nombre d'itérations pour RMSProp soit plus long, l'optimisation par RMSProp est beaucoup plus rapide.

Ensuite, nous examinons le résultat du clustering.

 cluster_asgmt_rmsprop = sess.run(tf.argmax(
    tf.reduce_mean(posterior, axis=1), axis=1), feed_dict={
        loc_for_posterior: loc_rmsprop_[tf.newaxis, :],
        precision_for_posterior: precision_rmsprop_[tf.newaxis, :],
        mix_probs_for_posterior: mix_probs_rmsprop_[tf.newaxis, :]})

idxs, count = np.unique(cluster_asgmt_rmsprop, return_counts=True)

print('Number of inferred clusters = {}\n'.format(len(count)))
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})

print('Number of elements in each cluster = {}\n'.format(count))

cmap = plt.get_cmap('tab10')
plt.scatter(
    observations[:, 0], observations[:, 1],
    1,
    c=cmap(convert_int_elements_to_consecutive_numbers_in(
        cluster_asgmt_rmsprop)))
plt.axis([-10, 10, -10, 10])
plt.show()
 
Number of inferred clusters = 4

Number of elements in each cluster = [ 1644 15267 16647 16442]


png

Le nombre de clusters n'a pas été correctement déduit par l'optimisation RMSProp dans notre expérience. Nous regardons également le poids du mélange.

 plt.ylabel('MAP inferece of mixture weight')
plt.xlabel('Component')
plt.bar(range(0, max_cluster_num), mix_probs_rmsprop_)
plt.show()
 

png

Nous pouvons voir que le nombre incorrect de composants a des poids de mélange significatifs.

Bien que l'optimisation prenne plus de temps, pSGLD, qui a un schéma d'échantillonnage de Monte Carlo, a mieux performé dans notre expérience.

5. Conclusion

Dans ce cahier, nous avons décrit comment regrouper un grand nombre d'échantillons ainsi que pour déduire le nombre de clusters simultanément en ajustant un mélange de processus Dirichlet de distribution gaussienne à l'aide de pSGLD.

L'expérience a montré que le modèle a réussi à regrouper les échantillons et à déduire le nombre correct de groupes. De plus, nous avons montré que le schéma d'échantillonnage de Monte Carlo de pSGLD nous permet de visualiser l'incertitude du résultat. Non seulement en regroupant les échantillons, nous avons également vu que le modèle pouvait déduire les paramètres corrects des composants du mélange. Sur la relation entre les paramètres et le nombre de clusters inférés, nous avons étudié comment le modèle apprend le paramètre pour contrôler le nombre de clusters effectifs en visualisant la corrélation entre la convergence de 𝛼 et le nombre de clusters inférés. Enfin, nous avons examiné les résultats de l'ajustement du modèle à l'aide de RMSProp. Nous avons vu RMSProp, qui est l'optimiseur sans schéma d'échantillonnage Monte Carlo, fonctionne beaucoup plus rapidement que pSGLD mais a produit moins de précision dans le clustering.

Bien que l'ensemble de données de jouets ne contienne que 50 000 échantillons avec seulement deux dimensions, l'optimisation de la manière par mini-lots utilisée ici est évolutive pour des ensembles de données beaucoup plus volumineux.