-
Notifications
You must be signed in to change notification settings - Fork 0
/
lhoods_david.asv
41 lines (38 loc) · 1.05 KB
/
lhoods_david.asv
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
function f=lhoods_david(params,ydata)
%Input data - hospitalisations (i.e. without multipliers)
xdata=4:7;%Match MLE
lx=length(xdata);
ydata=ydata(xdata,:);
ydata=reshape(ydata,4*lx,1);
%
%Scale incidence back up to absolute:
NNbar=[19169690;62121035;184015269;39570590];
%NNbar=[19037307;62045041;182377351;38799891];
NNmat=repmat(NNbar',lx,1);
NNvec=kron(NNbar,ones(lx,1));
ydata=ydata.*NNvec;
%Scale ymean also
%}
%Input data without multipliers
na=4;%Number of age groups output - 1 or 4 - ****match with byAge****
qmean=1./[143.4400,364.7100,148.2000,64.6700];
qvec=repmat(qmean,lx,1);
qvec=reshape(qvec,4*lx,1);
%ydata=ydata./qvec;
ysim=pandemic1DallV(params,xdata,0,0,0).*NNmat;
ymean=reshape(ysim,na*lx,1).*qvec;
ysd=ymean.*(1-qvec);
ysd=sqrt(ysd);
%Normal likelihood:
%L=-log(normpdf(ydata,ymean,ysd));
L=-log(mvnpdf(ydata,ymean,ysd));
f=sum(L);
end
%{
function f=Lqtheta(qvec,params,xdata,ydata,lx)
y=pandemic1DallV(params,xdata,0,1,0);
ymean=reshape(y,4*lx,1);
ysd=sqrt(ymean.*(1-qvec));
f=normpdf(ydata,ymean,ysd);
end
%}