<jupyter_start><jupyter_text>MPPCA demoMixtures of Probabilistic Principal Component Analysers is a simple yet powerful algorithm used to cluster data into linear subspaces. Its applications cover clustering, density estimation and classification.Reference : "Mixtures of Probabilistic Principal Component Analysers", <NAME> and <NAME>, Neural Computation 11(2), pp 443–482, MIT Press, 1999<jupyter_code>import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
import pickle
import gzip
import time
import collections

import mppca

%matplotlib inline
def print_clusters(X, clusters, mu):
    N, d = X.shape
    
    if d == 2:
        for i in range(max(clusters)+1):
            plt.scatter(X[clusters==i, 0], X[clusters==i, 1], c=np.random.rand(3,1))
            
        plt.scatter(mu[:, 0], mu[:, 1], c='k', s=100)
            
    elif d==3:
        fig = plt.figure()
        ax = Axes3D(fig)

        for i in range(max(clusters)+1):
            ax.scatter(X[clusters==i, 0], X[clusters==i, 1], X[clusters==i, 2], c=np.random.rand(3,1))
        # ax.scatter(mu[:, 0], mu[:, 1], mu[:, 2], c='k', s=100) does not work
        
    else:
        raise ValueError('Incorrect data dimensionality')<jupyter_output><empty_output><jupyter_text>Clustering on a cross  Construct the dataset<jupyter_code>def dataset_cross(N, d, p):
    """
    N : number of desired points
    d : space dimension
    p : proportion of the 1st branch of the cross (0<p<1)

    X : dataset
    v1, v2 : directions of the branches
    center : center of the cross
    """

    t = np.random.randn(N, 1)
    v1 = np.random.randn(1, d)
    v2 = 1.5*np.random.randn(1, d)
    center = 5.*np.random.randn(1, d)

    separation = int(round(N*p))

    X = np.zeros((N, d))
    X[:separation, :] = v1*t[:separation, :]
    X[separation:, :] = v2*t[separation:, :]
    X += 0.1*np.random.randn(N, d)
    X += center

    return X, v1, v2, center
X, v1, v2, center = dataset_cross(N=500, d=3, p=0.4)

fig = plt.figure()
ax = Axes3D(fig)

ax.scatter(X[:,0], X[:,1], X[:,2])<jupyter_output><empty_output><jupyter_text>Parameter initialization with a K-means<jupyter_code>p = 2 # We expect 2 clusters
q = 1 # We expect 1-dimensional subspaces
pi, mu, W, sigma2, clusters = mppca.initialization_kmeans(X, p, q)
print_clusters(X, clusters, mu)<jupyter_output><empty_output><jupyter_text>K-means cannot cluster correctly a cross  Run the MPPCA algotrithm<jupyter_code>niter = 500
pi, mu, W, sigma2, R, L, sigma2hist = mppca.mppca_gem(X, pi, mu, W, sigma2, niter)<jupyter_output>....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................<jupyter_text>Print the clusters<jupyter_code>clusters = R.argmax(axis=1)
print_clusters(X, clusters, mu)
plt.plot(L)<jupyter_output><empty_output><jupyter_text>Loss curve  Decomposition of a swiss roll into linear subspaces  Construct the dataset<jupyter_code>def dataset_suissroll(N, d):
    assert (d==2)

    t = np.random.randn(N, 2)
    X = np.zeros((N, 2))
    X[:, 0] = ((t[:, 0]+4)*3+t[:, 1])*np.cos((t[:, 0]+4)*1.5)
    X[:, 1] = ((t[:, 0]+4)*3+t[:, 1])*np.sin((t[:, 0]+4)*1.5)

    return X
X = dataset_suissroll(N=1000, d=2)

