<jupyter_start><jupyter_text>Posterior Sampling for Gaussian Mixture Model with CRP using Gibbs samplerReference: https://pdfs.semanticscholar.org/9ece/0336316d78837076ef048f3d07e953e38072.pdf<jupyter_code>import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
# Lets generate some data
X1 = np.random.multivariate_normal([5, 5], np.diag([0.5, 0.5]), size=20)
X2 = np.random.multivariate_normal([8, 8], np.diag([0.5, 0.5]), size=20)
X3 = np.random.multivariate_normal([20, 20], np.diag([0.5, 0.5]), size=10)

X = np.vstack([X1, X2, X3])

X
N, D = X.shape
# GMM paramters
mus = [] # This lists a 2x1 vector (mean vector of each gaussian)
sigma = np.eye(D)
precision = np.linalg.inv(sigma)
zs = np.zeros([N], dtype=int)
C = [] # Cluster, bu=inary matrix of KxM
Ns = [] # Count of each cluster
sigma
precision
zs
# CRP prior
alpha = 200
# Base distribution prior: N(mutheta, precisiontheta)
mu_theta = np.ones(D)
sigma_theta = np.eye(D)
precision_theta = np.linalg.inv(np.eye(D))
G_theta = stats.multivariate_normal(mean=mu_theta, cov=np.eye(D))
# Initiate the ONE cluster
C.append(np.ones(N, dtype=int))
zs[:] = 0
Ns.append(N)
mus.append(G_theta.rvs())

K = 1
mvn = stats.multivariate_normal
# Gibbs sampler
for it in range(20):
    # here we sample from full conditional of assignment from CRP prior
    # z ~ GEM(alpha)
    # Now, for each data point, draw a cluster assignment
    for i in range(N):
         # Remove assignment from cluster
        # ------------------------------

        zi = zs[i]
        C[zi][i] = 0
        Ns[zi] -= 1

        # If empty, remove cluster
        if Ns[zi] == 0:
            # Fix indices
            zs[zs > zi] -= 1

            # Delete cluster
            del C[zi]
            del Ns[zi]
            del mus[zi]

            # Decrement cluster count
            K -= 1

        # Draw new assignment zi weighted by CRP prior
        # --------------------------------------------

        probs = np.zeros(K + 1)
        zs_minus_i = zs[np.arange(len(zs)) != i]

        # Probs of joining existing cluster
        for k in range(K):
            nk_minus = zs_minus_i[zs_minus_i == k].shape[0]
            crp = nk_minus / (N + alpha - 1)
            probs[k] = crp * mvn.pdf(X[i], mus[k], sigma)

        # Prob of creating new cluster
        crp = alpha / (N + alpha - 1)
        lik = mvn.pdf(X[i], mu_theta, sigma_theta + sigma)  # marginal dist. of x
        probs[K] = crp*lik

        # Normalize
        probs /= np.sum(probs)

        # Draw new assignment for i
        z = np.random.multinomial(n=1, pvals=probs).argmax()

        # Update assignment trackers
        if z == K:
            C.append(np.zeros(N, dtype=int))
            Ns.append(0)
            mus.append(G_theta.rvs())
            K += 1

        zs[i] = z
        C[z][i] = 1
        Ns[z] += 1

    # -------------------------------------------------
    # Sample from full conditional of cluster parameter
    # -------------------------------------------------

    # Assume fixed covariance => posterior is Normal
    # mu ~ N(mu, sigma)
    for k in range(K):
        Xk = X[zs == k]
        Ns[k] = Xk.shape[0]

        # Covariance of posterior
        lambda_post = precision_theta + Ns[k]*precision
        cov_post = np.linalg.inv(lambda_post)

        # Mean of posterior
        left = cov_post
        right = precision_theta @ mu_theta + Ns[k]*precision @ np.mean(Xk, axis=0)
        mus_post = left @ right

        # Draw new mean sample from posterior
        mus[k] = mvn.rvs(mus_post, cov_post)<jupyter_output><empty_output><jupyter_text>Even though we only initialize with one cluster, the result should be: **Expected output:**``` python 20 data in cluster-0, mean: [ 5  5 ]20 data in cluster-1, mean: [ 8  8 ]10 data in cluster-2, mean: [ 20  20 ]``` Note: cluster label is exchangeable<jupyter_code>for k in range(K):
    print('{} data in cluster-{}, mean: {}'.format(Ns[k], k, mus[k]))<jupyter_output>19 data in cluster-0, mean: [ 4.5457113   4.93233594]
21 data in cluster-1, mean: [ 8.25548463  7.50082634]
10 data in cluster-2, mean: [ 18.60459256  18.30911234]
{'fraction_non_alphanumeric': 0.16175760502172865, 'fraction_numerical': 0.04056011588604539, 'mean_word_length': 2.1338880484114977, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 10, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '382278', 'license': 'Apache-2.0', 'n_tokens_mistral': 1554, 'n_tokens_neox': 1401, 'n_words': 399, 'path': 'deep-learning/Models_Algorithms_More/Bayesian Non-parametric Part 2.ipynb', 'repo_name': 'HiteshDhola/Artificial-Intelligence-Deep-Learning-Machine-Learning-Tutorials', 'size': '4142'}