1C -------------------------------------------------------------------
2C MENGWONG2 (no missing values) computes the Marginal Lilkelihood estimates as deteiled
3C by Meng and Wong, Statistica Sinica, 1996
4C Developed by A.Rossi, C.Planas and G.Fiorentini
5C
6C OUTPUT:
7C  MLMW(:,1) all parameters,
8C  MLMW(:,2) non-var params
9C  MLMW(1,:) no iteration,
10C  MLMW(2,:) SD,
11C  MLMW(3,:) 10 iterations
12C
13C Remarks:
14C  NPAR is total # of params,
15C  NPARD = NPAR - #Variances
16C
17C Copyright (C) 2010-2014 European Commission
18C
19C This file is part of Program DMM
20C
21C DMM is free software developed at the Joint Research Centre of the
22C European Commission: you can redistribute it and/or modify it under
23C the terms of the GNU General Public License as published by
24C the Free Software Foundation, either version 3 of the License, or
25C (at your option) any later version.
26C
27C DMM is distributed in the hope that it will be useful,
28C but WITHOUT ANY WARRANTY; without even the implied warranty of
29C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
30C GNU General Public License for more details.
31C
32C You should have received a copy of the GNU General Public License
33C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
34C -------------------------------------------------------------------
35	SUBROUTINE MENGWONG2(G,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
36	1                     INFOS,yk,gibpar,gibZ,thetaprior,psiprior,
37     2                     tipo,pdll,MLSTART,MLMW)
38
39C INPUT
40      INTEGER G,nobs,d(2),ny,nz,nx,nu,nv,ns(6),nstot,nt,np(3),
41	1 INFOS(9,6),gibZ(G,nobs)
42	DOUBLE PRECISION yk(nobs,ny+nz),gibpar(G,nt+np(1)),
43	1 thetaprior(nt,4),psiprior(np(2),np(3)),MLSTART
44	CHARACTER*2 tipo(nt)
45	POINTER (pdll,fittizia)
46
47C OUTPUT
48	DOUBLE PRECISION MLMW(2,2)
49
50C LOCALS
51	INTEGER NPAR,I,J,K,IG,NPOS(nt+np(1)),IFAIL,NQ,ISEQ,ISEQ0,SEQ(nv),
52	1 IS(nobs,6),NIM,NI,IND(1),NPARTH,NN,NSI,II,JJ
53	DOUBLE PRECISION,ALLOCATABLE::MAT(:,:),VQN(:,:),VQD(:,:),
54	1 VHN(:,:),VHD(:,:)
55	DOUBLE PRECISION parm(nt),SIGM(nt,nt),
56	1 COM(nt+1,nt),ISIGM(nt,nt),par(nt+np(1)),SEGA(nt+np(1)),
57     2 ub(nt),lb(nt),R3((nt+1)*(nt+2)/2),WORK(nt)
58	DOUBLE PRECISION,ALLOCATABLE:: PTR(:,:,:),PMAT(:,:),PE(:),GAM(:),
59	1 ALPHA(:,:),MOM(:,:)
60	DOUBLE PRECISION P1(INFOS(8,1),INFOS(8,1)),
61	1 P2(INFOS(8,2),INFOS(8,2)),P3(INFOS(8,3),INFOS(8,3)),
62     2 P4(INFOS(8,4),INFOS(8,4)),P5(INFOS(8,5),INFOS(8,5)),
63     3 P6(INFOS(8,6),INFOS(8,6))
64	DOUBLE PRECISION Ppar(nt+np(1)),Fpar,PS,QS,QPSI,C,DET,TRC,A0
65	DOUBLE PRECISION ERRM,ERR,U,AUX,INDC(1),MUC,SS(2,2),MWNUM,MWDEN
66	DOUBLE PRECISION ZERO,ONE,PI
67	DATA ZERO/0.0D0/,ONE/1.0D0/,PI/3.141592653589793D0/
68
69C EXTERNAL SUBROUTINES
70      EXTERNAL NEWEYWESTCOV,NEWEYWESTCOV2,mvncdf,DPOTRF,DPOTRI,setgmn,
71     1 genmn,gengam,DESIGNZ,PPROD,ERGODIC,INT2SEQ
72
73C EXTERNAL FUNCTIONS
74      DOUBLE PRECISION PTHETA2,PRIOR,PRIORDIR,genunf,gengam
75
76      PAR(:) = GIBPAR(1,:) ! set constant values
77	NPARTH = 0
78	DO I = 1,nt
79	 IF (GIBPAR(1,I).NE.GIBPAR(2,I)) THEN
80	  NPARTH = NPARTH + 1
81        NPOS(NPARTH) = I
82	 ENDIF
83      ENDDO
84	DO I = 1,np(1)
85       NPOS(NPARTH+I) = nt+I
86      ENDDO
87	NPAR=NPARTH+np(1)
88      Ppar(:) = 0.D0
89	parm(:) = ZERO
90	DO I = 1,NPARTH
91	 parm(I) = SUM(gibpar(:,NPOS(I)))/DFLOAT(G)
92      ENDDO
93
94	NQ = 0
95	CALL NEWEYWESTCOV2(G,NPARTH,NQ,gibpar(:,NPOS(1:NPARTH)),
96	1  parm(1:NPARTH),SIGM(1:NPARTH,1:NPARTH)) ! THETA Var-covar
97
98	IF (nv.GT.0) THEN
99	 ALLOCATE(PTR(nobs,nstot,nstot),PMAT(nstot,nstot),PE(nstot))
100
101C Transition prob for QS
102	 DO I = 1,nstot-1
103	  PTR(1,I,1) = SUM(ABS(gibZ(1:G,1).EQ.I))/DFLOAT(G)
104       ENDDO
105       PTR(1,nstot,1) = ONE-SUM(PTR(1,1:nstot-1,1))
106
107	 DO 50 K = 2,nobs
108	  DO 50 I = 1,nstot-1
109	   DO 50 J = 1,nstot
110	    COM(1,1) = SUM(ABS(gibZ(1:G,K-1).EQ.J))
111	    IF (COM(1,1).GT.ZERO) THEN
112	     PTR(K,I,J) = SUM(ABS((gibZ(1:G,K).EQ.I).AND.
113     #		          (gibZ(1:G,K-1).EQ.J)))/COM(1,1)
114	    ELSE
115	     PTR(K,I,J) = ONE/DFLOAT(nstot)
116	    ENDIF
11750	   PTR(K,nstot,J) = ONE-SUM(PTR(K,1:nstot-1,J))
118
119C Mean and Var of PSI
120	 ALLOCATE (ALPHA(np(2),np(3)),MOM(np(1),2))
121	 DO I=1,np(1)
122	  MOM(I,1) = SUM(gibpar(:,nt+I))/DFLOAT(G)
123	  MOM(I,2) = SUM(gibpar(:,nt+I)**2)/DFLOAT(G)
124	  MOM(I,2) = MOM(I,2)-MOM(I,1)**2
125	 ENDDO
126C Hyperparameters of Dirichelt for Q(PSI)
127C Mothod of Moments: a0 = m1(1-m1)/V1+1, ai = mi*a0, i=1,2,..,N
128	 NN = 0
129	 K  = 0
130	 DO I = 1,nv
131	  NSI = INFOS(8,I)           ! # of states for S
132	  IF (INFOS(9,I).EQ.1) THEN  ! S~IID
133	   A0 = MOM(NN+1,1)*(1.D0-MOM(NN+1,1))/MOM(NN+1,2)+1.D0 !alpha0
134	   DO  ii = 1,NSI-1
135          ALPHA(K+1,ii) = MOM(NN+ii,1)*A0
136         ENDDO
137	   ALPHA(K+1,NSI) = A0-SUM(ALPHA(K+1,1:NSI-1))
138	   K  = K + 1
139         NN = NN + NSI-1
140	  ELSEIF (INFOS(9,I).EQ.2) THEN  ! S~Markov
141	   DO jj = 1,NSI
142	    A0 = MOM(NN+1,1)*(1.D0-MOM(NN+1,1))/MOM(NN+1,2)+1.D0 !alpha0
143	    DO  ii = 1,NSI-1
144           ALPHA(K+1,ii) = MOM(NN+ii,1)*A0
145          ENDDO
146	    ALPHA(K+1,NSI) = A0-SUM(ALPHA(K+1,1:NSI-1))
147          K  = K + 1
148	    NN = NN + NSI-1
149         ENDDO
150	  ENDIF
151	 ENDDO
152	ENDIF
153C Importance sampling
154C Sample THIS from N(THHAT,SIGHAT) with boundaries
155C Evaluate p(THIS) ~ N(THHAT,SIGHAT)
156	NIM  = 1000000
157	ERRM = 1.D-8
158
159C Normalization constants TRC and TRCD
160      lb(1:NPARTH) = thetaprior(NPOS(1:NPARTH),3)
161      ub(1:NPARTH) = thetaprior(NPOS(1:NPARTH),4)
162	CALL mvncdf(lb(1:NPARTH),ub(1:NPARTH),parm(1:NPARTH),
163	1   SIGM(1:NPARTH,1:NPARTH),NPARTH,ERRM,NIM,TRC,ERR,NI)
164
165C Inverse SIGM  & det for NPARTH
166	COM(1:NPARTH,1:NPARTH) = SIGM(1:NPARTH,1:NPARTH)
167	IFAIL = -1
168C	CALL F01ADF(NPARTH,COM(1:NPARTH+1,1:NPARTH),NPARTH+1,IFAIL) ! Inverse var-covar
169      CALL DPOTRF('L',NPARTH,COM(1:NPARTH,1:NPARTH),NPARTH,IFAIL) ! COM = L*L'
170      DET = 1.D0 ! det(SIGM)
171      DO I=1,NPARTH
172	 DET = DET*COM(I,I)**2
173      ENDDO
174      CALL DPOTRI('L',NPARTH,COM(1:NPARTH,1:NPARTH),NPARTH,IFAIL) ! COM = VV^-1
175
176      DO 60 I=1,NPARTH
177	ISIGM(I,I) = COM(I,I)
178	DO 60 J=1,I-1
179	ISIGM(I,J) = COM(I,J)
18060	ISIGM(J,I) = ISIGM(I,J)
181
182c	COM(1:NPARTH,1:NPARTH) = SIGM(1:NPARTH,1:NPARTH)
183c      IFAIL = -1
184c      CALL F03ABF(COM(1:NPARTH,1:NPARTH),NPARTH,NPARTH,DET,
185c	1 WORK(1:NPARTH),IFAIL)
186
187      C = (2.D0*PI)**(-.5D0*NPARTH)/DSQRT(DET) ! constant
188
189	ALLOCATE (MAT(G,2),VHN(G,2),VHD(G,2),VQN(G,2),VQD(G,2))
190	QS = ONE
191	PS = ONE
192	IS(:,:) = 1
193	DO 200 IG = 1,G
194C SAMPLING THETA
195       SEGA(:)  = -1.D0
196 	 IFAIL    = -1
197	 IND(1)   =  0
198	 INDC(1)  = -1.D0
199	 DO WHILE (INDC(1).LT.ZERO)
200	  INDC(1) = ZERO
201	  IND(1)  = IND(1) + 1
202	  IF (IND(1).GT.G) EXIT
203c	  CALL G05EAF(parm(1:NPARTH),NPARTH,SIGM(1:NPARTH,1:NPARTH),
204c	1              NPARTH,EPS,R3,(NPARTH+1)*(NPARTH+2)/2,IFAIL)
205c	  CALL G05EZF(SEGA(1:NPARTH),NPARTH,R3,(NPARTH+1)*(NPARTH+2)/2,
206c	1              IFAIL)
207        COM(1:NPARTH,1:NPARTH) = SIGM(1:NPARTH,1:NPARTH)
208        CALL setgmn(parm(1:NPARTH),COM(1:NPARTH,1:NPARTH),NPARTH,
209     1              NPARTH,R3(1:(NPARTH+2)*(NPARTH+1)/2))
210        CALL genmn(R3(1:(NPARTH+2)*(NPARTH+1)/2),SEGA(1:NPARTH),
211     1             WORK(1:NPARTH))
212	  DO  I=1,NPARTH
213	   IF (SEGA(I).LT.thetaprior(NPOS(I),3)) INDC(1)=-1
214         IF (SEGA(I).GT.thetaprior(NPOS(I),4)) INDC(1)=-2
215        ENDDO
216       END DO
217C SAMPLING PSI from Dirichlet(ALPHA)
218	 NN = NPARTH
219	 K  = 0
220	 DO 70 I = 1,nv
221	  NSI = INFOS(8,I) ! # of states for SI
222	  ALLOCATE(GAM(NSI))
223	  IF (INFOS(9,I).EQ.1) THEN ! S~IID
224	   DO  ii = 1,NSI
225	    IFAIL = -1
226C	    CALL G05FFF(ALPHA(K+1,ii),1.D0,1,GAM(ii),IFAIL)
227          GAM(ii) = gengam(1.D0,ALPHA(K+1,ii))
228         ENDDO
229	   SEGA(NN+1:NN+NSI-1) = GAM(1:NSI-1)/SUM(GAM(1:NSI))
230	   K  = K + 1
231         NN = NN + NSI-1
232	  ELSEIF (INFOS(9,I).EQ.2) THEN  ! S~Markov
233	   DO jj = 1,NSI
234	    DO ii = 1,NSI
235		 IFAIL = -1
236C	     CALL G05FFF(ALPHA(K+1,ii),1.D0,1,GAM(ii),IFAIL)
237           GAM(ii) = gengam(1.D0,ALPHA(K+1,ii))
238          ENDDO
239	    SEGA(NN+1:NN+NSI-1) = GAM(1:NSI-1)/SUM(GAM(1:NSI))
240          K  = K + 1
241	    NN = NN + NSI-1
242	   ENDDO
243        ENDIF
24470     DEALLOCATE(GAM)
245
246C SAMPLING S
247	 IF (nv.GT.0) THEN
248	  CALL DESIGNZ(nv,np(1),SEGA(NPARTH+1:NPAR),INFOS,
249	1               P1,P2,P3,P4,P5,P6)
250C PMAT(i,j) = Pr[Z(t+1)=i|Z(t)=j], Z = S1 x S2 x ... x Snv
251	  CALL PPROD(nv,nstot,INFOS,P1,P2,P3,P4,P5,P6,PMAT)
252C ERGODIC solves PE: PE*(I-P') = 0
253        CALL ERGODIC(nstot,PMAT,PE)
254C S(1)
255C	  U = G05CAF(U) ! Sampling from U(0,1)
256        U = genunf(0.D0,1.D0)
257	  ISEQ = 1
258	  AUX  = PTR(1,ISEQ,1)
259	  DO 80 WHILE (AUX.LT.U)
260	   ISEQ = ISEQ + 1
26180	  AUX  = AUX  + PTR(1,ISEQ,1)
262	  CALL INT2SEQ(ISEQ,nv,INFOS,SEQ,IS(1,:))
263	  QS = PTR(1,ISEQ,1)
264	  PS = PE(ISEQ) ! P(S1)
265	  ISEQ0 = ISEQ
266C S(2),...,S(nobs)
267	  DO 90 K = 2,nobs
268C	   U = G05CAF(U) ! Sampling from U(0,1)
269         U = genunf(0.D0,1.D0)
270	   ISEQ = 1
271	   AUX  = PTR(K,ISEQ,ISEQ0)
272	   DO 85 WHILE (AUX.LT.U)
273	   ISEQ = ISEQ + 1
27485	   AUX  = AUX  + PTR(K,ISEQ,ISEQ0)
275	   CALL INT2SEQ(ISEQ,nv,INFOS,SEQ,IS(K,:))
276	   QS = QS*PTR(K,ISEQ,ISEQ0)
277	   PS = PS*PMAT(ISEQ,ISEQ0)
27890	  ISEQ0 = ISEQ
279	 ENDIF
280
281C QUADRATIC FORM FOR for THETA
282       DO 91 I = 1,NPARTH
28391     COM(I,1) = SUM((SEGA(1:NPARTH)-parm(1:NPARTH))*ISIGM(1:NPARTH,I))
284	 MUC = SUM(COM(1:NPARTH,1)*(SEGA(1:NPARTH)-parm(1:NPARTH)))
285
286C	 VQN(IG,1) = QS*C*DEXP(-.5D0*MUC)/TRC
287
288	 par(NPOS(1:NPARTH+np(1))) = SEGA(1:NPARTH+np(1)) ! (THETA,PSI)
289
290C PRIOR for THETA
291	 DO 92 I = 2,NPARTH
29292	 Ppar(I) = PRIOR(par(NPOS(I)),thetaprior(NPOS(I),:),tipo(NPOS(I)))
293
294C PRIOR for PSI and Q(PSI)~Dirichlet(a1,a2,...,aN)
295	 QPSI = 0.D0
296	 NN = NPARTH
297	 K  = 0
298	 DO 100 J = 1,nv
299	  NSI = INFOS(8,J)
300	  IF(INFOS(9,J).EQ.1) THEN ! S~IID
301	   Ppar(NPARTH+K+1) = PRIORDIR(par(NPOS(NN+1:NN+NSI-1)),
302	1                               psiprior(K+1,1:NSI),NSI)
303	   QPSI = QPSI+PRIORDIR(par(NPOS(NN+1:NN+NSI-1)),
304	1                        ALPHA(K+1,1:NSI),NSI)
305         K  = K + 1
306         NN = NN + NSI-1
307	  ELSEIF(INFOS(9,J).EQ.2) THEN ! S~Markov
308	   DO 99 I = 1,NSI
309	    Ppar(NPARTH+K+1) = PRIORDIR(par(NPOS(NN+1:NN+NSI-1)),
310	1                                psiprior(K+1,1:NSI),NSI)
311  	    QPSI = QPSI+PRIORDIR(par(NPOS(NN+1:NN+NSI-1)),
312	1                         ALPHA(K+1,1:NSI),NSI)
313		K  = K + 1
31499	    NN = NN + NSI-1
315	  ENDIF
316100    CONTINUE
317
318	 Fpar = PTHETA2(NPOS(1),nobs,d,ny,nz,nx,nu,ns,nt,IS,yk,
319	1                par(1:nt),thetaprior(NPOS(1),:),
320     2                tipo(NPOS(1)),pdll)
321       Fpar = Fpar + SUM(Ppar(2:NPARTH+K)) ! log f(y|par,S)f(par,S)
322
323	 VQN(IG,1) = DEXP(QPSI)*QS*C*DEXP(-.5D0*MUC)/TRC
324
325200	 VHN(IG,1) = Fpar + DLOG(PS)
326
327C ---------------------
328C Meng-Wong denominator
329C ---------------------
330	QS = ONE
331	PS = ONE
332	DO 400 IG = 1,G
333   	 IF (nv.GT.0) THEN
334	  CALL DESIGNZ(nv,np(1),gibpar(IG,nt+1:nt+np(1)),INFOS,
335	1               P1,P2,P3,P4,P5,P6)
336C PMAT(i,j) = Pr[Z(t+1)=i|Z(t)=j], Z = S1 x S2 x ... x Snv
337	  CALL PPROD(nv,nstot,INFOS,P1,P2,P3,P4,P5,P6,PMAT)
338C ERGODIC solves PE: PE*(I-P') = 0
339        CALL ERGODIC(nstot,PMAT,PE)
340
341	  QS = PTR(1,gibZ(IG,1),1)
342	  PS = PE(gibZ(IG,1))
343	  CALL INT2SEQ(gibZ(IG,1),nv,INFOS,SEQ,IS(1,:))
344	  DO 210 K = 2,nobs
345	   QS = QS*PTR(K,gibZ(IG,K),gibZ(IG,K-1))
346	   PS = PS*PMAT(gibZ(IG,K),gibZ(IG,K-1))
347210	   CALL INT2SEQ(gibZ(IG,K),nv,INFOS,SEQ,IS(K,:))
348       ENDIF
349
350C PRIOR for THETA
351       DO 310 I = 2,NPARTH
352310	 Ppar(I) = PRIOR(gibpar(IG,NPOS(I)),thetaprior(NPOS(I),:),
353     1                 tipo(NPOS(I)))
354
355C PRIOR for PSI and Q(PSI)~Dirichlet(a1,a2,...,aN)
356	 QPSI = 0.D0
357	 NN = NPARTH
358	 K  = 0
359	 DO 305 J = 1,nv
360	  NSI = INFOS(8,J)
361	  IF(INFOS(9,J).EQ.1) THEN ! S~IID
362	   Ppar(NPARTH+K+1) = PRIORDIR(gibpar(IG,NPOS(NN+1:NN+NSI-1)),
363	1                               psiprior(K+1,1:NSI),NSI)
364	   QPSI = QPSI+PRIORDIR(gibpar(IG,NPOS(NN+1:NN+NSI-1)),
365	1                        ALPHA(K+1,1:NSI),NSI)
366         K  = K + 1
367         NN = NN + NSI-1
368	  ELSEIF(INFOS(9,J).EQ.2) THEN ! S~Markov
369	   DO 304 I = 1,NSI
370	    Ppar(NPARTH+K+1) = PRIORDIR(gibpar(IG,NPOS(NN+1:NN+NSI-1)),
371	1                                psiprior(K+1,1:NSI),NSI)
372	    QPSI = QPSI+PRIORDIR(gibpar(IG,NPOS(NN+1:NN+NSI-1)),
373	1                         ALPHA(K+1,1:NSI),NSI)
374		K  = K + 1
375304	    NN = NN + NSI-1
376	  ENDIF
377305    CONTINUE
378
379	 Fpar = PTHETA2(NPOS(1),nobs,d,ny,nz,nx,nu,ns,nt,IS,yk,
380	1                gibpar(IG,1:nt),thetaprior(NPOS(1),:),
381     2                tipo(NPOS(1)),pdll)
382       Fpar = Fpar + SUM(Ppar(2:NPARTH+K)) ! log f(y|par,S)f(par,S)
383
384	 VHD(IG,1) = Fpar + DLOG(PS)
385
386	 COM(:,1) = ZERO
387	 DO 320 I = 1,NPARTH
388320	 COM(I,1) = SUM((gibpar(IG,NPOS(1:NPARTH))-parm(1:NPARTH))
389     #          * ISIGM(1:NPARTH,I))
390	 MUC = SUM(COM(1:NPARTH,1)*(gibpar(IG,NPOS(1:NPARTH))
391     #	 - parm(1:NPARTH)))
392
393	 VQD(IG,1) = DEXP(QPSI)*QS*DEXP(-.5D0*MUC)*C/TRC
394400   CONTINUE
395
396	IND = MAXLOC(VHN(:,1))
397	DET = VHN(IND(1),1)
398
399      MAT(:,1) = DEXP(VHN(:,1)-DET)/(DEXP(VHN(:,1)-MLSTART)+VQN(:,1))
400      MAT(:,2) = VQD(:,1)/(DEXP(VHD(:,1)-MLSTART)+VQD(:,1))
401
402      CALL NEWEYWESTCOV(G,2,1,MAT(:,1:2),SS)
403      MLMW(2,1)   = SUM(MAT(:,1))/SUM(MAT(:,2))
404      MLMW(1,1)   = DLOG(MLMW(2,1)) + DET
405      MLMW(1:2,2) = SS(1,1)*G/SUM(MAT(:,1))**2 +
406     +              SS(2,2)*G/SUM(MAT(:,2))**2 +
407     +     	    - 2.D0*SS(1,2)*G/(SUM(MAT(:,1))*SUM(MAT(:,2)))
408
409      MLMW(2,1) = MLMW(1,1) ! log scale
410      DO 500 I=1,10
411       MWNUM     = SUM(DEXP(VHN(:,1)-DET)
412	1           / (DEXP(VHN(:,1)-MLMW(2,1))+VQN(:,1)))
413	 MWDEN     = SUM(VQD(:,1)/(DEXP(VHD(:,1)-MLMW(2,1))+VQD(:,1)))
414	 MLMW(2,1) = DLOG(MWNUM/MWDEN) + DET	! log-scale
415500   CONTINUE
416
417	DEALLOCATE (MAT,VHN,VHD,VQN,VQD)
418	IF (nv.GT.0) DEALLOCATE (PTR,PMAT,PE,ALPHA,MOM)
419
420	RETURN
421	END
422