clearvars;
if ~exist('first','var')
run('custom_settings')
load([getenv('datadir') filesep 'simulated'])
first = 1;
end
K = 1;
[N W] = size(X);
% find the weakest spectra and use them as background estimation
Xpower = sum(X'.^2);
[val inx] = sort(Xpower);
B0 = median(X(inx(1:50),:));
b0 = ones(1,N);
M = 10000;
Mest = 1000;
tic
rng('default');
[Am,Vm,bm,Bm,taum,loglikm,mhjumps,logpm] = gibbsvoigt(X',K,M, ...
'alpha',1e1, ...
'muc',gendata.theta(:,1), ...
'tauc',10*ones(K,1), ...
'mug',10*ones(K,1), ...
'taug',10*ones(K,1), ...
'V',[gendata.theta(:,1)+20*randn(K,1) 5+10*rand(K,1) rand(K,1)], ...
'tau',1, ...
'taub', 1000, ...
'tauB', 100, ...
'b',b0, ...
'B',B0, ...
'nV',false(opt.gendata.K,1), ...
'nA',false(opt.gendata.K,1), ...
'nb',false, ...
'nB',false, ...
'mhstep',[0.4 0.2 0.01], ...
'nt',false);
toc
save([strMATdir 'icassp2017' filesep 'simulated_sbr_4_voigt.mat'] , ...
'Am','Vm','bm','Bm','taum','loglikm','mhjumps','logpm')