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