plt.scatter(X[:,0], X[:,1])<jupyter_output><empty_output><jupyter_text>Parameter initialization with a k-means<jupyter_code>p = 7 # We expect 7 clusters
q = 1 # We expect 1-dimensional subspaces
pi, mu, W, sigma2, clusters = mppca.initialization_kmeans(X, p, q)
print_clusters(X, clusters, mu)<jupyter_output><empty_output><jupyter_text>Run the MPPCA algotrithm<jupyter_code>niter = 500
pi, mu, W, sigma2, R, L, sigma2hist = mppca.mppca_gem(X, pi, mu, W, sigma2, niter)<jupyter_output>....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................<jupyter_text>Print the clusters<jupyter_code>clusters = R.argmax(axis=1)
print_clusters(X, clusters, mu)
plt.plot(L)<jupyter_output><empty_output><jupyter_text>Classification on the MNIST dataset  Load the dataset<jupyter_code>f = gzip.open('mnist.pkl.gz', 'r')
train_set, valid_set, test_set = pickle.load(f, encoding='latin1')
f.close()

train_data, train_labels = train_set
valid_data, valid_labels = valid_set
test_data, test_labels = test_set

d = train_data.shape[1]
K = 10 # Number of different numbers
len(train_set[0]), len(valid_set[0]), len(test_set[0])
i = 52
plt.imshow(train_data[i].reshape((28, 28)))
print('label:%d'%train_labels[i])<jupyter_output>label:7<jupyter_text>Train the model<jupyter_code>n_train = 100  # Faster computations at the expense of results

# Optimal p, q values depend on n_train.
# NB: On real applications, these hyperparameters should be trained
# using the validation set!
p = 3 # Number of expected clusters (per number)
q = 2 # Dimension of the subspaces

niter = 100

pi = np.zeros(p*K)
mu = np.zeros((p*K, d))
W = np.zeros((p*K, d, q))
sigma2 = np.zeros(p*K)
cluster_labels = np.zeros(p*K, dtype=np.int32)

for i in range(K):
    X = train_data[train_labels == i][:n_train, :]

    pi_i, mu_i, W_i, sigma2_i, clusters = mppca.initialization_kmeans(X, p, q)

    t = time.time()
    pi_i, mu_i, W_i, sigma2_i, R, L, sigma2hist = mppca.mppca_gem(
        X, pi_i, mu_i, W_i, sigma2_i, niter)
    print('GEM took %f seconds.'%(time.time()-t))
    
    pi[p*i:p*(i+1)] = pi_i
    mu[p*i:p*(i+1), :] = mu_i
    W[p*i:p*(i+1), :, :] = W_i
    sigma2[p*i:p*(i+1)] = sigma2_i
    cluster_labels[p*i:p*(i+1)] = i<jupyter_output>....................................................................................................GEM took 21.422765 seconds.
....................................................................................................GEM took 17.008500 seconds.
....................................................................................................GEM took 18.368042 seconds.
....................................................................................................GEM took 16.966522 seconds.
....................................................................................................GEM took 18.250899 seconds.
....................................................................................................GEM took 23.181151 seconds.
....................................................................................................GEM took 20.651325 seconds.
....................................................................................................GEM [...]<jupyter_text>Visualize some parameters<jupyter_code># cluster repartition
h = plt.hist(pi, bins=100)
i = 12
# center of the i-th cluster
plt.imshow(mu[i, :].reshape((28, 28)))
print('label=%d'%cluster_labels[i])
# Plot the 1st principal component of the i-th cluster
plt.imshow(W[i, :, 0].reshape((28, 28)))<jupyter_output><empty_output><jupyter_text>Predict on the test set<jupyter_code>R = mppca.mppca_predict(test_data, pi, mu, W, sigma2)
clusters = R.argmax(axis=1)
predicted_numbers = cluster_labels[clusters]
error_rate = (test_labels!=predicted_numbers).sum()/len(test_labels)
print(error_rate)<jupyter_output>0.1323
{'fraction_non_alphanumeric': 0.3429945576509303, 'fraction_numerical': 0.030122769269712694, 'mean_word_length': 1.7285911602209945, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 50, 'https://': 0, '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': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '302429', 'license': 'WTFPL', 'n_tokens_mistral': 2608, 'n_tokens_neox': 2435, 'n_words': 581, 'path': 'mppca_demo.ipynb', 'repo_name': 'michelbl/MPPCA', 'size': '7928'}