Bayesian গাউসিয়ান মিশ্রণ মডেল এবং হ্যামিল্টোনিয়ান MCMC

TensorFlow.org এ দেখুন Google Colab-এ চালান GitHub-এ উৎস দেখুন নোটবুক ডাউনলোড করুন

এই কোল্যাবে আমরা শুধুমাত্র TensorFlow সম্ভাব্যতা আদিম ব্যবহার করে একটি Bayesian Gaussian Mixture Model (BGMM) এর পশ্চাৎভাগ থেকে নমুনা অনুসন্ধান করব।

মডেল

জন্য \(k\in\{1,\ldots, K\}\) মিশ্রণ উপাদান মাত্রা প্রতিটি \(D\), আমরা মডেল চাই \(i\in\{1,\ldots,N\}\) IID নিম্নলিখিত Bayesian গসিয়ান মিশ্রণ মডেল ব্যবহার করে নমুনা:

\[\begin{align*} \theta &\sim \text{Dirichlet}(\text{concentration}=\alpha_0)\\ \mu_k &\sim \text{Normal}(\text{loc}=\mu_{0k}, \text{scale}=I_D)\\ T_k &\sim \text{Wishart}(\text{df}=5, \text{scale}=I_D)\\ Z_i &\sim \text{Categorical}(\text{probs}=\theta)\\ Y_i &\sim \text{Normal}(\text{loc}=\mu_{z_i}, \text{scale}=T_{z_i}^{-1/2})\\ \end{align*}\]

মনে রাখবেন, scale আর্গুমেন্ট সব আছে cholesky শব্দার্থবিদ্যা। আমরা এই কনভেনশনটি ব্যবহার করি কারণ এটি টিএফ ডিস্ট্রিবিউশনের (যা নিজেই এই কনভেনশনটি আংশিকভাবে ব্যবহার করে কারণ এটি গণনাগতভাবে সুবিধাজনক)।

আমাদের লক্ষ্য হল পশ্চাৎভাগ থেকে নমুনা তৈরি করা:

\[p\left(\theta, \{\mu_k, T_k\}_{k=1}^K \Big| \{y_i\}_{i=1}^N, \alpha_0, \{\mu_{ok}\}_{k=1}^K\right)\]

লক্ষ করুন যে, \(\{Z_i\}_{i=1}^N\) উপস্থিত নেই - শুধুমাত্র ঐ র্যান্ডম ভেরিয়েবল যা দিয়ে স্কেল না আগ্রহী ভাবেই আমরা সবাই \(N\)। (এবং সৌভাগ্য যে সেখানে একটি TF বন্টন যা আউট খর্ব করা পরিচালনা করেন \(Z_i\)।)

গণনাগতভাবে জটিল স্বাভাবিককরণের মেয়াদের কারণে এই বিতরণ থেকে সরাসরি নমুনা করা সম্ভব নয়।

মেট্রোপলিস-হেস্টিংস আলগোরিদিম অবাধ্য-টু-স্বাভাবিক ডিস্ট্রিবিউশন থেকে স্যাম্পলিং জন্য কৌশল আছে।

TensorFlow সম্ভাব্যতা মেট্রোপলিস-হেস্টিংসের উপর ভিত্তি করে বেশ কয়েকটি সহ MCMC বিকল্পের একটি সংখ্যা অফার করে। এই নোটবুক, আমরা ব্যবহার করব হ্যামিল্টনিয়ান মন্টে কার্লো ( tfp.mcmc.HamiltonianMonteCarlo )। এইচএমসি প্রায়শই একটি ভাল পছন্দ কারণ এটি দ্রুত একত্রিত হতে পারে, রাষ্ট্রীয় স্থানকে যৌথভাবে নমুনা করে (সমন্বয় অনুসারে বিরোধিতা করে), এবং টিএফ-এর অন্যতম গুণাবলী: স্বয়ংক্রিয় পার্থক্য। তাই বলা হয়, একটি BGMM অবর থেকে স্যাম্পলিং আসলে ভাল অন্যান্য পন্থা, যেমন, দ্বারা সম্পন্ন করা যেতে পারে গিব এর স্যাম্পলিং ।

