Clustering with Dirichlet process mixtures

Practical session at the Applied Bayesian Statistics School, Como, June 2014.

In this course we will consider Dirichlet process mixture of Gaussians with a conjugate normal-inverse Wishart base distribution.

The practical is divided into five parts:

Part 1: Investigate the properties of the Chinese restaurant process partition through simulations

Part 2: Investigate the properties of the Pitman-Yor CRP partition through simulations

Part 3: Simulation from the Dirichlet process mixture of Gaussians

Part 4: Posterior inference via Gibbs sampling for BNP clustering on simulated data. Four algorithms will be considered: (a) Gibbs sampler based on the Blackwell-MacQueen urn (b) Gibbs sampler based on the Chinese Restaurant process (c) Collapsed Gibbs sampler based on the Chinese Restaurant process (d) Slice sampler

Part 5: BNP clustering of gene expression data

Author: François Caron, University of Oxford

Tested on Matlab R2014a and Octave 6.4.2.

Contents

close all
clear all

PART 1: Properties of the Chinese Restaurant Process random partition

This part will make use of the function crprnd, which generates a partition from the Chinese Restaurant process. We will also use the function dpstickrnd that generates weights from the stick-breaking construction of the Dirichlet Process. Have a look to these functions.

type('crprnd')
type('dpstickrnd')

function [partition, m, K, partition_bin] = crprnd(alpha, n)

%
% CRPRND samples a partition from the Chinese restaurant process
%   [partition, m, K, partition_bin] = crprnd(alpha, n)
%
%--------------------------------------------------------------------------
% INPUTS
%   - alpha:    scale parameter of the CRP
%   - n:        number of objects in the partition
%
% OUTPUTS
%   - partition: Vector of length n. partition(i) is the cluster membership
%               of object i. Clusters are numbered by order of appearance
%   - m:        Vector of length n. m(j) is the size of cluster j.
%               m(j)=0 for j>K
%   - K:        Integer. Number of clusters
%   - partition_bin: locical matrix of size n*n. partition(i,j)=true if object
%               i is in cluster j, false otherwise
%--------------------------------------------------------------------------
% EXAMPLE
% alpha = 3; n= 100;
% [partition, m, K, partition_bin] = crprnd(alpha, n);
%--------------------------------------------------------------------------


% Check parameter
if alpha<=0
    error('Parameter alpha must be a positive scalar')
end

m = zeros(1, n);
partition = zeros(1, n);
partition_bin = false(n);

% Initialization
partition(1) = 1;
partition_bin(1,1) = true;
m(1) = 1;
K = 1;
% Iterations
for i=2:n
    % Compute the probability of joining an existing cluster or a new one
    proba = [m(1:K), alpha]/(alpha+i-1);
    % Sample from a discrete distribution w.p. proba
    u = rand;
    partition(i) = find(u<=cumsum(proba), 1);
    partition_bin(i,partition(i)) = true;
    % Increment the size of the cluster
    m(partition(i)) = m(partition(i)) + 1;
    % Increment the number of clusters if new
    K = K + isequal(partition(i), K+1);    
end

function weights = dpstickrnd(alpha, K)

weights = zeros(K, 1);
weights(1) = betarnd(1, alpha);
for i=2:K
    b = betarnd(1, alpha);
    weights(i) = b * (1.0 - sum(weights(1:i)));
    if sum(weights)>=1.0        
        break;
    end
end

PART 2: Properties of the Pitman Yor random partition

This part will make use of the function pycrprnd, which generates a partition from the Pitman-Yor Chinese Restaurant process. Have a look at this function.

type('pycrprnd')

function [partition, m, K, partition_bin] = pycrprnd(alpha, sigma, n)

