*! mixlwtp_gf0 1.1.0 05Apr2016 *! author arh version 11.1 mata: void mixlwtp_gf0(M, todo, b, lnfj, S, H) { external mixlwtp_X external mixlwtp_Y external mixlwtp_T external mixlwtp_CSID external mixlwtp_nrep external mixlwtp_np external mixlwtp_kfix external mixlwtp_krnd external mixlwtp_krln external mixlwtp_burn external mixlwtp_corr external mixlwtp_user nrep = mixlwtp_nrep np = mixlwtp_np kfix = mixlwtp_kfix krnd = mixlwtp_krnd krln = mixlwtp_krln burn = mixlwtp_burn corr = mixlwtp_corr user = mixlwtp_user B = b' kall = kfix + krnd if (kfix > 0) { MFIX = B[|1,1\kfix,1|] MFIX = MFIX :* J(kfix,nrep,1) } MRND = B[|(kfix+1),1\kall,1|] if (corr == 1) { external mixlwtp_ncho ncho = mixlwtp_ncho SRND = invvech(B[|(kall+1),1\(kall+ncho),1|]) :* lowertriangle(J(krnd,krnd,1)) } else { SRND = diag(B[|(kall+1),1\(kfix+2*krnd),1|]) } if (user == 1) external mixlwtp_USERDRAWS P = J(np,1,0) i = 1 for (n=1; n<=np; n++) { if (user == 1) { ERR = invnormal(mixlwtp_USERDRAWS[|1,(1+nrep*(n-1))\krnd,(nrep*n)|]) } else { ERR = invnormal(halton(nrep,krnd,(1+burn+nrep*(n-1)))') } if (kfix > 0) BETA = MFIX \ (MRND :+ (SRND*ERR)) else BETA = MRND :+ (SRND*ERR) if (krln > 0) { if ((kall-krln) > 0) { BETA = BETA[|1,1\(kall-krln),nrep|] \ exp(BETA[|(kall-krln+1),1\kall,nrep|]) } else { BETA = exp(BETA) } } R = J(1,nrep,1) t = 1 nc = mixlwtp_T[i,1] for (t=1; t<=nc; t++) { YMAT = mixlwtp_Y[|i,1\(i+mixlwtp_CSID[i,1]-1),cols(mixlwtp_Y)|] XMAT = mixlwtp_X[|i,1\(i+mixlwtp_CSID[i,1]-1),cols(mixlwtp_X)|] B1 = BETA[1..(rows(BETA)-1),.] \ J(1,cols(BETA),1) B2 = J(rows(XMAT),cols(BETA),1) :* BETA[rows(BETA),.] V = XMAT * B1 V = B2 :* V EV = exp(V) :/ colsum(exp(V),1) R = R :* colsum(YMAT :* EV,1) i = i + mixlwtp_CSID[i,1] } P[n,1] = mean(R',1) } lnfj = ln(P) } end exit