-
Notifications
You must be signed in to change notification settings - Fork 2
Expand file tree
/
Copy pathRunCalcParEstsAndDICdiffs2.m
More file actions
executable file
·51 lines (46 loc) · 1.76 KB
/
RunCalcParEstsAndDICdiffs2.m
File metadata and controls
executable file
·51 lines (46 loc) · 1.76 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
function DICminMdl=RunCalcParEstsAndDICdiffs2(rslts,np,burnin,IPD,str,plotMltpl,doPlots,savePlots,runName)
% set(0,'DefaultFigureVisible','off') % uncomment to stop figures being plotted to screen
nMdls=numel(rslts);
%% Calculate parameter estimates and CIs for different models
thin=1; % no thinning currently to match DIC calcs
mode_p=NaN(nMdls,np);
HPDI=NaN(np,2,nMdls);
mode_p1=NaN(nMdls,1);
HPDI1=NaN(nMdls,2);
pcorr=cell(nMdls,1);
pIPtAcorr=cell(nMdls,1);
mode_sptl=NaN(nMdls,1);
HPDI_sptl=NaN(nMdls,2);
mode_bckgrnd=NaN(nMdls,1);
HPDI_bckgrnd=NaN(nMdls,2);
mode_d_half=NaN(nMdls,1);
HPDI_d_half=NaN(nMdls,2);
mode_d_half_out=NaN(nMdls,1);
HPDI_d_half_out=NaN(nMdls,2);
mode_WHHRI=NaN(nMdls,1);
HPDI_WHHRI=NaN(nMdls,2);
mean_IP=NaN(nMdls,1);
mode_IP=NaN(nMdls,1);
HPDI_IP=NaN(nMdls,2);
MdlParEsts=NaN(nMdls,3*(np+1));
for i=1:nMdls
fprintf([rslts{i} '\n'])
[mode_p(i,:),HPDI(:,:,i),mode_p1(i),HPDI1(i,:),pcorr{i},pIPtAcorr{i},mode_sptl(i),HPDI_sptl(i,:),mode_bckgrnd(i),HPDI_bckgrnd(i,:),mode_d_half(i),HPDI_d_half(i,:),mode_d_half_out(i),HPDI_d_half_out(i,:),mode_WHHRI(i),HPDI_WHHRI(i,:),mean_IP(i),mode_IP(i),HPDI_IP(i,:)]=ProcessOutput2(rslts{i},burnin,thin,plotMltpl,doPlots,savePlots);
tmp=[];
for j=1:np
tmp=[tmp,mode_p(i,j),HPDI(j,:,i)];
end
MdlParEsts(i,:)=[tmp,mode_p1(i),HPDI1(i,:)];
end
save(['MdlParEstsFinal' IPD runName])
%% Calculate DIC differences from best-fitting model
[DICs,DICmin,DICdiffs,RMLs,DICminMdl]=RunCalcDICdiffs(rslts,str);
%% Output parameter estimates and DICs to file
MdlParEstsAndDICs=[MdlParEsts,DICs,DICdiffs];
save(['MdlParEstsAndDICs' IPD runName],'MdlParEstsAndDICs')
%% Output acceptance rates
acc_rate=NaN(nMdls,11);
for i=1:nMdls
acc_rate(i,:)=GetAccRates(rslts{i});
end
save(['AccRates' IPD runName],'acc_rate')