%
% PYCRPRND samples a partition from the Pitman-Yor Chinese restaurant process
%   [partition, m, K, partition_bin] = pycrprnd(alpha, sigma, n)
%
%--------------------------------------------------------------------------
% INPUTS
%   - alpha:    scale parameter of the PY CRP
%   - sigma:    discount parameter of the PY CRP
%   - n:        number of objects in the partition
%
% OUTPUTS
%   - partition: Vector of length n. partition(i) is the cluster membership
%               of object i. Clusters are numbered by order of appearance
%   - m:        Vector of length n. m(j) is the size of cluster j.
%               m(j)=0 for j>K
%   - K:        Integer. Number of clusters
%   - partition_bin: locical matrix of size n*n. partition(i,j)=true if object
%               i is in cluster j, false otherwise
%--------------------------------------------------------------------------
% EXAMPLE
% alpha = 3; n= 100;
% [partition, m, K, partition_bin] = crprnd(alpha, n);
%--------------------------------------------------------------------------


m = zeros(n, 1);
partition = zeros(n, 1);
partition_bin = false(n);

% Initialization
partition(1) = 1;
partition_bin(1,1) = true;
m(1) = 1;
K = 1;
% Iterations
for i=2:n
    % Compute the probability of joining an existing cluster or a new one
    proba = [m(1:K)-sigma; alpha + sigma*K]/(alpha+i-1);
    % Sample from a discrete distribution w.p. proba
    u = rand;
    partition(i) = find(u<=cumsum(proba), 1);
    partition_bin(i,partition(i)) = true;
    % Increment the size of the cluster
    m(partition(i)) = m(partition(i)) + 1;
    % Increment the number of clusters if new
    K = K + isequal(partition(i), K+1);    
end
partition_bin = partition_bin(:, 1:K);

PART 3: Simulation from a Dirichlet process mixture of Gaussians

Create a function dpmgaussrnd that simulates data from a Dirichlet process mixture of Gaussians where the base distribution is normal inverse Wishart. Use the function normalinvwishrnd.m to simulate from the normal inverse Wishart distribution. Plot samples from the DPM of Gaussians, for various values of alpha and of the parameters of the normal inverse Wishart distribution.

PART 4: Clustering with Dirichlet process mixtures of Gaussians on simulated data

In this part, we'll consider various MCMC algorithms to get posterior estimates of the partition.

Load and plot the simulated data

load ./data/simudata
figure('name', 'simulated data')
plot(y(1,:), y(2,:), '.')
xlabel('X')
ylabel('Y')

% Parameters of the base distribution G0
hyperG0.mu = [0;0];
hyperG0.kappa = 1;
hyperG0.nu = 4;
hyperG0.lambda = eye(2);
% Scale parameter of DPM
alpha = 3;
% Number of iterations
niter = 20;
% do some plots
doPlot = 2;
type_algo = 'BMQ'; % other algorithms: 'CRP', 'collapsesCRP', 'slicesampler'

[c_st, c_est, similarity] = gibbsDPM(y, hyperG0, alpha, niter, type_algo, doPlot);

Iteration 2/20
10 clusters

Iteration 3/20
14 clusters

Iteration 4/20
9 clusters

Iteration 5/20
9 clusters

Iteration 6/20
9 clusters

Iteration 7/20
9 clusters

Iteration 8/20
13 clusters

Iteration 9/20
15 clusters

Iteration 10/20
11 clusters

Iteration 11/20
10 clusters

Iteration 12/20
12 clusters

Iteration 13/20
8 clusters

Iteration 14/20
5 clusters

Iteration 15/20
5 clusters

Iteration 16/20
6 clusters

Iteration 17/20
9 clusters

Iteration 18/20
5 clusters

Iteration 19/20
5 clusters

Iteration 20/20
5 clusters

PART 5: BNP clustering of gene expression data

We now consider BNP clustering on gene expression data. The data are composed of the expression of 384 genes under 5 different experimental conditions. The objective is to cluster those genes in order to identify genes with the same function.

% load and plot gene expression data
load ./data/genedata

figure('name', 'data')
imagesc(y')
xlabel('Experiments')
ylabel('Genes')
colorbar

figure('name', 'data')
plot(y)
xlabel('Experiments')
ylabel('Expression')
set(gca, 'Xtick', 1:5)

Published with MATLAB® R2013a