%matplotlib inline


import functools

import matplotlib.pyplot as plt; plt.style.use('ggplot')
import numpy as np
import seaborn as sns; sns.set_context('notebook')

import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
import tensorflow_probability as tfp

tfd = tfp.distributions
tfb = tfp.bijectors

physical_devices = tf.config.experimental.list_physical_devices('GPU')
if len(physical_devices) > 0:
  tf.config.experimental.set_memory_growth(physical_devices[0], True)

আসলে মডেল তৈরি করার আগে, আমাদের একটি নতুন ধরনের ডিস্ট্রিবিউশন সংজ্ঞায়িত করতে হবে। উপরের মডেল স্পেসিফিকেশন থেকে, এটা স্পষ্ট যে আমরা একটি বিপরীত কোভেরিয়েন্স ম্যাট্রিক্স, অর্থাৎ [প্রিসিশন ম্যাট্রিক্স](https://en.wikipedia.org/wiki/Precision_(statistics%29) দিয়ে MVN-কে প্যারামিটারাইজ করছি। মেমরি, আমরা আমাদের আনছি করতে হবে Bijector করেছেন। Bijector এগিয়ে রূপান্তর ব্যবহার করব:

  • Y = tf.linalg.triangular_solve((tf.linalg.matrix_transpose(chol_precision_tril), X, adjoint=True) + loc ।

আর log_prob হিসাব শুধু বিপরীত, অর্থাত্:

  • X = tf.linalg.matmul(chol_precision_tril, X - loc, adjoint_a=True) ।

যেহেতু ক্ষেত্রে HMC জন্য আমাদের সকলের দরকার log_prob , এর মানে হল আমরা কখনও কলিং এড়ানোর tf.linalg.triangular_solve (যেমন কেনার ক্ষেত্রে দেখা হবে tfd.MultivariateNormalTriL )। এই সুবিধাজনক যেহেতু হয় tf.linalg.matmul সাধারণত দ্রুত ভাল ক্যাশে এলাকায় বকেয়া করা হয়।

class MVNCholPrecisionTriL(tfd.TransformedDistribution):
  """MVN from loc and (Cholesky) precision matrix."""

  def __init__(self, loc, chol_precision_tril, name=None):
    super(MVNCholPrecisionTriL, self).__init__(
        distribution=tfd.Independent(tfd.Normal(tf.zeros_like(loc),
                                                scale=tf.ones_like(loc)),
                                     reinterpreted_batch_ndims=1),
        bijector=tfb.Chain([
            tfb.Shift(shift=loc),
            tfb.Invert(tfb.ScaleMatvecTriL(scale_tril=chol_precision_tril,
                                           adjoint=True)),
        ]),
        name=name)

tfd.Independent বন্টন করিয়া স্বাধীন এক বিতরণের স্বপক্ষে, পরিসংখ্যানগত স্বাধীন স্থানাঙ্ক সঙ্গে একটি বহুচলকীয় বন্টন করে। কম্পিউটিং নিরিখে log_prob , ঘটনা মাত্রা (গুলি) উপর একটি সহজ যোগফল হিসাবে এই "মেটা-ডিস্ট্রিবিউশান" টেপা।

এছাড়াও লক্ষ্য করা যে আমরা গ্রহণ adjoint স্কেল ম্যাট্রিক্স ( "TRANSPOSE")। এর কারণ স্পষ্টতা বিপরীত সহভেদাংক, অর্থাত্, যদি \(P=C^{-1}\) এবং যদি \(C=AA^\top\), তারপর \(P=BB^{\top}\) যেখানে \(B=A^{-\top}\)।

যেহেতু এই বন্টন চতুর ধরনের, এর দ্রুত যাচাই করুন যে আমাদের দিন MVNCholPrecisionTriL কাজ করে যেমন আমরা মনে করি এটি করা উচিত।

def compute_sample_stats(d, seed=42, n=int(1e6)):
  x = d.sample(n, seed=seed)
  sample_mean = tf.reduce_mean(x, axis=0, keepdims=True)
  s = x - sample_mean
  sample_cov = tf.linalg.matmul(s, s, adjoint_a=True) / tf.cast(n, s.dtype)
  sample_scale = tf.linalg.cholesky(sample_cov)
  sample_mean = sample_mean[0]
  return [
      sample_mean,
      sample_cov,
      sample_scale,
  ]

dtype = np.float32
true_loc = np.array([1., -1.], dtype=dtype)
true_chol_precision = np.array([[1., 0.],
                                [2., 8.]],
                               dtype=dtype)
true_precision = np.matmul(true_chol_precision, true_chol_precision.T)
true_cov = np.linalg.inv(true_precision)

d = MVNCholPrecisionTriL(
    loc=true_loc,
    chol_precision_tril=true_chol_precision)

[sample_mean, sample_cov, sample_scale] = [
    t.numpy() for t in compute_sample_stats(d)]

print('true mean:', true_loc)
print('sample mean:', sample_mean)
print('true cov:\n', true_cov)
print('sample cov:\n', sample_cov)

true mean: [ 1. -1.]
sample mean: [ 1.0002806 -1.000105 ]
true cov:
 [[ 1.0625   -0.03125 ]
 [-0.03125   0.015625]]
sample cov:
 [[ 1.0641273  -0.03126175]
 [-0.03126175  0.01559312]]

যেহেতু নমুনা গড় এবং সহভঙ্গি প্রকৃত গড় এবং সহভরিতার কাছাকাছি, তাই মনে হচ্ছে বিতরণটি সঠিকভাবে বাস্তবায়িত হয়েছে৷ এখন আমরা ব্যবহার করব MVNCholPrecisionTriL tfp.distributions.JointDistributionNamed BGMM মডেল নির্দিষ্ট করতে। পর্যবেক্ষণ মডেল জন্য, আমরা ব্যবহার করব tfd.MixtureSameFamily স্বয়ংক্রিয়ভাবে আউট সংহত করতে \(\{Z_i\}_{i=1}^N\) স্বপক্ষে।

dtype = np.float64
dims = 2
components = 3
num_samples = 1000

bgmm = tfd.JointDistributionNamed(dict(
  mix_probs=tfd.Dirichlet(
    concentration=np.ones(components, dtype) / 10.),
  loc=tfd.Independent(
    tfd.Normal(
        loc=np.stack([
            -np.ones(dims, dtype),
            np.zeros(dims, dtype),
            np.ones(dims, dtype),
        ]),
        scale=tf.ones([components, dims], dtype)),
    reinterpreted_batch_ndims=2),
  precision=tfd.Independent(
    tfd.WishartTriL(
        df=5,
        scale_tril=np.stack([np.eye(dims, dtype=dtype)]*components),
        input_output_cholesky=True),
    reinterpreted_batch_ndims=1),
  s=lambda mix_probs, loc, precision: tfd.Sample(tfd.MixtureSameFamily(
      mixture_distribution=tfd.Categorical(probs=mix_probs),
      components_distribution=MVNCholPrecisionTriL(
          loc=loc,
          chol_precision_tril=precision)),
      sample_shape=num_samples)
))

def joint_log_prob(observations, mix_probs, loc, chol_precision):
  """BGMM with priors: loc=Normal, precision=Inverse-Wishart, mix=Dirichlet.

  Args:
    observations: `[n, d]`-shaped `Tensor` representing Bayesian Gaussian
      Mixture model draws. Each sample is a length-`d` vector.
    mix_probs: `[K]`-shaped `Tensor` representing random draw from
      `Dirichlet` prior.
    loc: `[K, d]`-shaped `Tensor` representing the location parameter of the
      `K` components.
    chol_precision: `[K, d, d]`-shaped `Tensor` representing `K` lower
      triangular `cholesky(Precision)` matrices, each being sampled from
      a Wishart distribution.

  Returns:
    log_prob: `Tensor` representing joint log-density over all inputs.
  """
  return bgmm.log_prob(
      mix_probs=mix_probs, loc=loc, precision=chol_precision, s=observations)

"প্রশিক্ষণ" ডেটা তৈরি করুন

এই ডেমোর জন্য, আমরা কিছু র্যান্ডম ডেটার নমুনা দেব।

true_loc = np.array([[-2., -2],
                     [0, 0],
                     [2, 2]], dtype)
random = np.random.RandomState(seed=43)

true_hidden_component = random.randint(0, components, num_samples)
observations = (true_loc[true_hidden_component] +
                random.randn(num_samples, dims).astype(dtype))

HMC ব্যবহার করে Bayesian Inference

এখন যেহেতু আমরা আমাদের মডেল নির্দিষ্ট করতে TFD ব্যবহার করেছি এবং কিছু পর্যবেক্ষিত ডেটা পেয়েছি, আমাদের কাছে HMC চালানোর জন্য প্রয়োজনীয় সমস্ত অংশ রয়েছে।

এই কাজের জন্য, আমরা একটি ব্যবহার করব আংশিক আবেদন "নিচে পিন" জিনিস আমরা নমুনা চাই না করতে। এই ক্ষেত্রে উপায়ে আমরা কেবল নিচে পিন প্রয়োজন যে observations । (অধি-পরামিতি ইতিমধ্যে পূর্বে ডিস্ট্রিবিউশন এবং অংশ নয় করার জন্য বেকড হয় joint_log_prob ফাংশন স্বাক্ষর।)

unnormalized_posterior_log_prob = functools.partial(joint_log_prob, observations)

initial_state = [
    tf.fill([components],
            value=np.array(1. / components, dtype),
            name='mix_probs'),
    tf.constant(np.array([[-2., -2],
                          [0, 0],
                          [2, 2]], dtype),
                name='loc'),
    tf.linalg.eye(dims, batch_shape=[components], dtype=dtype, name='chol_precision'),
]

সীমাহীন প্রতিনিধিত্ব

হ্যামিলটোনিয়ান মন্টে কার্লো (HMC) এর আর্গুমেন্টের সাপেক্ষে টার্গেট লগ-সম্ভাব্যতা ফাংশন আলাদা হতে হবে। তদ্ব্যতীত, HMC নাটকীয়ভাবে উচ্চ পরিসংখ্যানগত দক্ষতা প্রদর্শন করতে পারে যদি রাষ্ট্র-স্থান সীমাবদ্ধ না হয়।

এর মানে হল BGMM পোস্টেরিয়র থেকে নমুনা নেওয়ার সময় আমাদের দুটি প্রধান সমস্যা নিয়ে কাজ করতে হবে:

  1. \(\theta\) অর্থাত্, একটি বিযুক্ত সম্ভাব্যতা ভেক্টর প্রতিনিধিত্ব করে, যেমন হওয়া আবশ্যক যে \(\sum_{k=1}^K \theta_k = 1\) এবং \(\theta_k>0\)।
  2. \(T_k\) প্রতিনিধিত্ব করে একটি বিপরীত সহভেদাংক ম্যাট্রিক্স, অর্থাত, যেমন যে হতে হবে \(T_k \succ 0\), অর্থাত্, হয় ইতিবাচক নির্দিষ্ট ।

এই প্রয়োজনীয়তা মোকাবেলা করার জন্য আমাদের প্রয়োজন হবে:

  1. সীমাবদ্ধ ভেরিয়েবলগুলিকে একটি অনিয়ন্ত্রিত স্থানে রূপান্তর করুন
  2. অনিয়ন্ত্রিত স্থানে MCMC চালান
  3. অনিয়ন্ত্রিত ভেরিয়েবলগুলিকে সীমাবদ্ধ স্থানে ফিরিয়ে আনুন।

সঙ্গে MVNCholPrecisionTriL , আমরা ব্যবহার করব Bijector গুলি সংকোচহীন স্থান থেকে র্যান্ডম ভেরিয়েবল রুপান্তর।

  • Dirichlet মাধ্যমে সংকোচহীন স্থান থেকে রুপান্তরিত করা হয় softmax ফাংশন ।

  • আমাদের নির্ভুলতা র্যান্ডম ভেরিয়েবল হল পোস্টিভ সেমিডেফিনিট ম্যাট্রিক্সের উপর একটি বন্টন। এই unconstrain আমরা ব্যবহার করব FillTriangular এবং TransformDiagonal bijectors। এগুলো ভেক্টরকে নিম্ন-ত্রিভুজাকার ম্যাট্রিসে রূপান্তর করে এবং নিশ্চিত করে যে তির্যকটি ধনাত্মক। কারণ এটি সক্ষম শুধুমাত্র স্যাম্পলিং সাবেক দরকারী \(d(d+1)/2\) বদলে floats \(d^2\)।

unconstraining_bijectors = [
    tfb.SoftmaxCentered(),
    tfb.Identity(),
    tfb.Chain([
        tfb.TransformDiagonal(tfb.Softplus()),
        tfb.FillTriangular(),
    ])]

@tf.function(autograph=False)
def sample():
  return tfp.mcmc.sample_chain(
    num_results=2000,
    num_burnin_steps=500,
    current_state=initial_state,
    kernel=tfp.mcmc.SimpleStepSizeAdaptation(
        tfp.mcmc.TransformedTransitionKernel(
            inner_kernel=tfp.mcmc.HamiltonianMonteCarlo(
                target_log_prob_fn=unnormalized_posterior_log_prob,
                 step_size=0.065,
                 num_leapfrog_steps=5),
            bijector=unconstraining_bijectors),
         num_adaptation_steps=400),
    trace_fn=lambda _, pkr: pkr.inner_results.inner_results.is_accepted)

[mix_probs, loc, chol_precision], is_accepted = sample()

আমরা এখন চেইনটি কার্যকর করব এবং পোস্টেরিয়র মানে মুদ্রণ করব।

acceptance_rate = tf.reduce_mean(tf.cast(is_accepted, dtype=tf.float32)).numpy()
mean_mix_probs = tf.reduce_mean(mix_probs, axis=0).numpy()
mean_loc = tf.reduce_mean(loc, axis=0).numpy()
mean_chol_precision = tf.reduce_mean(chol_precision, axis=0).numpy()
precision = tf.linalg.matmul(chol_precision, chol_precision, transpose_b=True)

print('acceptance_rate:', acceptance_rate)
print('avg mix probs:', mean_mix_probs)
print('avg loc:\n', mean_loc)
print('avg chol(precision):\n', mean_chol_precision)

acceptance_rate: 0.5305
avg mix probs: [0.25248723 0.60729516 0.1402176 ]
avg loc:
 [[-1.96466753 -2.12047249]
 [ 0.27628865  0.22944732]
 [ 2.06461244  2.54216122]]
avg chol(precision):
 [[[ 1.05105032  0.        ]
  [ 0.12699955  1.06553113]]

 [[ 0.76058015  0.        ]
  [-0.50332767  0.77947431]]

 [[ 1.22770457  0.        ]
  [ 0.70670027  1.50914164]]]

loc_ = loc.numpy()
ax = sns.kdeplot(loc_[:,0,0], loc_[:,0,1], shade=True, shade_lowest=False)
ax = sns.kdeplot(loc_[:,1,0], loc_[:,1,1], shade=True, shade_lowest=False)
ax = sns.kdeplot(loc_[:,2,0], loc_[:,2,1], shade=True, shade_lowest=False)
plt.title('KDE of loc draws');

png

উপসংহার

এই সাধারণ কোল্যাব দেখিয়েছে কিভাবে TensorFlow সম্ভাব্যতা আদিম শ্রেণীবিন্যাস বায়েসিয়ান মিশ্রণের মডেল তৈরি করতে ব্যবহার করা যেতে পারে।