이 페이지는 Cloud Translation API를 통해 번역되었습니다.
Switch to English

사전 조건화 된 확률 적 기울기 Langevin Dynamics를 사용하여 Dirichlet 공정 혼합물 모델 피팅

이 노트북에서는 Dirichlet Process Mixture of Gaussian 분포를 피팅하여 많은 수의 샘플을 클러스터링하고 동시에 클러스터 수를 추론하는 방법을 보여줍니다. 추론을 위해 pSGLD (Preconditioned Stochastic Gradient Langevin Dynamics)를 사용합니다.

목차

  1. 샘플

  2. 모델

  3. 최적화

  4. 결과 시각화

    4.1. 클러스터링 된 결과

    4.2. 불확실성 시각화

    4.3. 선택한 혼합물 성분의 평균 및 규모

    4.4. 각 혼합물 성분의 혼합물 중량

    4.5. $ \ alpha $의 수렴

    4.6. 반복에 대한 추론 된 클러스터 수

    4.7. RMSProp을 사용하여 모델 피팅

  5. 결론


1. 샘플

먼저 장난감 데이터 세트를 설정합니다. 세 개의 이변 량 가우스 분포에서 50,000 개의 무작위 샘플을 생성합니다.

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. 모델

여기에서는 대칭 디리클레 사전을 사용하여 가우스 분포의 디리클레 프로세스 혼합물을 정의합니다. 노트북 전체에서 벡터 수량은 굵게 표시됩니다. $ i \ in {1, \ ldots, N} $ 샘플에 대해 $ j \ in {1, \ ldots, K} $ 가우시안 분포가 혼합 된 모델은 다음과 같이 공식화됩니다.

$$\begin{align*} p(\boldsymbol{x}_1,\cdots, \boldsymbol{x}_N) &=\prod_{i=1}^N \text{GMM}(x_i), \\ &\,\quad \text{with}\;\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*}$$

어디:

$$\begin{align*} x_i&\sim \text{Normal}(\text{loc}=\boldsymbol{\mu}_{z_i},\,\text{scale}=\boldsymbol{\sigma}_{z_i}) \\ z_i &= \text{Categorical}(\text{prob}=\boldsymbol{\pi}),\\ &\,\quad \text{with}\;\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*}$$

우리의 목표는 클러스터의 유추 된 인덱스를 나타내는 $ z_i $를 통해 $ j $ 번째 클러스터에 각 $ x_i $를 할당하는 것입니다.

이상적인 Dirichlet 혼합물 모델의 경우 $ K $는 $ \ infty $로 설정됩니다. 그러나 충분히 큰 $ K $로 Dirichlet 혼합물 모델을 근사화 할 수있는 것으로 알려져 있습니다. 초기 값을 $ K $로 임의로 설정했지만, 단순한 가우스 혼합 모델과 달리 최적화를 통해 최적의 클러스터 수를 추론 할 수도 있습니다.

이 노트북에서는 이변 량 가우스 분포를 혼합 성분으로 사용하고 $ 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. 최적화

사전 조건화 된 확률 적 기울기 Langevin Dynamics (pSGLD)를 사용하여 모델을 최적화하여 미니 배치 기울기 하강 방식으로 많은 수의 샘플에 대해 모델을 최적화 할 수 있습니다.

$ t \, $에서 $ \ boldsymbol {\ theta} \ equiv {\ boldsymbol {\ pi}, \, \ alpha, \, \ boldsymbol {\ mu_j}, \, \ boldsymbol {\ sigma_j}} $ 매개 변수를 업데이트하려면 미니 배치 크기 $ M $로 반복하면 업데이트는 다음과 같이 샘플링됩니다.

$$\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*}$$

위의 방정식에서 $ \ epsilon _ {t} $는 $ t \, $ th 반복에서의 학습률이고 $ \ log p (\ theta_t) $는 $ \ theta $의 로그 사전 분포의 합계입니다. $ G (\ boldsymbol {\ theta} _ {t}) $는 각 매개 변수의 그라데이션 배율을 조정하는 전제 조건입니다.

# 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

pSGLD에 대한 손실 함수로 $ \ text {GMM} (x_ {t_k}) $ 및 사전 확률 $ p (\ theta_t) $의 공동 로그 확률을 사용합니다.

pSGLDAPI에 지정된대로 이전 확률의 합을 샘플 크기 $ 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. 결과 시각화

4.1. 클러스터링 된 결과

먼저 클러스터링 결과를 시각화합니다.

각 샘플 $ x_i $를 클러스터 $ j $에 할당하기 위해 $ z_i $의 사후를 다음과 같이 계산합니다.

$$\begin{align*} j = \underset{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

거의 동일한 수의 샘플이 적절한 클러스터에 할당되고 모델이 올바른 클러스터 수를 성공적으로 추론 한 것을 볼 수 있습니다.

4.2. 불확실성 시각화

여기에서는 각 샘플에 대해 시각화하여 클러스터링 결과의 불확실성을 살펴 봅니다.

엔트로피를 사용하여 불확실성을 계산합니다.

$$\begin{align*} \text{Uncertainty}_\text{entropy} = -\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*}$$

