-
Notifications
You must be signed in to change notification settings - Fork 0
/
cdcLhoodsTest.m
81 lines (79 loc) · 1.75 KB
/
cdcLhoodsTest.m
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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
function f=cdcLhoodsTest(params,ydataNX)
%Input data - hospitalisations (i.e. without multipliers)
totalInc=0;
xdata=1:36;%[1:15,24:30];%Match MLE
lx=length(xdata);
NNbar=[1464566;3790730;10236474;3984200;2412129];
nbar=length(NNbar);
ydataNX=ydataNX(xdata,:);
mu=[393.0256 999.3054 702.6867 406.0680 177.1958];
sig=[129.6984 329.7708 231.8866 134.0024 58.4746];
ysim=cdcPandemicSimulationW5(params,xdata,0,0,0,243);
%{
%Eliminate age group(s):
ageOut=5;
nbar=nbar-length(ageOut);
ydataNX(:,ageOut)=[];
ysim(:,ageOut)=[];
mu(ageOut)=[];
sig(ageOut)=[];
%}
%ysim=ysim(xdata,:);%simCut in cdcPandemicSimulationW5
if totalInc==1
muvec=sum(mu);
sigvec=sqrt(sum(sig.^2));
hosp=sum(ydataNX,2);
y=sum(ysim,2);
else
muvec=reshape(repmat(mu,lx,1),nbar*lx,1);
sigvec=reshape(repmat(sig,lx,1),nbar*lx,1);
hosp=reshape(ydataNX,nbar*lx,1);
y=reshape(ysim,nbar*lx,1);
end
%}
y=y./hosp;
y(hosp==0)=0;
%
thresh=100;
x=zeros(lx*nbar,1);
%x(y<thresh)=1;
x(hosp<thresh)=1;
y(x==1)=NaN;
if totalInc==0
muvec(x==1)=NaN;
sigvec(x==1)=NaN;
end
%}
%
%Normal likelihood:
L=-log(normpdf(y,muvec,sigvec));
f=nansum(L);
%}
%{
figure
L2=-log(normpdf(zeros(length(y),1),muvec,sigvec));
if totalInc==0
L=reshape(L,lx,nbar);
L2=reshape(L2,lx,nbar);
end
X=[L;L2];
X(isinf(X)==1)=0;
maxy=max(max(X));
fs=12; lw=2;
hold on
plot(xdata,L,'x','markersize',4,'linewidth',lw);
plot(xdata,L2,':','linewidth',lw);
xlabel('Time (weeks)','FontSize',fs);
ylabel('Likelihood','FontSize',fs);
set(gca,'FontSize',fs);
axis([1,xdata(end),0,maxy]);
if totalInc==1
legend('L(sim)','L(0)','location','SE')
end
%legend('0-4','5-17','18-49','50-64','65+','location',NE)
grid on
grid minor
box on
hold off
%}
end