Bayesian nonparametric latent feature models: Beta process, Indian Buffet process and related models

This Matlab/Octave script reproduces the figures of the lecture on BNP latent features models given at the Applied Bayesian Statistics School, Como, June 2014.

Author: François Caron, University of Oxford

Tested on Matlab R2014a and Octave 6.4.2.

One-parameter Indian Buffet process
Parametric beta-Bernoulli model
Beta Bernoulli process
The three parameter stable Indian buffet process

set(0, 'defaultaxesfontsize', 16)
rng('default')

One-parameter Indian Buffet process

Content of the function sample_ibp.m

type('ibprnd.m')

function Z = ibprnd(alpha, n)

%
% IBPRND samples from the one-parameter Indian buffet process
%   Z = IBPRND(alpha, n)
%
% Requires the statistics toolbox
%--------------------------------------------------------------------------
% INPUT
%   - alpha: scale parameter of the IBP
%   - n:     number of objects
%
% OUTPUT
%   - Z:    logical matrix of size n*K
%--------------------------------------------------------------------------


Jmax = 100000; 
Z = false(n, Jmax);
m = zeros(1, Jmax);

% First object
K = poissrnd(alpha);
Z(1,1:K) = true;
m(1:K) = 1;
for i=2:n
    % Each object picks feature j w.p. m_j/i
    Z(i, 1:K) = (rand(1,K)<m(1:K)/i);
    % New features 
    Knew = poissrnd(alpha/i);
    Z(i, K+1:K+Knew) = true;
    K = K + Knew;    
    % Update the counts
    m(1:K) = m(1:K) + Z(i,1:K);
end
Z = Z(:, 1:K);

n = 50; % Number of objects
alpha_all = [1, 5, 10]; % Scale parameter

for i=1:length(alpha_all)
    % Sample from the IBP
    Z = ibprnd(alpha_all(i), n);

    % Plot the feature allocations
    figure('name', 'IBP')
    imagesc(zeros(n, 60))
    hold on
    imagesc(double(Z))
    colormap('gray')
    ylabel('Objects')
    xlabel('Features')
    title(['alpha=' num2str(alpha_all(i))]);
    xlim([0,60])
    saveas(gca, ['ibp' num2str(alpha_all(i))], 'epsc2')
end

Illustration of some properties of the IBP

alpha = 10;
n = 10000;
Z = ibprnd(alpha, n);

% Distribution of the number of features per objects
figure
[freq, bins] = hist(sum(Z,2), 1:20);
plot(bins, freq, 'o', 'MarkerSize', 5,  'MarkerEdgeColor' , 'none',...
    'MarkerFaceColor' , [1 .6 .6] );
xlabel('Degree of objects')
ylabel('Number of occurences')
saveas(gca, 'ibpdegreeobjects', 'epsc2')

% Distribution of the proportion of objects having a given number of feature
figure
h = plot_degree(Z, 2);
set(h, 'Marker', 'o', 'MarkerSize', 5,  'MarkerEdgeColor' , 'none',...
    'MarkerFaceColor' , [1 .6 .6] );
xlabel('Degree of features')
saveas(gca, 'ibpdegreefeatures', 'epsc2')

% Number of features versus number of objects
figure
h = plot(sum(cumsum(Z, 1)>0, 2));
set(h, 'linewidth', 3, 'color' , [1 .6 .6] );
xlabel('Number of objects')
ylabel('Number of features')
saveas(gca, 'ibpnbfeatures', 'epsc2')

Parametric beta-Bernoulli model

Content of the function betaberrnd.m

type('betaberrnd.m')

function z = betaberrnd(alpha, n, p)

% BETABER samples from a parametric beta-Bernoulli model
%   


pi = betarnd(alpha/p, 1, 1, p);
z = (rand(n, p)<repmat(pi, n, 1));

Simulation

p = 100; alpha = 5; n = 50;
Z = betaberrnd(alpha, n, p);
figure
imagesc(double(Z))
colormap('gray')
ylabel('Objects')
xlabel('Features')
saveas(gca, 'betaber100', 'epsc2')

