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')