pSGLD에서 우리는 각 반복에서 훈련 매개 변수의 값을 사후 분포의 표본으로 취급합니다. 따라서 각 매개 변수에 대해 $ O $ 반복의 값에 대한 엔트로피를 계산합니다. 최종 엔트로피 값은 모든 클러스터 할당의 엔트로피를 평균하여 계산됩니다.

# 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

위 그래프에서 휘도가 낮을수록 불확실성이 커집니다. 클러스터 경계 근처의 샘플이 특히 더 높은 불확실성을 가지고 있음을 알 수 있습니다. 이러한 샘플은 클러스터링하기 어렵 기 때문에 직관적으로 사실입니다.

4.3. 선택한 혼합물 성분의 평균 및 규모

다음으로 선택한 클러스터의 $ \ mu_j $ 및 $ \ 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]


다시 말하지만 $ \ boldsymbol {\ mu_j} $ 및 $ \ boldsymbol {\ sigma_j} $는 실측 값에 가깝습니다.

4.4 각 혼합물 성분의 혼합물 중량

추론 된 혼합물 가중치도 살펴 봅니다.

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

png

우리는 단지 몇 (3) 개의 혼합물 성분 만이 상당한 가중치를 갖고 나머지 가중치는 0에 가까운 값을 가지고 있음을 알 수 있습니다. 이것은 또한 모델이 샘플 분포를 구성하는 정확한 수의 혼합물 성분을 성공적으로 추론했음을 보여줍니다.

4.5. $ \ alpha $의 수렴

Dirichlet 분포의 농도 매개 변수 $ \ alpha $의 수렴을 살펴 봅니다.

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

$ \ alpha $가 작 으면 Dirichlet 혼합 모델에서 예상되는 군집 수가 적다는 사실을 고려할 때이 모델은 반복을 통해 최적의 군집 수를 학습하는 것으로 보입니다.

4.6. 반복에 대한 추론 된 클러스터 수

추론 된 클러스터 수가 반복에 따라 어떻게 변하는 지 시각화합니다.

이를 위해 반복에 대한 클러스터 수를 추론합니다.

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

반복을 통해 클러스터 수가 3 개에 가까워지고 있습니다. 반복을 통해 $ \ alpha $를 더 작은 값으로 수렴 한 결과, 모델이 최적의 클러스터 수를 추론하는 매개 변수를 성공적으로 학습하고 있음을 알 수 있습니다.

흥미롭게도, 추론이 훨씬 이후의 반복에서 수렴 된 $ \ alpha $와 달리 초기 반복에서 올바른 수의 클러스터로 수렴되었음을 알 수 있습니다.

4.7. RMSProp을 사용하여 모델 피팅

이 섹션에서는 pSGLD의 Monte Carlo 샘플링 체계의 효과를 확인하기 위해 RMSProp를 사용하여 모델에 맞 춥니 다. 샘플링 방식이없고 pSGLD가 RMSProp를 기반으로하기 때문에 비교를 위해 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)

pSGLD에 비해 RMSProp에 대한 반복 횟수는 더 길지만 RMSProp에 의한 최적화는 훨씬 빠릅니다.

다음으로 클러스터링 결과를 살펴 봅니다.

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

실험에서 RMSProp 최적화에 의해 클러스터 수가 올바르게 추론되지 않았습니다. 우리는 또한 혼합물 무게를 봅니다.

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

png

잘못된 수의 성분에 상당한 혼합물 중량이 있음을 알 수 있습니다.

최적화에 더 오랜 시간이 걸리지 만 Monte Carlo 샘플링 방식을 사용하는 pSGLD가 실험에서 더 잘 수행되었습니다.

5. 결론

이 노트북에서는 pSGLD를 사용하여 Dirichlet Process Mixture of Gaussian 분포를 피팅하여 많은 수의 샘플을 클러스터링하고 동시에 클러스터 수를 추론하는 방법을 설명했습니다.

이 실험은 모델이 샘플을 성공적으로 클러스터링하고 정확한 클러스터 수를 추론 한 것으로 나타났습니다. 또한 pSGLD의 Monte Carlo 샘플링 방식을 통해 결과의 불확실성을 시각화 할 수 있음을 보여주었습니다. 샘플을 클러스터링 할뿐만 아니라 모델이 혼합물 성분의 올바른 매개 변수를 추론 할 수 있음을 확인했습니다. 매개 변수와 추론 된 클러스터 수 간의 관계에 대해 모델이 𝛼 수렴과 추론 된 클러스터 수 간의 상관 관계를 시각화하여 효과적인 클러스터 수를 제어하는 ​​매개 변수를 학습하는 방법을 조사했습니다. 마지막으로 RMSProp을 사용하여 모델을 피팅 한 결과를 살펴 보았습니다. Monte Carlo 샘플링 체계가없는 최적화 프로그램 인 RMSProp는 pSGLD보다 훨씬 빠르게 작동하지만 클러스터링의 정확도는 떨어집니다.

장난감 데이터 세트에는 2 차원 만있는 샘플이 50,000 개 뿐이지 만 여기에 사용 된 미니 배치 방식 최적화는 훨씬 더 큰 데이터 세트에 맞게 확장 가능합니다.