p = 500;
alpha = 5;
n = 50;
Z = betaberrnd(alpha, n, p);
figure
imagesc(double(Z))
colormap('gray')
ylabel('Objects')
xlabel('Features')
saveas(gca, 'betaber500', 'epsc2')

Beta Bernoulli process

Content of the function ibpstickrnd

type('ibpstickrnd')

function [Z, weights] = ibpstickrnd(alpha, n)

%
% IBPSTICKRND samples the feature allocation and stick weights from the one-parameter Indian buffet process
%   [Z, weights] = ibpstickrnd(alpha, n)
%
% Requires the statistics toolbox
%--------------------------------------------------------------------------
% INPUT
%   - alpha: scale parameter of the IBP
%   - n:     number of objects
%
% OUTPUT
%   - Z:    logical matrix of size n*K
%   - weights: stick weights of each feature
%--------------------------------------------------------------------------


Jmax = 100000; 
Z = false(n, Jmax);
weights = zeros(1, Jmax);

% First object
K = poissrnd(alpha);
Z(1,1:K) = true;
weights(1:K) = betarnd(1, 1, 1, K);
for i=2:n
    % Each object picks feature j w.p. m_j/i
    Z(i, 1:K) = (rand(1,K)<weights(1:K));
    % New features 
    Knew = poissrnd(alpha/i);
    Z(i, K+1:K+Knew) = true;
    
    % Update the stick weights
    weights(K+1:K+Knew) = betarnd(1, i, 1, Knew);
    K = K + Knew;    
    
    
end
Z = Z(:, 1:K);
weights = weights(1:K);

Simulation

alpha = 3;
n = 50;
[Z, weights] = ibpstickrnd(alpha, n);
locations = rand(size(weights));% feature locations are U(0,1)


figure('name', 'Lévy measure')
stepsize = 1e-1;
[A, B] = meshgrid( 0:.01:1, stepsize:stepsize:1);
C = alpha.*B.^(-1).*(1-B).^(alpha-1);
surf(A,B,C)
shading interp
view(31, 26)
hold on
h = stem(locations, weights, 'or');
set(h,'linewidth', 2, 'markersize', 6, 'MarkerFaceColor', 'r');
xlabel('Feature space \Theta')
ylabel('Stick weights')
zlabel('Lévy intensity')
saveas(gca, 'betaprocess', 'epsc2')

figure('name', 'beta Bernoulli process')
subplot(2,1,1)
stem(locations, weights, 'Marker', 'none','color', 'r', 'linewidth', 3);
legend('B', 'location', 'EastOutside')
legend boxoff
ylim([0,1])
ylabel('Stick weights')
subplot(2,1,2)
for i=1:n
    h = plot(locations(Z(i,:)), i*ones(sum(Z(i,:)), 1), 'o');
    set(h, 'Marker', 'o', 'MarkerSize', 5,  'MarkerEdgeColor' , 'none',...
    'MarkerFaceColor' , [1 .6 .6] );
    hold on
end
ylim([.5, n+.5])
xlabel('Feature space \Theta')
ylabel('Objects')
legend('Z', 'location', 'EastOutside')
legend boxoff
saveas(gca, 'betaberprocess', 'epsc2')

The three parameter stable Indian buffet process

Content of the function stableibprnd.m

type('stableibprnd.m')

function Z = stableibprnd(alpha, sigma, c, n)

%
% STABLEIBPRND samples from the three-parameter stable Indian buffet process
% Z = stableibprnd(alpha, sigma, c, n)
%
% Requires the statistics toolbox
%--------------------------------------------------------------------------
% INPUT
%   - alpha: mass parameter of the IBP >=0
%   - sigma: stability exponent in [0,1)
%   - c:     concentration parameter c > -sigma
%   - n:     number of objects
%
% OUTPUT
%   - Z:    logical matrix of size n*K
%--------------------------------------------------------------------------


Jmax = 100000; 
Z = false(n, Jmax);
m = zeros(1, Jmax);

