function [exp_r,xp,r_samp,g_post] = spm_BMS_gibbs (lme, alpha0, Nsamp)
% Bayesian model selection for group studies using Gibbs sampling
% FORMAT [exp_r,xp,r_samp,g_post] = spm_BMS_gibbs (lme, alpha0, Nsamp)
%
% INPUT:
% lme      - array of log model evidences 
%              rows: subjects
%              columns: models (1..Nk)
% alpha0   - [1 x Nk] vector of prior model counts
% Nsamp    - number of samples (default: 1e6)
% 
% OUTPUT:
% exp_r    - [1 x  Nk] expectation of the posterior p(r|y)
% xp       - exceedance probabilities
% r_samp   - [Nsamp x Nk] matrix of samples from posterior
% g_post   - [Ni x Nk] matrix of posterior probabilities with 
%            g_post(i,k) being post prob that subj i used model k
%__________________________________________________________________________
% Copyright (C) 2009-2013 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_BMS_gibbs.m 6381 2015-03-17 17:55:09Z will $


if nargin < 3 || isempty(Nsamp)
    Nsamp = 1e4;
end

Ni = size(lme,1);  % number of subjects
Nk = size(lme,2);  % number of models

% prior observations
%--------------------------------------------------------------------------
if nargin < 2 || isempty(alpha0)
    alpha0 = ones(1,Nk);    
end
alpha0     = alpha0(:)';

% Initialise; sample r from prior
r  = zeros(1,Nk);
for k = 1:Nk
    r(:,k) = spm_gamrnd(alpha0(k),1);
end
sr = sum(r,2);
for k = 1:Nk
    r(:,k) = r(:,k)./sr;
end
        
% Subtract max evidence for subject 
lme = lme - max(lme,[],2)*ones(1,Nk);

% Gibbs sampling 
r_samp = zeros(Nsamp,Nk);
g_post = zeros(Ni,Nk);

for samp = 1:2*Nsamp
    
    mod_vec = sparse(Ni,Nk);
    % Sample m's given y, r
    for i = 1:Ni
        % Pick a model for this subject
        u         = exp(lme(i,:) + log(r)) + eps;
        g         = u / sum(u);
        gmat(i,:) = g;
        modnum    = spm_multrnd(g,1);
        mod_vec(i,modnum) = 1;
    end
    
    % Sample r's given y, m
    beta          = sum(mod_vec,1);
    alpha         = alpha0+beta;
    for k = 1:Nk
        r(:,k)    = spm_gamrnd(alpha(k),1);
    end
    sr = sum(r,2);
    for k = 1:Nk
        r(:,k)    = r(:,k) ./ sr;
    end

    % Only keep last Nsamp samples
    if samp > Nsamp
        r_samp(samp-Nsamp,:) = r;
        g_post    = g_post+gmat;
    end
    
    if mod(samp,1e4)==0
        fprintf('%d samples out of %d\n',samp,2*Nsamp);
    end
    
end
g_post = g_post/Nsamp;

% Posterior mean
exp_r = mean(r_samp,1);

% Exceedence probs
xp    = zeros(1,Nk);
[y,j] = max(r_samp,[],2);
tmp   = histc(j,1:Nk)';
xp    = tmp / Nsamp;