https://github.com/cran/HLSM
Tip revision: dcd2fa6f1994b6500075e0b6f0310a45f112b45f authored by Samrachana Adhikari on 18 June 2014, 00:00:00 UTC
version 0.2
version 0.2
Tip revision: dcd2fa6
HLSM.cpp
#include <R.h>
#include <R_ext/Print.h>
#include <iostream>
#include <math.h>
#include <Rmath.h>
//small functions used later in MCMC
double square(double x){
double y = x*x;
return y ;
}
double logitInverse(double x){
return 1.0/(1.0 + exp(-1.0 * x));
}
//read row ii of a matrix with dd columns
void ReadRow(double *XX, int ii, double *Row, int *nn, int *dd)
{
for(int kk =0; kk < dd[0]; kk++)
{
Row[kk] = XX[ii + kk*nn[0]];
}
}
//write row ii in a matrix
void WriteRow(double *XX, int ii, double *Xvec, int *nn, int *dd)
{
for(int kk = 0; kk < dd[0]; kk ++)
{
XX[ii + kk*nn[0]] = Xvec[kk];
}
}
//extract Y and Z matrix for KK network with nn rows and pp columns
void getY(double *Mat,double *NewMat, int *nn,int kk){
int ss = 0.0;
if(kk > 0){
for(int ll = 0;ll < kk; ll++){
ss += nn[ll]*nn[ll];
}
}
for(int xx= 0;xx < nn[kk]; xx++){
for(int yy =0; yy < nn[kk]; yy++){
NewMat[xx+yy*nn[kk]] = Mat[xx+yy*nn[kk]+ss];
}
}
}
void getZ(double *Mat,double *NewMat, int *nn,int dd,int kk){
int ss = 0.0;
if(kk > 0){
for(int ll = 0;ll < kk; ll++){
ss += (dd*nn[ll]);
}
}
for(int xx= 0;xx < nn[kk]; xx++){
for(int yy =0; yy < dd; yy++){
NewMat[xx+yy*nn[kk]] = Mat[xx+yy*nn[kk]+ss];
}
}
}
//ss = sum(nn[1:kk] for kk groups
void readX(double *X,double *newX, int *nn, int pp, int kk){
int ss = 0.0;
if(kk > 0){
for(int ll = 0;ll < kk; ll++){
ss += nn[ll]*nn[ll];
}
}
for(int xx = 0; xx < nn[kk]; xx++){
for(int yy = 0; yy < nn[kk];yy ++){
for(int zz =0; zz <pp; zz++){
newX[xx+yy*nn[kk]+zz*nn[kk]*nn[kk]] = X[xx+yy*nn[kk]+zz*nn[kk]*nn[kk]+pp*ss];
}
}
}
}
//compute distance matrix
void distMat(int *nn, int *dd, double *ZZ, double *dMat){
int ii,jj,kk;
double tmp;
for(ii = 0 ; ii <= (nn[0]-1) ; ii++){
for(jj = 0 ; jj <= ii ; jj++){
tmp = 0.0;
for(kk = 0 ; kk < dd[0] ; kk++){
double VV = ZZ[ii+kk*nn[0]] - ZZ[jj+kk*nn[0]];
tmp = tmp + square(VV);
}
dMat[jj*nn[0]+ii] = sqrt(tmp);
dMat[ii*nn[0]+jj] = sqrt(tmp);
}
}
}
//Compute logpriors
void LogpriorBeta(double *Beta, double *MuBeta, double *SigmaBeta, double *val)
{
double sigma = sqrt(SigmaBeta[0]);
val[0] = dnorm(Beta[0], MuBeta[0], sigma, 1);
}
void LogpriorAlpha(double *Alpha, double *MuAlpha, double *SigmaAlpha, double *val)
{
double sigma = sqrt(SigmaAlpha[0]);
val[0] = dnorm(Alpha[0],MuAlpha[0],sigma,1);
}
//Z is treated as independent normals in each coordinate
void LogpriorZZ(double *ZZ, double *Mu, double *Var, int *dd, double *val){
int ii;
double total = 0.0;
for(ii = 0 ; ii < dd[0] ; ii++){
double sigma = sqrt(Var[ii]);
total = total + dnorm(ZZ[ii],Mu[ii],sigma,1);
}
val[0] = total;
}
//compute loglikelihood for a single network
void FullLogLik(double *beta, double *YY, double *XX, double *ZZ, double *alpha, int *Tr, double *intercept,int *nn, int *pp, int *dd, double *Val){
double* dMat = 0;
dMat = new double[nn[0]*nn[0]];
distMat(nn,dd,ZZ,dMat);
double total = 0.0;
double tmp1,tmp2, tmpVal1, tmpVal2;
for(int ii = 1; ii < nn[0]; ii++){
for(int jj = 0; jj < ii; jj++){
tmpVal1 = 0.0;
tmpVal2 = 0.0;
//compute t(X)*beta for all P
for(int kk = 0; kk < pp[0]; kk++){
tmpVal1 = tmpVal1 + beta[kk]*XX[ii+jj*nn[0]+kk*nn[0]*nn[0]];
tmpVal2 = tmpVal2 + beta[kk]*XX[jj + ii*nn[0] + kk*nn[0]*nn[0]];
}
double vv1 = (intercept[0] + tmpVal1 + alpha[0]*Tr[0] - dMat[jj*nn[0] + ii]);
double vv2 = (intercept[0] + tmpVal2 + alpha[0]*Tr[0] - dMat[jj+ii*nn[0]]);
tmp1 = logitInverse(vv1);
tmp2 = logitInverse(vv2);
if(YY[jj*nn[0] + ii] == 1){
total = total + log(tmp1);
}else if(YY[jj*nn[0] + ii] == 0){
total = total + log(1.0 - tmp1);
}
if(YY[jj + ii*nn[0]] == 1){
total = total + log(tmp2);
}else if(YY[jj + ii*nn[0]] == 0){
total = total + log(1.0 - tmp2);
}
}
}
Val[0] = total;
delete[] dMat;
}
//sum of loglikelihoods for all KK networks
void AllLogLik(double *X, double *Y, double *Z, int *T, int *nn, int *pp, int *dd, int *KK, double *beta, double *intercept,double *alpha, double *lliknew)
{
double Val;
// int slen = KK[0]+1;
// double* sumn = 0; // need this to store the updated parameters
// sumn = new double[slen];
// sumn[0] = 0;
// for(int iz =0; iz < KK[0]; iz++){
// sumn[iz+1] = sumn[iz] + nn[iz];
// }
lliknew[0] = 0.0;
for(int kk = 0; kk < KK[0]; kk++){
// int ss = sumn[kk];
double* XMat = 0;
XMat = new double[nn[kk]*nn[kk]*pp[0]];
double* YMat = 0;
YMat = new double[nn[kk]*nn[kk]];
double* ZMat = 0;
ZMat = new double[nn[kk]*dd[0]];
getY(Y,YMat,nn,kk);
readX(X,XMat,nn,pp[0],kk);
getZ(Z,ZMat,nn,dd[0],kk);
//use the X Y and Z for that group.
FullLogLik(beta, YMat, XMat, ZMat, alpha, &T[kk], intercept,&nn[kk],pp,dd,&Val);
lliknew[0] = lliknew[0] + Val;
delete[] XMat;
delete[] YMat;
delete[] ZMat;
}
// delete[] sumn;
}
//update Z as an independent normal
void updateZ(double *XX,double *YY,double *ZZ,int *TT,int *nn, int *pp,int *dd,double *beta,double *intercept,double *alpha,double *zzPrior,double *Var,double *tuneZ,double *llOld,double *accZ)
{
double draw, logRR;
double* Znew = 0;
Znew = new double[nn[0]*dd[0]];
double lpNew;
double llNew;
double lpOld;
for(int ss = 0 ; ss < nn[0]*dd[0] ; ss++){
Znew[ss] = ZZ[ss];
}
for(int ii = 0 ; ii < nn[0] ; ii++){ //ii denotes row
//note that if we were to use the constraints commented below
// we will keep the first entry of Z fixed. ie start from ii = 1
double* ZZsm = 0;
ZZsm = new double[dd[0]];
// double ZZsm[dd[0]];
ReadRow(ZZ,(ii),ZZsm,nn,dd);
LogpriorZZ(ZZsm,zzPrior,Var,dd,&lpOld);
double* ZnewSm = 0;
ZnewSm = new double[dd[0]];
// double ZnewSm[dd[0]];
for(int kk = 0 ; kk < dd[0] ; kk++){ //kk denotes number of columns
ZnewSm[kk] = ZZsm[kk] + tuneZ[ii]*rnorm(0.0,1.0);
// if(ii == 1 && kk == 0){ //Positivity constraint on Znew[2,1]
// if(ZnewSm[kk] < Znew[0]){
// ZnewSm[kk] = -1.0*(ZnewSm[kk]-Znew[0])+1.0;
// }
// }
// if(ii == 1 && kk == 1){ //set Znew[2,2] = 0
// ZnewSm[kk] = 0.0;
// }
// if((kk == 1 && ii == 2)){ //Positivity constraint on Znew[3,2]
// if(ZnewSm[kk] < 0){
// ZnewSm[kk] = -1.0 * ZnewSm[kk];
// }
// }
}
WriteRow(Znew, (ii), ZnewSm,nn,dd);
FullLogLik(beta, YY, XX, Znew, alpha, TT, intercept,nn,pp,dd,&llNew);
LogpriorZZ(ZnewSm,zzPrior,Var,dd,&lpNew);
logRR = llNew - llOld[0] + lpNew - lpOld;
draw = runif(0.0,1.0);
if(log(draw) < logRR){
WriteRow(ZZ, (ii), ZnewSm,nn,dd);
llOld[0] = llNew;
accZ[ii] = accZ[ii] + 1;
}else{
WriteRow(Znew,ii,ZZsm,nn,dd);
}
delete[] ZnewSm;
delete[] ZZsm;
}
delete[] Znew;
}
//UPDATES FOR RANDOM EFFECT MODEL
//update beta
void updateBetamulti(double *X,double *Y,double *Z,int *T,int *nn,int *pp,int *dd, double *beta, double *intercept,double *alpha, double *mu,double *sigmasq,double *tune,double *llik,double *acc)
{
double* betanew = 0;
betanew = new double[pp[0]];
// double betanew[pp[0]];
for(int k =0;k<pp[0];k++){
betanew[k] = beta[k];
}
double lpBetanew;
double llikBetanew;
double lpBeta;
for(int ii = 0; ii < pp[0]; ii++){
LogpriorBeta(&beta[ii], &mu[ii], &sigmasq[ii],&lpBeta);
betanew[ii] = beta[ii] + tune[ii]*rnorm(0.0,1.0);
LogpriorBeta(&betanew[ii],&mu[ii],&sigmasq[ii],&lpBetanew);
FullLogLik(betanew, Y, X, Z, alpha, T, intercept,nn, pp, dd, &llikBetanew);
double logratio = lpBetanew - lpBeta + llikBetanew - llik[0];
if(log(runif(0.0,1.0)) < logratio){
beta[ii] = betanew[ii];
llik[0] = llikBetanew;
// lpBeta[ii] = lpBetanew;
acc[ii] = acc[ii] + 1;
}else{ //if not accepted reset betanew[ii] with old beta
betanew[ii] = beta[ii];
}
}
delete[] betanew;
}
//update intercept
void updateIntercept(double *X,double *Y,double *Z,int *T,int *nn,int *pp,int *dd,double *beta,double *intercept,double *alpha,double *mu,double *sigmasq,double *tuneInt,double *lpInt,double *llik,double *acc){
double intnew;
double lpIntNew;
double llikIntnew;
intnew = intercept[0] + tuneInt[0]*rnorm(0.0,1.0);
LogpriorBeta(&intnew, mu,sigmasq,&lpIntNew);
FullLogLik(beta,Y,X,Z,alpha,T,&intnew,nn,pp,dd,&llikIntnew);
double logratio = lpIntNew-lpInt[0]+llikIntnew-llik[0];
if(log(runif(0.0,1.0)) < logratio){
intercept[0] = intnew;
llik[0] = llikIntnew;
lpInt[0] = lpIntNew;
acc[0] = acc[0] + 1;
}
}
//UPDATES FOR FIXED EFFECT MODEL
//update beta
void updateBetamultiFixedEF(double *X,double *Y,double *Z,int *T,int *nn,int *pp,int *dd, int *KK, double *beta, double *intercept,double *alpha, double *mu,double *sigmasq,double *tune,double *llik,double *acc)
{
double* betanew = 0;
betanew = new double[pp[0]];
// double betanew[pp[0]];
for(int k =0;k<pp[0];k++){
betanew[k] = beta[k];
}
double lpBetanew;
double llikBetanew;
double lpBeta;
for(int ii = 0; ii < pp[0]; ii++){
LogpriorBeta(&beta[ii], &mu[ii], &sigmasq[ii],&lpBeta);
betanew[ii] = beta[ii] + tune[ii]*rnorm(0.0,1.0);
LogpriorBeta(&betanew[ii],&mu[ii],&sigmasq[ii],&lpBetanew);
AllLogLik(X, Y, Z, T, nn, pp, dd,KK, betanew, intercept,alpha, &llikBetanew);
double logratio = lpBetanew - lpBeta + llikBetanew - llik[0];
if(log(runif(0.0,1.0)) < logratio){
beta[ii] = betanew[ii];
llik[0] = llikBetanew;
// lpBeta[ii] = lpBetanew;
acc[ii] = acc[ii] + 1;
}else{ //if not accepted reset betanew[ii] with old beta
betanew[ii] = beta[ii];
}
}
delete[] betanew;
}
//update intercept
void updateInterceptFixedEF(double *X,double *Y,double *Z,int *T,int *nn,int *pp,int *dd,int *KK, double *beta,double *intercept,double *alpha,double *mu,double *sigmasq,double *tuneInt,double *lpInt,double *llik,double *acc){
double intnew;
double lpIntnew;
double llikIntnew;
intnew = intercept[0] + tuneInt[0]*rnorm(0.0,1.0);
LogpriorBeta(&intnew, mu,sigmasq,&lpIntnew);
AllLogLik(X, Y, Z, T, nn, pp, dd,KK, beta, &intnew,alpha,&llikIntnew);
double logratio = lpIntnew-lpInt[0]+llikIntnew-llik[0];
if(log(runif(0.0,1.0)) < logratio){
intercept[0] = intnew;
llik[0] = llikIntnew;
lpInt[0] = lpIntnew;
acc[0] = acc[0] + 1;
}
}
/////
//update alpha
void updateAlpha(double *X,double *Y,double *Z,int *T,int *nn,int *pp,int *dd,int *KK,double *beta, double *intercept,double *alpha,double *mu,double *sigmasq,double *tuneAlpha,double *lpAlpha,double *llikAll,double *accalpha)
{
double Alphanew = alpha[0] + tuneAlpha[0]*rnorm(0.0,1.0);
double llikAlphanew;
double lpAlphanew;
AllLogLik(X, Y, Z, T, nn, pp, dd,KK, beta, intercept,&Alphanew,&llikAlphanew);
LogpriorAlpha(&Alphanew, mu, sigmasq,&lpAlphanew);
double logratio = llikAlphanew - llikAll[0] + lpAlphanew -lpAlpha[0];
if(log(runif(0.0,1.0)) < logratio){
alpha[0] = Alphanew;
lpAlpha[0] = lpAlphanew;
accalpha[0] = accalpha[0] + 1;
}
}
/////////////////////////////////
//SAMPLER FOR FIXED EFFECT MODEL
extern "C" {
void sampleFixedIntervention(int *niter, double *XX,double *YY,double *ZZ,int *TT,int *nn,int *PP,int *dd,int *KK,double *beta,double *intercept,double *alpha,double *MuAlpha,double *SigmaAlpha,double *MuBeta, double *SigmaBeta, double *MuZ,double *VarZ, double *tuneBetaAll,double *tuneInt,double *tuneAlpha,double *tuneZAll,double *accBetaAll,double *accAlpha, double *accIntAll,double *accZAll, double *betaFinal, double *AlphaFinal, double *ZZFinal, double *InterceptFinal,double *Zvar1,double *Zvar2,double *likelihood,double *PriorA,double *PriorB, int *intervention){
int slen = KK[0]+1;
double* sumn = 0; // need this to store the updated parameters
sumn = new double[slen];
sumn[0] = 0;
for(int iz =0; iz < KK[0]; iz++){
sumn[iz+1] = sumn[iz] + nn[iz];
}
int ll = KK[0];
int sumAll = sumn[ll];
// int lenY = sumn[KK[0]+1];
double lpAlpha,lpInt,llik;
// double lpZ, lpBeta;
double* D = 0;
D = new double[dd[0]];
// double D[dd[0]];
// double C;
double muInt = MuBeta[PP[0]];
double sigmaInt = SigmaBeta[PP[0]];
for(int ii = 0; ii < niter[0]; ii++){
GetRNGstate();
LogpriorAlpha(alpha, MuAlpha, SigmaAlpha,&lpAlpha);
//Recursion for KK groups
double llikall = 0.0;
for(int ss = 0.0;ss < dd[0];ss++){
D[ss] = 0.0;
}
for(int kk = 0; kk < KK[0]; kk++){
int ss = sumn[kk];
// int ss2 = sumn[kk + 1];
double* accZ = 0;
accZ = new double[nn[kk]];
double* tuneZ = 0;
tuneZ = new double[nn[kk]];
double* XMat = 0;
XMat = new double[nn[kk]*nn[kk]*PP[0]];
double* YMat = 0;
YMat = new double[nn[kk]*nn[kk]];
double* ZMat = 0;
ZMat = new double[nn[kk]*dd[0]];
//read X, Y, X for a network
getY(YY,YMat,nn,kk);
readX(XX,XMat,nn,PP[0],kk);
getZ(ZZ,ZMat,nn,dd[0],kk);
for(int zz =0;zz < nn[kk];zz++){
accZ[zz] = accZAll[zz+ss];
tuneZ[zz] = tuneZAll[zz+ss];
}
//loglikelihood for the group
FullLogLik(beta, YMat, XMat, ZMat,alpha, &TT[kk],intercept, &nn[kk],PP,dd,&llik);
updateZ(XMat,YMat,ZMat,&TT[kk],&nn[kk],PP,dd,beta,intercept,alpha,MuZ,VarZ,tuneZ,&llik,accZ);
//store
for(int nZ=0; nZ < nn[kk]; nZ++){
for(int dZ=0; dZ < dd[0]; dZ++){
ZZ[nZ+dZ*nn[kk]+dd[0]*ss] = ZMat[nZ+dZ*nn[kk]];
// ZZFinal[nZ+ss+ss2*dZ+ii*sumAll*dd[0]] = ZMat[nZ+dZ*nn[kk]];
ZZFinal[nZ+dZ*nn[kk]+dd[0]*ss+ii*sumAll*dd[0]] = ZMat[nZ+dZ*nn[kk]];
D[dZ] += square(ZMat[nZ+dZ*nn[kk]] - 0.0);
}
accZAll[nZ+ss] = accZ[nZ];
}
llikall += llik;
delete[] XMat;
delete[] YMat;
delete[] ZMat;
delete[] tuneZ;
delete[] accZ;
}
// Rprintf("%f",llikall);
LogpriorBeta(intercept, &muInt, &sigmaInt,&lpInt); //logprior of intercept
updateInterceptFixedEF(XX,YY,ZZ,TT,nn,PP,dd,KK,beta,intercept,alpha,&muInt,&sigmaInt,tuneInt,&lpInt,&llikall, accIntAll);
updateBetamultiFixedEF(XX,YY,ZZ,TT,nn,PP,dd,KK, beta,intercept,alpha,MuBeta,SigmaBeta,tuneBetaAll,&llikall, accBetaAll);
//store the updated values
InterceptFinal[ii] = intercept[0];
for(int pp = 0; pp < PP[0]; pp++){
betaFinal[ii+pp*niter[0]] = beta[pp];
}
//update alpha
if(intervention[0] == 1){
updateAlpha(XX,YY,ZZ,TT,nn,PP,dd,KK,beta,intercept,alpha,MuAlpha,SigmaAlpha,tuneAlpha,&lpAlpha,&llikall,accAlpha);
AlphaFinal[ii] = alpha[0];
}
//Update varaince for ZZ
/*
for(int vv = 0; vv < dd[0]; vv++){
D[vv] = D[vv]/2.0 + PriorB[0];
C = PriorA[0] + (sumAll)/2.0;
VarZ[vv] = 1.0/rgamma(C,1.0/D[vv]);
}
*/
Zvar1[ii] = VarZ[0];
Zvar2[ii] = VarZ[1];
likelihood[ii] = llikall;
PutRNGstate();
}
delete[] sumn;
delete[] D;
}
/////////////////////////////////////////////////////
///////////////SAMLPE FOR RANDOM EFFECT MODEL////////
//////////////////////////////////////////////////////
void sampleRandomIntervention(int *niter, double *XX,double *YY,double *ZZ,int *TT,int *nn,int *PP,int *dd,int *KK,double *beta,double *intercept,double *alpha,double *MuAlpha,double *SigmaAlpha,double *MuBeta, double *SigmaBeta, double *MuZ,double *VarZ, double *tuneBetaAll,double *tuneInt,double *tuneAlpha,double *tuneZAll,double *accBetaAll,double *accAlpha, double *accIntAll,double *accZAll, double *betaFinal, double *AlphaFinal, double *ZZFinal, double *InterceptFinal,double *Zvar1,double *Zvar2,double *postVar, double *postMu, double *likelihood,double *PriorA,double *PriorB, int *intervention){
//Set starting points
// InitialVal(XX,YY,ZZ,nn,TT,dd,PP,KK,beta,intercept,alpha,intervention);
int slen = KK[0]+1;
double* sumn = 0; // need this to store the updated parameters
sumn = new double[slen];
sumn[0] = 0;
for(int iz =0; iz < KK[0]; iz++){
sumn[iz+1] = sumn[iz] + nn[iz];
}
int ll = KK[0];
int sumAll = sumn[ll];
// int lenY = sumn[KK[0]+1];
double lpAlpha,lpInt,llik;
//double lpZ, lpBeta;
double* accbeta = 0;
accbeta = new double[PP[0]];
double accInt = 0;
double* tuneBeta = 0;
tuneBeta = new double[PP[0]];
double* D = 0;
D = new double[dd[0]];
// double D[dd[0]];
// double C;
for(int ii = 0; ii < niter[0]; ii++){
double muInt = MuBeta[PP[0]];
double sigmaInt = SigmaBeta[PP[0]];
GetRNGstate();
LogpriorAlpha(alpha, MuAlpha, SigmaAlpha,&lpAlpha);
//Recursion for KK groups
double llikall = 0.0;
for(int ss = 0.0;ss < dd[0];ss++){
D[ss] = 0.0;
}
for(int kk = 0; kk < KK[0]; kk++){
int ss = sumn[kk];
// int ss2 = sumn[kk + 1];
double* accZ = 0;
accZ = new double[nn[kk]];
double* tuneZ = 0;
tuneZ = new double[nn[kk]];
double* XMat = 0;
XMat = new double[nn[kk]*nn[kk]*PP[0]];
double* YMat = 0;
YMat = new double[nn[kk]*nn[kk]];
double* ZMat = 0;
ZMat = new double[nn[kk]*dd[0]];
double* betaKK = 0;
betaKK = new double[PP[0]];
//read X, Y, X for a network
getY(YY,YMat,nn,kk);
readX(XX,XMat,nn,PP[0],kk);
getZ(ZZ,ZMat,nn,dd[0],kk);
//make sure beta has pp rows and KK columns
for(int ll = 0; ll<PP[0];ll++){
betaKK[ll] = beta[ll+kk*PP[0]];
accbeta[ll] = accBetaAll[ll+kk*PP[0]];
tuneBeta[ll] = tuneBetaAll[ll+kk*PP[0]];
//betaKK[ll] = beta[kk+ll*KK[0]];
//accbeta[ll] = accBetaAll[kk+ll*KK[0]];
//tuneBeta[ll] = tuneBetaAll[kk+ll*KK[0]];
}
accInt = accIntAll[kk];
for(int zz =0;zz < nn[kk];zz++){
accZ[zz] = accZAll[zz+ss];
tuneZ[zz] = tuneZAll[zz+ss];
}
//loglikelihood for the group
FullLogLik(betaKK, YMat, XMat, ZMat,alpha, &TT[kk], &intercept[kk],&nn[kk], PP, dd,&llik);
LogpriorBeta(&intercept[kk], &muInt, &sigmaInt,&lpInt); //logprior of intercept
updateBetamulti(XMat,YMat,ZMat,&TT[kk],&nn[kk],PP,dd,betaKK,&intercept[kk],alpha,MuBeta,SigmaBeta,tuneBeta,&llik, accbeta);
updateZ(XMat,YMat,ZMat,&TT[kk],&nn[kk],PP,dd,betaKK,&intercept[kk],alpha,MuZ, VarZ,tuneZ,&llik,accZ);
updateIntercept(XMat,YMat,ZMat,&TT[kk],&nn[kk],PP,dd,betaKK,&intercept[kk],alpha,&muInt,&sigmaInt,&tuneInt[kk],&lpInt,&llik, &accInt);
//store the updated values
InterceptFinal[ii+kk*niter[0]] = intercept[kk];
for(int pp = 0; pp < PP[0]; pp++){
beta[pp+kk*PP[0]] = betaKK[pp];
betaFinal[ii+pp*niter[0]+kk*PP[0]*niter[0]] = betaKK[pp];
accBetaAll[pp+kk*PP[0]] = accbeta[pp];
}
for(int nZ=0; nZ < nn[kk]; nZ++){
for(int dZ=0; dZ < dd[0]; dZ++){
ZZ[nZ+dZ*nn[kk]+dd[0]*ss] = ZMat[nZ+dZ*nn[kk]];
// ZZFinal[nZ+ss+ss2*dZ+ii*sumAll*dd[0]] = ZMat[nZ+dZ*nn[kk]];
ZZFinal[nZ+dZ*nn[kk]+dd[0]*ss+ii*sumAll*dd[0]] = ZMat[nZ+dZ*nn[kk]];
D[dZ] += square(ZMat[nZ+dZ*nn[kk]] - 0.0);
}
accZAll[nZ+ss] = accZ[nZ];
}
accIntAll[kk] = accInt;
llikall += llik;
delete[] XMat;
delete[] YMat;
delete[] ZMat;
delete[] betaKK;
delete[] tuneZ;
delete[] accZ;
}
//update alpha
// Rprintf("%f",llikall);
if(intervention[0] == 1){
updateAlpha(XX,YY,ZZ,TT,nn,PP,dd,KK,beta,intercept,alpha,MuAlpha,SigmaAlpha,tuneAlpha,&lpAlpha,&llikall,accAlpha);
AlphaFinal[ii] = alpha[0];
}
//update mu and sigmasq for beta
for(int ee = 0;ee< PP[0];ee++){
double denom = SigmaBeta[ee] + KK[0];
double sumbeta = 0.0;
for(int ww = 0;ww < KK[0];ww++){
sumbeta += beta[ee+ww*PP[0]];
}
double tempmu = (sumbeta + SigmaBeta[ee]*MuBeta[ee])/denom;
double tempssq = SigmaBeta[ee]/denom;
MuBeta[ee] = rnorm(tempmu,tempssq);
postMu[ii+niter[0]*ee] = MuBeta[ee];
}
//Also update mu for intercept
double denom = SigmaBeta[PP[0]] + KK[0];
double sumbeta = 0.0;
for(int ww = 0; ww < KK[0]; ww++){
sumbeta += intercept[ww];
}
double tempmu = (sumbeta + SigmaBeta[PP[0]]*MuBeta[PP[0]])/denom;
double tempssq = SigmaBeta[PP[0]]/denom;
MuBeta[PP[0]] = rnorm(tempmu,tempssq);
postMu[ii+niter[0]*PP[0]] = MuBeta[PP[0]];
//Update variance for slopes
for(int ee = 0; ee < PP[0]; ee++){
double B = 0.0;
for(int vv = 0;vv<KK[0];vv++){
B += square(beta[ee+vv*PP[0]] - MuBeta[ee]);
}
B = B/2.0 + PriorB[0];
double A = PriorA[0] + KK[0]/2.0;
SigmaBeta[ee] = 1.0/rgamma(A,1/B);
postVar[ii+niter[0]*ee] = SigmaBeta[ee];
}
//Also update variance for intercept
double B = 0.0;
for(int vv = 0; vv < KK[0]; vv++){
B += square(intercept[vv] - MuBeta[PP[0]]);
}
B = B/2.0 + PriorB[0];
double A = PriorA[0] + KK[0]/2.0;
SigmaBeta[PP[0]] = 1.0/rgamma(A,1/B);
postVar[ii+niter[0]*PP[0]] = SigmaBeta[PP[0]];
//Update varaince for ZZ
/*
for(int vv = 0; vv < dd[0]; vv++){
D[vv] = D[vv]/2.0 + PriorB[0];
C = PriorA[0] + (sumAll)/2.0;
VarZ[vv] = 1.0/rgamma(C,1.0/D[vv]);
}
*/
Zvar1[ii] = VarZ[0];
Zvar2[ii] = VarZ[1];
likelihood[ii] = llikall;
PutRNGstate();
}
delete[] accbeta;
delete[] tuneBeta;
delete[] sumn;
delete[] D;
}
}