% First object
K = poissrnd(alpha);
Z(1,1:K) = true;
m(1:K) = 1;
for i=2:n
    % Each object picks feature j w.p. (m_j-sigma)/(i-1 + c)
    Z(i, 1:K) = ( rand(1,K) < ((m(1:K)-sigma)/(i-1 + c)) );
    % New features 
    Knew = poissrnd(alpha* exp(gammaln(1 + c) + gammaln(i - 1 + c + sigma)...
        - gammaln(i + c) - gammaln(c+sigma)));
    Z(i, K+1:K+Knew) = true;
    % Update the counts
    m(1:K) = m(1:K) + Z(i, 1:K);
    m(K+1:K+Knew) = 1;
    
    K = K + Knew;     
end
Z = Z(:, 1:K);

Simulations

alpha = 5; % Scale parameter
sigma_all = [0, .5, .9];
c = 1;
n = 40;
for i=1:length(sigma_all)
    % Sample from the stable IBP
    Z = stableibprnd(alpha, sigma_all(i), c, n);

    % Plot the feature allocations
    figure('name', 'stable IBP')
    imagesc(zeros(n, 160))
    hold on
    imagesc(double(Z))
    ylabel('Objects')
    xlabel('Features')
    colormap('gray')
    ylabel('Objects')
    xlabel('Features')
    xlim([0,160])
    title(['sigma=' num2str(sigma_all(i))]);
    saveas(gca, ['stableibp' num2str(10*sigma_all(i))], 'epsc2')
end

Illustration of the properties of the model

n = 1000;
% Parameters of the stable IBP
alpha = 10;
sigma_all = [0, .5, .9];
c = 1;

h1 = figure('name', 'Degree of objects');
h2 = figure('name', 'Degree of features');
h3 = figure('name', 'Nb of features vs nb objects');
colors = [1 .6 .6; .6, 1, .6; .6, .6, 1];
for i=1:length(sigma_all)
    % Sample from the stable IBP
    Z = stableibprnd(alpha, sigma_all(i), c, n);

    figure(h1)
    [freq, bins] = hist(sum(Z,2), 1:20);
    h = plot(bins, freq, '-o') ;
    set(h, 'Marker', 'o', 'MarkerSize', 5,  'MarkerEdgeColor' , 'none',...
        'MarkerFaceColor' , colors(i,:), 'color', colors(i,:));
    hold on

    figure(h2)
    h = plot_degree(Z, 2);
    set(h, 'Marker', 'o', 'MarkerSize', 5,  'MarkerEdgeColor' , 'none',...
        'MarkerFaceColor' , colors(i,:) );
    hold on

    temp = sum(cumsum(Z, 1)>0, 2);
    figure(h3)
    h = plot(sum(cumsum(Z, 1)>0, 2));
    set(h, 'linewidth', 3, 'color' , colors(i,:) );
    hold on
end

Each object has marginally a Poisson number of features, with the same rate

figure(h1)
xlabel('Degree of objects')
ylabel('Number of occurences')
legend('\sigma=0', '\sigma=0.5', '\sigma=0.9')
saveas(gca, 'stableibpdegreeobjects', 'epsc2')

For sigma>0 the distribution of the degree of features asymptotically follows a distribution with a power-law behavior

figure(h2)
xlabel('Degree of features')
ylabel('Distribution')
legend('\sigma=0', '\sigma=0.5', '\sigma=0.9')
ylim([1e-6, 1])
saveas(gca, 'stableibpdegreefeatures', 'epsc2')

The number of features grows in O(log(n)) for sigma=0, and O(n^sigma) for sigma>0.

figure(h3)
xlabel('Number of objects')
ylabel('Number of features')
legend({'\sigma=0', '\sigma=0.5', '\sigma=0.9'}, 'location', 'northwest')
xlim([0, 100])
saveas(gca, 'stableibpnbfeatures', 'epsc2')

Bayesian nonparametric latent feature models: Beta process, Indian Buffet process and related models

Contents

One-parameter Indian Buffet process

Parametric beta-Bernoulli model

Beta Bernoulli process

The three parameter stable Indian buffet process