123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100 |
- 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;
|