/* MICROARRAY NORMALIZATION-FINDING THE CONTROLS-used in the manucript ``A method for normalizing microarrays using the genes that are not differentially expresed'' Cavan Reilly, Ph.D. Assistant Professor Division of Biostatistics University of Minnesota e-mail: cavanr@biostat.umn.edu STATISTICAL MODEL This fits the gene normalization model with replicated arrays and uses hierarchical distributions for the indicators and the scan specific effects 1st condition x_1ij ~ N(phi_1+theta_1i+psi_1i(eta_i-theta_1i),sigma^2_1i) y_1ij ~ N(phi_2+eta_i,tau^2_1i) x_2ij ~ N(phi_3+theta_1i+psi_1i(eta_i-theta_1i),sigma^2_2i) y_2ij ~ N(phi_4+eta_i,tau^2_2i) 2nd condition x_3ij ~ N(phi_5+theta_2i+psi_2i(eta_i-theta_2i),sigma^2_3i) y_3ij ~ N(phi_6+eta_i,tau^2_3i) x_4ij ~ N(phi_7+theta_2i+psi_2i(eta_i-theta_2i),sigma^2_4i) y_4ij ~ N(phi_8+eta_i,tau^2_4i) 3rd condition x_5ij ~ N(phi_9+theta_3i+psi_3i(eta_i-theta_3i),sigma^2_5i) y_5ij ~ N(phi_10+eta_i,tau^2_5i) x_6ij ~ N(phi_11+theta_3i+psi_3i(eta_i-theta_2i),sigma^2_6i) y_6ij ~ N(phi_12+eta_i,tau^2_6i) and the structure for channel 1 scans phi_i ~N(Phi1,zeta^2) for channel 2 scans phi_i ~N(Phi2,zeta^2) psi_ki ~ Bern(pi_k) then flat priors for everything else via the Gibbs sampler COMPUTATIONAL CONSIDERATIONS I have an S-plus (or R) function that takes the lists that store the data (one list for each array), finds initial values, and writes data and initial values to file. These lists have original intensity, background correction, background corrected intensities, etc. The 8th element of a list is the experimental sample channel and the 7th element is the reference sample. Here is the S-plus (or R) function > geneNrmRepInit function(xx1, xx2, xx3, xx4, xx5, xx6) { x1 <- log(matrix(xx1[[8]], byrow = T, 139, 3) + 3, base = 2) y1 <- log(matrix(xx1[[7]], byrow = T, 139, 3) + 3, base = 2) x2 <- log(matrix(xx2[[8]], byrow = T, 139, 3) + 3, base = 2) y2 <- log(matrix(xx2[[7]], byrow = T, 139, 3) + 3, base = 2) x3 <- log(matrix(xx3[[8]], byrow = T, 139, 3) + 3, base = 2) y3 <- log(matrix(xx3[[7]], byrow = T, 139, 3) + 3, base = 2) x4 <- log(matrix(xx4[[8]], byrow = T, 139, 3) + 3, base = 2) y4 <- log(matrix(xx4[[7]], byrow = T, 139, 3) + 3, base = 2) x5 <- log(matrix(xx5[[8]], byrow = T, 139, 3) + 3, base = 2) y5 <- log(matrix(xx5[[7]], byrow = T, 139, 3) + 3, base = 2) x6 <- log(matrix(xx6[[8]], byrow = T, 139, 3) + 3, base = 2) y6 <- log(matrix(xx6[[7]], byrow = T, 139, 3) + 3, base = 2) pp <- 0.1 ##calculate variance matrix: sig1, sig2, tau1 and tau2 sig1 <- sqrt(apply(x1, 1, var)) sig2 <- sqrt(apply(x2, 1, var)) sig3 <- sqrt(apply(x3, 1, var)) sig4 <- sqrt(apply(x4, 1, var)) sig5 <- sqrt(apply(x5, 1, var)) sig6 <- sqrt(apply(x6, 1, var)) tau1 <- sqrt(apply(y1, 1, var)) tau2 <- sqrt(apply(y2, 1, var)) tau3 <- sqrt(apply(y3, 1, var)) tau4 <- sqrt(apply(y4, 1, var)) tau5 <- sqrt(apply(y5, 1, var)) tau6 <- sqrt(apply(y6, 1, var)) sig1 <- ifelse(sig1 < 0.01, 0.1, sig1) sig2 <- ifelse(sig2 < 0.01, 0.1, sig2) sig3 <- ifelse(sig3 < 0.01, 0.1, sig3) sig4 <- ifelse(sig4 < 0.01, 0.1, sig4) sig5 <- ifelse(sig5 < 0.01, 0.1, sig5) sig6 <- ifelse(sig6 < 0.01, 0.1, sig6) tau1 <- ifelse(tau1 < 0.01, 0.1, tau1) tau2 <- ifelse(tau2 < 0.01, 0.1, tau2) tau3 <- ifelse(tau3 < 0.01, 0.1, tau3) tau4 <- ifelse(tau4 < 0.01, 0.1, tau4) tau5 <- ifelse(tau5 < 0.01, 0.1, tau5) tau6 <- ifelse(tau6 < 0.01, 0.1, tau6) phi4 <- 0 phi2 <- mean(c(y1)) - mean(c(y2)) phi8 <- 0 phi6 <- mean(c(y3)) - mean(c(y4)) phi12 <- 0 phi10 <- mean(c(y5)) - mean(c(y6)) eta <- 0.5 * (apply(y1, 1, mean) + apply(y2, 1, mean) - phi2 + apply(y3, 1, mean) + apply(y4, 1, mean) - phi6 + apply( y5, 1, mean) + apply(y6, 1, mean) - phi12) psi1 <- getphi1(eta, 0.5 * (x1 + x2), pp) psi2 <- getphi1(eta, 0.5 * (x3 + x4), pp) psi3 <- getphi1(eta, 0.5 * (x5 + x6), pp) phi3 <- 0 phi1 <- mean(apply(x1, 1, mean)[psi1 == 1] - apply(x2, 1, mean)[psi1 == 1]) phi5 <- 0 phi7 <- mean(apply(x3, 1, mean)[psi2 == 1] - apply(x4, 1, mean)[psi2 == 1]) phi9 <- 0 phi11 <- mean(apply(x5, 1, mean)[psi3 == 1] - apply(x6, 1, mean)[psi3 == 1]) theta1 <- ifelse(psi1 == 1, eta, 0.5 * (apply(x1, 1, mean) - phi1 + apply(x2, 1, mean) - phi5)) theta2 <- ifelse(psi2 == 1, eta, 0.5 * (apply(x3, 1, mean) - phi3 + apply(x4, 1, mean) - phi7)) theta3 <- ifelse(psi3 == 1, eta, 0.5 * (apply(x5, 1, mean) - phi5 + apply(x6, 1, mean) - phi9)) phi <- mean(c(phi1, phi2, phi3, phi4, phi5, phi6, phi7, phi8, phi9, phi10, phi11, phi12)) zeta <- sqrt(var(c(phi1, phi2, phi3, phi4, phi5, phi6, phi7, phi8, phi9, phi10, phi11, phi12))) pi1 <- mean(psi1) pi2 <- mean(psi2) pi3 <- mean(psi3) init <- c(theta1, theta2, theta3, eta, psi1, psi2, psi3, sig1, sig2, sig3, sig4, sig5, sig6, tau1, tau2, tau3, tau4, tau5, tau6, phi1, phi2, phi3, phi4, phi5, phi6, phi7, phi8, phi9, phi10, phi11, phi12, phi, zeta, pi1, pi2, pi3) write(init, file = "initValu.dat", ncol = 1) write.table(x1, "arr1c1.dat", sep = " ") write.table(y1, "arr1c2.dat", sep = " ") write.table(x2, "arr2c1.dat", sep = " ") write.table(y2, "arr2c2.dat", sep = " ") write.table(x3, "arr3c1.dat", sep = " ") write.table(y3, "arr3c2.dat", sep = " ") write.table(x4, "arr4c1.dat", sep = " ") write.table(y4, "arr4c2.dat", sep = " ") write.table(x5, "arr5c1.dat", sep = " ") write.table(y5, "arr5c2.dat", sep = " ") write.table(x6, "arr6c1.dat", sep = " ") write.table(y6, "arr6c2.dat", sep = " ") } > to use just type > geneNrmRepInit(NE4h1,NE4h2,VAC4h1,VAC4h2,VR4h1,VR4h2) then run the following C program. While the C++ code should be portable, there are a few compiler dependencies (I have debugged this with Borland C++ Builder standard version 3.0). These are 1. the "pragma" statements 2. the "getch" commands 3. include statements for the files condefs and conio (I think) If running on a UNIX platform get rid of these statements. The author may be able to assist if further difficulties are encountered. Contact at above e-mail address. */ #include #include #include #include #include /* next 2 lines are compiler dependent */ #pragma hdrstop #pragma argsused #define FREE_ARG char* static double dsqrarg; #define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg) /* program specific routines */ void gibbsSmplr(long setSeed,int noChn,int numIter,int numburn,int numpar, double **initValu,double ***simres); void sqrtR(double *stat,int noChn,int numiter,int numpar,double ***tens); /* matrix routines */ double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl, long ndh); void free_d3tensor(double ***t,long nrl,long nrh,long ncl,long nch, long ndl,long ndh); double **dmatrix(long nrl,long nrh,long ncl,long nch); void free_dmatrix(double **m, long nrl,long nrh,long ncl,long nch); double *dvector(long nl,long nh); void free_dvector(double *v,long nl,long nh); void nrerror(char error_text[]); /* random number routines */ double invGamma(double scal,int df, long *setSeed); double gasdev(long *idum); double ran1(long *idum); double gamdev(double ia, long *idum); double betadev(double a,double b, long *idum); /* reads in initial values */ void initlz(double **initValu,int noChn); /* nGlob is number of transcripts nopar is number of parameters */ int nGlob=139,numRep=3; int noPar=19*nGlob+17; /* p1Glob adn p2Glob are bounds on the proportion of unexpresed gene, i.e. the proportion of control genes */ double alphaGlob=1.0,betaGlob=1.0; double **x1Glob,**x2Glob,**x3Glob,**x4Glob,**x5Glob,**x6Glob,**y1Glob, **y2Glob,**y3Glob,**y4Glob,**y5Glob,**y6Glob; int main() { /* this file holds responses for array 1 channel 1*/ ifstream filer1("c:\\Spluswin\\home\\arr1c1.dat"); if(!filer1){ cerr << "Failure to open filer1\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 1 channel 2*/ ifstream filer2("c:\\Spluswin\\home\\arr1c2.dat"); if(!filer2){ cerr << "Failure to open filer2\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 2 channel 1*/ ifstream filer3("c:\\Spluswin\\home\\arr2c1.dat"); if(!filer3){ cerr << "Failure to open filer3\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 2 channel 2*/ ifstream filer4("c:\\Spluswin\\home\\arr2c2.dat"); if(!filer4){ cerr << "Failure to open filer4\n"; getch(); exit(EXIT_FAILURE); } /* note: array 3 is a replicate of array 1 and 4 of 2 */ /* this file holds responses for array 3 channel 1*/ ifstream filer5("c:\\Spluswin\\home\\arr3c1.dat"); if(!filer5){ cerr << "Failure to open filer5\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 3 channel 2*/ ifstream filer6("c:\\Spluswin\\home\\arr3c2.dat"); if(!filer6){ cerr << "Failure to open filer6\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 4 channel 1*/ ifstream filer7("c:\\Spluswin\\home\\arr4c1.dat"); if(!filer7){ cerr << "Failure to open filer7\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 4 channel 2*/ ifstream filer8("c:\\Spluswin\\home\\arr4c2.dat"); if(!filer8){ cerr << "Failure to open filer8\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 5 channel 1*/ ifstream filer9("c:\\Spluswin\\home\\arr5c1.dat"); if(!filer9){ cerr << "Failure to open filer9\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 5 channel 2*/ ifstream filer10("c:\\Spluswin\\home\\arr5c2.dat"); if(!filer10){ cerr << "Failure to open filer10\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 6 channel 1*/ ifstream filer11("c:\\Spluswin\\home\\arr6c1.dat"); if(!filer11){ cerr << "Failure to open filer11\n"; getch(); exit(EXIT_FAILURE); } /* this file holds responses for array 6 channel 2*/ ifstream filer12("c:\\Spluswin\\home\\arr6c2.dat"); if(!filer12){ cerr << "Failure to open filer12\n"; getch(); exit(EXIT_FAILURE); } /* samples from posterior distribution */ ofstream filew1("c:\\Spluswin\\home\\arrOut1.dat"); if(!filew1){ cerr << "Failure to open filew1\n"; getch(); exit(EXIT_FAILURE); } /* samples from posterior distribution */ ofstream filew2("c:\\Spluswin\\home\\arrOut2.dat"); if(!filew2){ cerr << "Failure to open filew2\n"; getch(); exit(EXIT_FAILURE); } /* samples from posterior distribution */ ofstream filew3("c:\\Spluswin\\home\\arrOut3.dat"); if(!filew3){ cerr << "Failure to open filew3\n"; getch(); exit(EXIT_FAILURE); } /* numIter is numiter per chain in metrop, numburn is number of burn in iterations, noChn is number of chains to use in MCMC and trimfac inicates how much thinning of chain output use 1 to have none */ int numIter=2000,numburn=1000,noChn=1,thinfac=1,i,j,k; long setSeed=-1273; double ***simres; x1Glob=dmatrix(1,nGlob,1,numRep); x2Glob=dmatrix(1,nGlob,1,numRep); x3Glob=dmatrix(1,nGlob,1,numRep); x4Glob=dmatrix(1,nGlob,1,numRep); x5Glob=dmatrix(1,nGlob,1,numRep); x6Glob=dmatrix(1,nGlob,1,numRep); y1Glob=dmatrix(1,nGlob,1,numRep); y2Glob=dmatrix(1,nGlob,1,numRep); y3Glob=dmatrix(1,nGlob,1,numRep); y4Glob=dmatrix(1,nGlob,1,numRep); y5Glob=dmatrix(1,nGlob,1,numRep); y6Glob=dmatrix(1,nGlob,1,numRep); simres=d3tensor(1,noChn,1,numIter-numburn,1,3*nGlob); /* read in array 1 channel 1*/ filer1.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer1 >> x1Glob[i][j]; } } /* read in array 1 channel 2 */ filer2.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer2 >> y1Glob[i][j]; } } /* read in array 2 channel 1*/ filer3.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer3 >> x2Glob[i][j]; } } /* read in array 2 channel 2 */ filer4.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer4 >> y2Glob[i][j]; } } /* read in array 3 channel 1*/ filer5.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer5 >> x3Glob[i][j]; } } /* read in array 3 channel 2 */ filer6.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer6 >> y3Glob[i][j]; } } /* read in array 4 channel 1*/ filer7.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer7 >> x4Glob[i][j]; } } /* read in array 4 channel 2 */ filer8.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer8 >> y4Glob[i][j]; } } /* read in array 5 channel 1*/ filer9.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer9 >> x5Glob[i][j]; } } /* read in array 5 channel 2 */ filer10.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer10 >> y5Glob[i][j]; } } /* read in array 6 channel 1*/ filer11.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer11 >> x6Glob[i][j]; } } /* read in array 6 channel 2 */ filer12.seekg(0); for(i=1;i<=nGlob;i++){ for(j=1;j<=numRep;j++){ filer12 >> y6Glob[i][j]; } } double **initValu; initValu=dmatrix(1,noPar,1,noChn); /* initial values moved to the end of the program-read in from a file */ initlz(initValu,noPar); /* run the gibbs sampler */ gibbsSmplr(setSeed,noChn,numIter,numburn,noPar,initValu,simres); /* write simulations for fold changes to file */ for(i=1;i<=noChn;i++){ for(j=1;j<=numIter-numburn;j++){ if(j%thinfac==0){ for(k=1;k<=nGlob;k++) filew1 << simres[i][j][k] << " "; for(k=nGlob+1;k<=2*nGlob;k++) filew2 << simres[i][j][k] << " "; for(k=2*nGlob+1;k<=3*nGlob;k++) filew3 << simres[i][j][k] << " "; } } } free_dmatrix(initValu,1,noPar,1,noChn); free_dmatrix(x1Glob,1,nGlob,1,numRep); free_dmatrix(x2Glob,1,nGlob,1,numRep); free_dmatrix(x3Glob,1,nGlob,1,numRep); free_dmatrix(x4Glob,1,nGlob,1,numRep); free_dmatrix(x5Glob,1,nGlob,1,numRep); free_dmatrix(x6Glob,1,nGlob,1,numRep); free_dmatrix(y1Glob,1,nGlob,1,numRep); free_dmatrix(y2Glob,1,nGlob,1,numRep); free_dmatrix(y3Glob,1,nGlob,1,numRep); free_dmatrix(y4Glob,1,nGlob,1,numRep); free_dmatrix(y5Glob,1,nGlob,1,numRep); free_dmatrix(y6Glob,1,nGlob,1,numRep); free_d3tensor(simres,1,noChn,1,numIter-numburn,1,3*nGlob); } void gibbsSmplr(long setSeed,int noChn,int numIter,int numburn,int numpar, double **initValu,double ***simres) { /* runs the Gibbs' sampler-simres holds sim. results in form chain by iteration by parameter, initValu is matrix of initial values for the chain this finds convergence diagnostic stat and reports it-stops if converrgence taken place */ int ii,jj,i,j,k; /* set init=1 if you want all simulations on some iteration written to file */ int init=0; double RSTOP=1.01,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,mm,vv,temp, odds,maxRstat,p1,*pi1,*pi2,*pi3,*phiHier1,*phiHier2,*zeta,**theta1,**theta2,**theta3,**eta, **phi,**psi1,**psi2,**psi3,**sig1,**sig2,**sig3,**sig4,**sig5,**sig6, **tau1,**tau2,**tau3,**tau4,**tau5,**tau6,*Rstat,***fldchng; double nu0=1.0,sig0=2.5,tau0=2.5; /* samples from posterior distribution */ ofstream filew1g("c:\\Spluswin\\home\\arrOut1g.dat"); if(!filew1g){ cerr << "Failure to open filew1g\n"; getch(); exit(EXIT_FAILURE); } pi1=dvector(1,noChn); pi2=dvector(1,noChn); pi3=dvector(1,noChn); phiHier1=dvector(1,noChn); phiHier2=dvector(1,noChn); zeta=dvector(1,noChn); theta1=dmatrix(1,nGlob,1,noChn); theta2=dmatrix(1,nGlob,1,noChn); theta3=dmatrix(1,nGlob,1,noChn); eta=dmatrix(1,nGlob,1,noChn); psi1=dmatrix(1,nGlob,1,noChn); psi2=dmatrix(1,nGlob,1,noChn); psi3=dmatrix(1,nGlob,1,noChn); sig1=dmatrix(1,nGlob,1,noChn); sig2=dmatrix(1,nGlob,1,noChn); sig3=dmatrix(1,nGlob,1,noChn); sig4=dmatrix(1,nGlob,1,noChn); sig5=dmatrix(1,nGlob,1,noChn); sig6=dmatrix(1,nGlob,1,noChn); tau1=dmatrix(1,nGlob,1,noChn); tau2=dmatrix(1,nGlob,1,noChn); tau3=dmatrix(1,nGlob,1,noChn); tau4=dmatrix(1,nGlob,1,noChn); tau5=dmatrix(1,nGlob,1,noChn); tau6=dmatrix(1,nGlob,1,noChn); phi=dmatrix(1,12,1,noChn); Rstat=dvector(1,3*nGlob); fldchng=d3tensor(1,noChn,1,numIter-numburn,1,3*nGlob); /* map initial values to theta1,theta2,eta,psi1,etc- and put an index on all these quantities so we can run chains in parallel */ for(ii=1;ii<=noChn;ii++){ for(j=1;j<=numpar;j++){ if(j<=nGlob) theta1[j][ii]=initValu[j][ii]; if((nGlob+1 <= j) && (j <= 2*nGlob)) theta2[j-nGlob][ii]=initValu[j][ii]; if((2*nGlob+1 <= j) && (j <= 3*nGlob)) theta3[j-2*nGlob][ii]=initValu[j][ii]; if((3*nGlob+1 <= j) & (j <= 4*nGlob)) eta[j-3*nGlob][ii]=initValu[j][ii]; if((4*nGlob+1 <= j) && (j <= 5*nGlob)) psi1[j-4*nGlob][ii]=initValu[j][ii]; if((5*nGlob+1 <= j) && (j <= 6*nGlob)) psi2[j-5*nGlob][ii]=initValu[j][ii]; if((6*nGlob+1 <= j) && (j <= 7*nGlob)) psi3[j-6*nGlob][ii]=initValu[j][ii]; if((7*nGlob+1 <= j) && (j <= 8*nGlob)) sig1[j-7*nGlob][ii]=initValu[j][ii]; if((8*nGlob+1 <= j) && (j <= 9*nGlob)) sig2[j-8*nGlob][ii]=initValu[j][ii]; if((9*nGlob+1 <= j) && (j <= 10*nGlob)) sig3[j-9*nGlob][ii]=initValu[j][ii]; if((10*nGlob+1 <= j) && (j <= 11*nGlob)) sig4[j-10*nGlob][ii]=initValu[j][ii]; if((11*nGlob+1 <= j) && (j <= 12*nGlob)) sig5[j-11*nGlob][ii]=initValu[j][ii]; if((12*nGlob+1 <= j) && (j <= 13*nGlob)) sig6[j-12*nGlob][ii]=initValu[j][ii]; if((13*nGlob+1 <= j) && (j <= 14*nGlob)) tau1[j-13*nGlob][ii]=initValu[j][ii]; if((14*nGlob+1 <= j) && (j <= 15*nGlob)) tau2[j-14*nGlob][ii]=initValu[j][ii]; if((15*nGlob+1 <= j) && (j <= 16*nGlob)) tau3[j-15*nGlob][ii]=initValu[j][ii]; if((16*nGlob+1 <= j) && (j <= 17*nGlob)) tau4[j-16*nGlob][ii]=initValu[j][ii]; if((17*nGlob+1 <= j) && (j <= 18*nGlob)) tau5[j-17*nGlob][ii]=initValu[j][ii]; if((18*nGlob+1 <= j) && (j <= 19*nGlob)) tau6[j-18*nGlob][ii]=initValu[j][ii]; if((19*nGlob+1 <= j) && (j <= 19*nGlob+12)) phi[j-19*nGlob][ii]=initValu[j][ii]; } phiHier1[ii]=phiHier2[ii]=initValu[19*nGlob+13][ii]; zeta[ii]=initValu[19*nGlob+14][ii]; /* use small initial values for these parameters */ pi1[ii]=0.01;//initValu[19*nGlob+15][ii]; pi2[ii]=0.01;//initValu[19*nGlob+16][ii]; pi3[ii]=0.01;//initValu[19*nGlob+17][ii]; } for(i=1;i<=numIter;i++){ // cout << "starting iteration number " << i << endl; for(ii=1;ii<=noChn;ii++){ // cout << "starting chain number " << ii << endl; /* sample theta1 given rest */ for(j=1;j<=nGlob;j++){ if(psi1[j][ii]==1.0) theta1[j][ii]=eta[j][ii]; else{ m1=0.0; for(k=1;k<=numRep;k++) m1+=x1Glob[j][k]-phi[1][ii]; m1/=DSQR(sig1[j][ii]); m2=0.0; for(k=1;k<=numRep;k++) m2+=x2Glob[j][k]-phi[3][ii]; m2/=DSQR(sig2[j][ii]); vv=1.0/(3.0*(1.0/DSQR(sig1[j][ii])+1.0/DSQR(sig2[j][ii]))); mm=(m1+m2)*vv; theta1[j][ii]=sqrt(vv)*gasdev(&setSeed)+mm; } } /* sample theta2 given rest */ for(j=1;j<=nGlob;j++){ if(psi2[j][ii]==1.0) theta2[j][ii]=eta[j][ii]; else{ m1=0.0; for(k=1;k<=numRep;k++) m1+=x3Glob[j][k]-phi[5][ii]; m1/=DSQR(sig3[j][ii]); m2=0.0; for(k=1;k<=numRep;k++) m2+=x4Glob[j][k]-phi[7][ii]; m2/=DSQR(sig4[j][ii]); vv=1.0/(3.0*(1.0/DSQR(sig3[j][ii])+1.0/DSQR(sig4[j][ii]))); mm=(m1+m2)*vv; theta2[j][ii]=sqrt(vv)*gasdev(&setSeed)+mm; } } /* sample theta3 given rest */ for(j=1;j<=nGlob;j++){ if(psi3[j][ii]==1.0) theta3[j][ii]=eta[j][ii]; else{ m1=0.0; for(k=1;k<=numRep;k++) m1+=x5Glob[j][k]-phi[9][ii]; m1/=DSQR(sig5[j][ii]); m2=0.0; for(k=1;k<=numRep;k++) m2+=x6Glob[j][k]-phi[11][ii]; m2/=DSQR(sig6[j][ii]); vv=1.0/(3.0*(1.0/DSQR(sig5[j][ii])+1.0/DSQR(sig6[j][ii]))); mm=(m1+m2)*vv; theta3[j][ii]=sqrt(vv)*gasdev(&setSeed)+mm; } } /* sample eta given rest*/ for(j=1;j<=nGlob;j++){ m1=m2=m3=m4=m5=m6=m7=m8=m9=m10=m11=m12=0.0; if(psi1[j][ii]==1.0){ for(k=1;k<=numRep;k++) m1+=x1Glob[j][k]-phi[1][ii]; m1/=DSQR(sig1[j][ii]); for(k=1;k<=numRep;k++) m3+=x2Glob[j][k]-phi[3][ii]; m3/=DSQR(sig2[j][ii]); } if(psi2[j][ii]==1.0){ for(k=1;k<=numRep;k++) m5+=x3Glob[j][k]-phi[5][ii]; m5/=DSQR(sig3[j][ii]); for(k=1;k<=numRep;k++) m7+=x4Glob[j][k]-phi[7][ii]; m7/=DSQR(sig4[j][ii]); } if(psi3[j][ii]==1.0){ for(k=1;k<=numRep;k++) m9+=x5Glob[j][k]-phi[9][ii]; m9/=DSQR(sig5[j][ii]); for(k=1;k<=numRep;k++) m11+=x6Glob[j][k]-phi[11][ii]; m11/=DSQR(sig6[j][ii]); } for(k=1;k<=numRep;k++) m2+=y1Glob[j][k]-phi[2][ii]; m2/=DSQR(tau1[j][ii]); for(k=1;k<=numRep;k++) m4+=y2Glob[j][k]-phi[4][ii]; m4/=DSQR(tau2[j][ii]); for(k=1;k<=numRep;k++) m6+=y3Glob[j][k]-phi[6][ii]; m6/=DSQR(tau3[j][ii]); for(k=1;k<=numRep;k++) m8+=y4Glob[j][k]-phi[8][ii]; m8/=DSQR(tau4[j][ii]); for(k=1;k<=numRep;k++) m10+=y5Glob[j][k]-phi[10][ii]; m10/=DSQR(tau5[j][ii]); for(k=1;k<=numRep;k++) m12+=y6Glob[j][k]-phi[12][ii]; m12/=DSQR(tau6[j][ii]); vv=numRep*(psi1[j][ii]*(1.0/DSQR(sig1[j][ii])+1.0/DSQR(sig2[j][ii]))+ psi2[j][ii]*(1.0/DSQR(sig3[j][ii])+1.0/DSQR(sig4[j][ii]))+ psi3[j][ii]*(1.0/DSQR(sig5[j][ii])+1.0/DSQR(sig6[j][ii]))+ 1.0/DSQR(tau1[j][ii])+1.0/DSQR(tau2[j][ii])+1.0/DSQR(tau3[j][ii])+ 1.0/DSQR(tau4[j][ii])+1.0/DSQR(tau5[j][ii])+1.0/DSQR(tau6[j][ii])); mm=(psi1[j][ii]*(m1+m3)+psi2[j][ii]*(m5+m7)+psi3[j][ii]*(m9+m11)+ m2+m4+m6+m8+m10+m12)/vv; eta[j][ii]=gasdev(&setSeed)/sqrt(vv)+mm; } /* sample phi given the rest */ mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=x1Glob[j][k]-theta1[j][ii]-psi1[j][ii]*(eta[j][ii]-theta1[j][ii]); } vv+=1.0/DSQR(sig1[j][ii]); temp/=DSQR(sig1[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier1[ii]/DSQR(zeta[ii]); mm/=vv; phi[1][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=x2Glob[j][k]-theta1[j][ii]-psi1[j][ii]*(eta[j][ii]-theta1[j][ii]); } vv+=1.0/DSQR(sig2[j][ii]); temp/=DSQR(sig2[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier1[ii]/DSQR(zeta[ii]); mm/=vv; phi[3][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=y1Glob[j][k]-eta[j][ii]; } vv+=1.0/DSQR(tau1[j][ii]); temp/=DSQR(tau1[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier2[ii]/DSQR(zeta[ii]); mm/=vv; phi[2][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=y2Glob[j][k]-eta[j][ii]; } vv+=1.0/DSQR(tau2[j][ii]); temp/=DSQR(tau2[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier2[ii]/DSQR(zeta[ii]); mm/=vv; phi[4][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=x3Glob[j][k]-theta2[j][ii]-psi2[j][ii]*(eta[j][ii]-theta2[j][ii]); } vv+=1.0/DSQR(sig3[j][ii]); temp/=DSQR(sig3[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier1[ii]/DSQR(zeta[ii]); mm/=vv; phi[5][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=x4Glob[j][k]-theta2[j][ii]-psi2[j][ii]*(eta[j][ii]-theta2[j][ii]); } vv+=1.0/DSQR(sig4[j][ii]); temp/=DSQR(sig4[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier1[ii]/DSQR(zeta[ii]); mm/=vv; phi[7][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=y3Glob[j][k]-eta[j][ii]; } vv+=1.0/DSQR(tau3[j][ii]); temp/=DSQR(tau3[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier2[ii]/DSQR(zeta[ii]); mm/=vv; phi[6][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=y4Glob[j][k]-eta[j][ii]; } vv+=1.0/DSQR(tau4[j][ii]); temp/=DSQR(tau4[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier2[ii]/DSQR(zeta[ii]); mm/=vv; phi[8][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=x5Glob[j][k]-theta3[j][ii]-psi3[j][ii]*(eta[j][ii]-theta3[j][ii]); } vv+=1.0/DSQR(sig5[j][ii]); temp/=DSQR(sig5[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier1[ii]/DSQR(zeta[ii]); mm/=vv; phi[9][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=x6Glob[j][k]-theta3[j][ii]-psi3[j][ii]*(eta[j][ii]-theta3[j][ii]); } vv+=1.0/DSQR(sig6[j][ii]); temp/=DSQR(sig6[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier1[ii]/DSQR(zeta[ii]); mm/=vv; phi[11][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=y5Glob[j][k]-eta[j][ii]; } vv+=1.0/DSQR(tau5[j][ii]); temp/=DSQR(tau5[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier2[ii]/DSQR(zeta[ii]); mm/=vv; phi[10][ii]=gasdev(&setSeed)/sqrt(vv)+mm; mm=0.0; vv=0.0; for(j=1;j<=nGlob;j++){ temp=0.0; for(k=1;k<=numRep;k++){ temp+=y6Glob[j][k]-eta[j][ii]; } vv+=1.0/DSQR(tau6[j][ii]); temp/=DSQR(tau6[j][ii]); mm+=temp; } vv*=numRep; vv+=1.0/DSQR(zeta[ii]); mm+=phiHier2[ii]/DSQR(zeta[ii]); mm/=vv; phi[12][ii]=gasdev(&setSeed)/sqrt(vv)+mm; /* sample phiHier */ // mm=0.0; // for(j=1;j<=12;j++) mm+=phi[j][ii]; // mm/=12.0; // phiHier[ii]=gasdev(&setSeed)*zeta[ii]/sqrt(12.0)+mm; /* sample hierarchical mean for scan effects for channel 1 */ mm=0.0; for(j=1;j<=6;j++) mm+=phi[2*j-1][ii]; mm/=6.0; phiHier1[ii]=gasdev(&setSeed)*zeta[ii]/sqrt(6.0)+mm; /* sample hierarchical mean for scan effects for channel 2*/ mm=0.0; for(j=1;j<=6;j++) mm+=phi[2*j][ii]; mm/=6.0; phiHier2[ii]=gasdev(&setSeed)*zeta[ii]/sqrt(6.0)+mm; /* sample zeta-hierarchical variance */ vv=0.0; for(j=1;j<=6;j++) vv+=DSQR(phi[2*j-1][ii]-phiHier1[ii]); for(j=1;j<=6;j++) vv+=DSQR(phi[2*j][ii]-phiHier2[ii]); zeta[ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+12.0),12.0+nu0,&setSeed)); /* sample psi1 given the rest */ for(j=1;j<=nGlob;j++){ m1=m2=m3=m4=0.0; for(k=1;k<=numRep;k++){ m1+=DSQR(x1Glob[j][k]-phi[1][ii]-eta[j][ii]); m2+=DSQR(x1Glob[j][k]-phi[1][ii]-theta1[j][ii]); m3+=DSQR(x2Glob[j][k]-phi[3][ii]-eta[j][ii]); m4+=DSQR(x2Glob[j][k]-phi[3][ii]-theta1[j][ii]); } m1/=DSQR(sig1[j][ii]); m2/=DSQR(sig1[j][ii]); m3/=DSQR(sig2[j][ii]); m4/=DSQR(sig2[j][ii]); temp=-0.5*(m1-m2+m3-m4); if(temp>300.0) temp=300.0; odds=exp(temp)*pi1[ii]/(1-pi1[ii]); p1=odds/(1.0+odds); if(ran1(&setSeed)300.0) temp=300.0; odds=exp(temp)*pi2[ii]/(1-pi2[ii]); p1=odds/(1.0+odds); if(ran1(&setSeed)300.0) temp=300.0; odds=exp(temp)*pi3[ii]/(1-pi3[ii]); p1=odds/(1.0+odds); if(ran1(&setSeed) numburn){ // filew1g << pi1[ii] << " "; // } /* sample pi2 */ mm=0.0; for(j=1;j<=nGlob;j++) mm+=psi2[j][ii]; pi2[ii]=betadev(nGlob-mm+1.0,mm+1.0,&setSeed); // if(i > numburn){ // filew1g << pi2[ii] << " "; // } /* sample pi3 */ mm=0.0; for(j=1;j<=nGlob;j++) mm+=psi3[j][ii]; pi3[ii]=betadev(nGlob-mm+1.0,mm+1.0,&setSeed); // if(i > numburn){ // filew1g << pi3[ii] << " "; // } /* sample Sigma1 given the rest */ /* note how I have set up the program to easily change assumptions regarding the structure of the variances-right now it has one variance per gene for each channel-so scaled across channels but not arrays can modify by commenting and uncommenting */ for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(x1Glob[j][k]-phi[1][ii]-theta1[j][ii] -psi1[j][ii]*(eta[j][ii]-theta1[j][ii])); for(k=1;k<=numRep;k++) vv+=DSQR(x2Glob[j][k]-phi[3][ii]-theta1[j][ii] -psi1[j][ii]*(eta[j][ii]-theta1[j][ii])); // sig1[j][ii]=sig2[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+2*numRep),2*numRep+nu0,&setSeed)); // vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(x3Glob[j][k]-phi[5][ii]-theta2[j][ii] -psi2[j][ii]*(eta[j][ii]-theta2[j][ii])); for(k=1;k<=numRep;k++) vv+=DSQR(x4Glob[j][k]-phi[7][ii]-theta2[j][ii] -psi2[j][ii]*(eta[j][ii]-theta2[j][ii])); // sig3[j][ii]=sig4[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+2*numRep),2*numRep+nu0,&setSeed)); // vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(x5Glob[j][k]-phi[9][ii]-theta3[j][ii] -psi3[j][ii]*(eta[j][ii]-theta3[j][ii])); for(k=1;k<=numRep;k++) vv+=DSQR(x6Glob[j][k]-phi[11][ii]-theta3[j][ii] -psi3[j][ii]*(eta[j][ii]-theta3[j][ii])); // sig5[j][ii]=sig6[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+2*numRep),2*numRep+nu0,&setSeed)); sig1[j][ii]=sig2[j][ii]=sig3[j][ii]=sig4[j][ii]=sig5[j][ii]=sig6[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+6*numRep),6*numRep+nu0,&setSeed)); // if(i > numburn){ // filew1g << sig1[j][ii] << " "; // } // sig1[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Sigma2 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(x2Glob[j][k]-phi[3][ii]-theta1[j][ii] -psi1[j][ii]*(eta[j][ii]-theta1[j][ii])); sig2[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Sigma3 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(x3Glob[j][k]-phi[5][ii]-theta2[j][ii] -psi2[j][ii]*(eta[j][ii]-theta2[j][ii])); sig3[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Sigma4 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(x4Glob[j][k]-phi[7][ii]-theta2[j][ii] -psi2[j][ii]*(eta[j][ii]-theta2[j][ii])); sig4[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Sigma5 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(x5Glob[j][k]-phi[9][ii]-theta3[j][ii] -psi3[j][ii]*(eta[j][ii]-theta3[j][ii])); sig5[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Sigma6 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(x6Glob[j][k]-phi[11][ii]-theta3[j][ii] -psi3[j][ii]*(eta[j][ii]-theta3[j][ii])); sig6[j][ii]=sqrt(invGamma((nu0*DSQR(sig0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Tau1 given the rest */ for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(y1Glob[j][k]-phi[2][ii]-eta[j][ii]); for(k=1;k<=numRep;k++) vv+=DSQR(y2Glob[j][k]-phi[4][ii]-eta[j][ii]); for(k=1;k<=numRep;k++) vv+=DSQR(y3Glob[j][k]-phi[6][ii]-eta[j][ii]); for(k=1;k<=numRep;k++) vv+=DSQR(y4Glob[j][k]-phi[8][ii]-eta[j][ii]); for(k=1;k<=numRep;k++) vv+=DSQR(y5Glob[j][k]-phi[10][ii]-eta[j][ii]); for(k=1;k<=numRep;k++) vv+=DSQR(y6Glob[j][k]-phi[12][ii]-eta[j][ii]); tau1[j][ii]=tau2[j][ii]=tau3[j][ii]=tau4[j][ii]=tau5[j][ii]=tau6[j][ii]=sqrt(invGamma((nu0*DSQR(tau0)+vv)/(nu0+6*numRep),6*numRep+nu0,&setSeed)); // if(i > numburn){ // filew1g << tau1[j][ii] << " "; // } // tau1[j][ii]=sqrt(invGamma((nu0*DSQR(tau0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Tau2 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(y2Glob[j][k]-phi[4][ii]-eta[j][ii]); tau2[j][ii]=sqrt(invGamma((nu0*DSQR(tau0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Tau3 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(y3Glob[j][k]-phi[6][ii]-eta[j][ii]); tau3[j][ii]=sqrt(invGamma((nu0*DSQR(tau0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Tau4 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(y4Glob[j][k]-phi[8][ii]-eta[j][ii]); tau4[j][ii]=sqrt(invGamma((nu0*DSQR(tau0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Tau5 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(y5Glob[j][k]-phi[10][ii]-eta[j][ii]); tau5[j][ii]=sqrt(invGamma((nu0*DSQR(tau0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } /* sample Tau6 given the rest */ /* for(j=1;j<=nGlob;j++){ vv=0.0; for(k=1;k<=numRep;k++) vv+=DSQR(y6Glob[j][k]-phi[12][ii]-eta[j][ii]); tau6[j][ii]=sqrt(invGamma((nu0*DSQR(tau0)+vv)/(nu0+numRep),numRep+nu0,&setSeed)); } */ if(i > numburn){ /* save fold changes */ for(j=1;j<=nGlob;j++){ fldchng[ii][i-numburn][j]=simres[ii][i-numburn][j]=theta1[j][ii]-theta2[j][ii]; fldchng[ii][i-numburn][j+nGlob]=simres[ii][i-numburn][j+nGlob]=theta1[j][ii]-theta3[j][ii]; fldchng[ii][i-numburn][j+2*nGlob]=simres[ii][i-numburn][j+2*nGlob]=theta2[j][ii]-theta3[j][ii]; } } if(i==400 && init==1){ /* can use this to write all parameters to file for some iterations */ for(j=1;j<=numpar;j++){ if(j<=nGlob) filew1g << theta1[j][ii] << " "; if((nGlob+1 <= j) && (j <= 2*nGlob)) filew1g << theta2[j-nGlob][ii] << " "; if((2*nGlob+1 <= j) && (j <= 3*nGlob)) filew1g << theta3[j-2*nGlob][ii] << " "; if((3*nGlob+1 <= j) & (j <= 4*nGlob)) filew1g << eta[j-3*nGlob][ii] << " "; if((4*nGlob+1 <= j) && (j <= 5*nGlob)) filew1g << psi1[j-4*nGlob][ii] << " "; if((5*nGlob+1 <= j) && (j <= 6*nGlob)) filew1g << psi2[j-5*nGlob][ii] << " "; if((6*nGlob+1 <= j) && (j <= 7*nGlob)) filew1g << psi3[j-6*nGlob][ii] << " "; if((7*nGlob+1 <= j) && (j <= 8*nGlob)) filew1g << sig1[j-7*nGlob][ii] << " "; if((8*nGlob+1 <= j) && (j <= 9*nGlob)) filew1g << sig2[j-8*nGlob][ii] << " "; if((9*nGlob+1 <= j) && (j <= 10*nGlob)) filew1g << sig3[j-9*nGlob][ii] << " "; if((10*nGlob+1 <= j) && (j <= 11*nGlob)) filew1g << sig4[j-10*nGlob][ii] << " "; if((11*nGlob+1 <= j) && (j <= 12*nGlob)) filew1g << sig5[j-11*nGlob][ii] << " "; if((12*nGlob+1 <= j) && (j <= 13*nGlob)) filew1g << sig6[j-12*nGlob][ii] << " "; if((13*nGlob+1 <= j) && (j <= 14*nGlob)) filew1g << tau1[j-13*nGlob][ii] << " "; if((14*nGlob+1 <= j) && (j <= 15*nGlob)) filew1g << tau2[j-14*nGlob][ii] << " "; if((15*nGlob+1 <= j) && (j <= 16*nGlob)) filew1g << tau3[j-15*nGlob][ii] << " "; if((16*nGlob+1 <= j) && (j <= 17*nGlob)) filew1g << tau4[j-16*nGlob][ii] << " "; if((17*nGlob+1 <= j) && (j <= 18*nGlob)) filew1g << tau5[j-17*nGlob][ii] << " "; if((18*nGlob+1 <= j) && (j <= 19*nGlob)) filew1g << tau6[j-18*nGlob][ii] << " "; if((19*nGlob+1 <= j) && (j <= 19*nGlob+12)) filew1g << phi[j-19*nGlob][ii] << " "; } filew1g << phiHier1[ii] << " "; filew1g << phiHier2[ii] << " "; filew1g << zeta[ii] << " "; filew1g << pi1[ii] << " "; filew1g << pi2[ii] << " "; filew1g << pi3[ii] << " "; } } if(i%100==0 && i > numburn && noChn > 1){ /* check convergence of gibbs sampler */ sqrtR(Rstat,noChn,i-numburn,3*nGlob,fldchng); cout << "max convergence diagnostic statistics after " << i << " iterations" << endl; maxRstat=0.0; for(j=1;j<=3*nGlob;j++){ if(Rstat[j]>maxRstat) maxRstat=Rstat[j]; } cout << maxRstat << endl; } } free_dvector(pi1,1,noChn); free_dvector(pi2,1,noChn); free_dvector(pi3,1,noChn); free_dvector(phiHier1,1,noChn); free_dvector(phiHier2,1,noChn); free_dvector(zeta,1,noChn); free_dmatrix(theta1,1,nGlob,1,noChn); free_dmatrix(theta2,1,nGlob,1,noChn); free_dmatrix(theta3,1,nGlob,1,noChn); free_dmatrix(eta,1,nGlob,1,noChn); free_dmatrix(psi1,1,nGlob,1,noChn); free_dmatrix(psi2,1,nGlob,1,noChn); free_dmatrix(psi3,1,nGlob,1,noChn); free_dmatrix(sig1,1,nGlob,1,noChn); free_dmatrix(sig2,1,nGlob,1,noChn); free_dmatrix(sig3,1,nGlob,1,noChn); free_dmatrix(sig4,1,nGlob,1,noChn); free_dmatrix(sig5,1,nGlob,1,noChn); free_dmatrix(sig6,1,nGlob,1,noChn); free_dmatrix(tau1,1,nGlob,1,noChn); free_dmatrix(tau2,1,nGlob,1,noChn); free_dmatrix(tau3,1,nGlob,1,noChn); free_dmatrix(tau4,1,nGlob,1,noChn); free_dmatrix(tau5,1,nGlob,1,noChn); free_dmatrix(tau6,1,nGlob,1,noChn); free_dmatrix(phi,1,12,1,noChn); free_dvector(Rstat,1,3*nGlob); free_d3tensor(fldchng,1,noChn,1,numIter-numburn,1,3*nGlob); } void sqrtR(double *stat,int noChn,int numiter,int numpar,double ***tens) { /* this calculates Gelman-Rubin sqrtR for assessing convergence for MCMC using overdispersed starting points-set trans=0 to calculate on original scale otherwise it calculates it on unbounded scale - if you want this yhou must modify below where the comments are */ double mean(double *,int); double var(double *,int); int i,j,k,trans=1; double *s,*m,*temp,W,inter,BonN; temp=dvector(1,numiter); s=dvector(1,noChn); m=dvector(1,noChn); for(i=1;i<=numpar;i++){ for(j=1;j<=noChn;j++){ for(k=1;k<=numiter;k++){ /* can calculate on the original rather than the unbounded scale by changing the next few lines */ if(trans==0 && i>noPar-1) temp[k]=exp(tens[j][k][i]); else temp[k]=tens[j][k][i]; } m[j]=mean(temp,numiter); s[j]=var(temp,numiter); } W=mean(s,noChn); BonN=var(m,noChn); inter=W*(numiter-1.0)/numiter+BonN; if(inter<0.0 || W<=0.0) stat[i]=100.0; else stat[i]=sqrt((inter)/W); } free_dvector(temp,1,numiter); free_dvector(s,1,noChn); free_dvector(m,1,noChn); } double mean(double *data, int n) { int j; double ave=0.0; for(j=1;j<=n;j++) ave += data[j]; ave /= n; return ave; } double var(double *data, int n) { double ave,dev,variance; double s=0.0; int j; ave=mean(data,n); for(j=1;j<=n;j++){ dev = data[j]-ave; s += dev*dev; } variance=s/(n-1.0); return variance; } /* matrix routines */ double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl, long ndh) { /* allocate a double 3-tensor */ long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1; double ***t; /* allocate pointers ot rows */ t=(double ***) malloc((size_t)((nrow+1)*sizeof(double**))); if(!t) nrerror("allocation failure 1 in d3tensor()"); t += 1; t -= nrl; /* allocate pointers to rows and set pointers to them */ t[nrl]=(double **) malloc((size_t)((nrow*ncol+1)*sizeof(double*))); if(!t[nrl]) nrerror("allocation failure 2 in d3tensor()"); t[nrl] += 1; t[nrl] -= ncl; /* allocate rows and set pointers to them */ t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+1)*sizeof(double))); if(!t[nrl][ncl]) nrerror("allocation failure 3 in d3tensor()"); t[nrl][ncl] += 1; t[nrl][ncl] -= ndl; for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep; for(i=nrl+1;i<=nrh;i++){ t[i]=t[i-1]+ncol; t[i][ncl]=t[i-1][ncl]+ncol*ndep; for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep; } /* return pointer to array of pointers to rows */ return t; } void free_d3tensor(double ***t,long nrl,long nrh,long ncl,long nch, long ndl,long ndh) { /* free a double tensor declared by d3tensor */ free((FREE_ARG) (t[nrl][ncl]+ndl-1)); free((FREE_ARG) (t[nrl]+ncl-1)); free((FREE_ARG) (t+nrl-1)); } double **dmatrix(long nrl,long nrh,long ncl,long nch) { /* allocate a double matrix with subscript range m[nrl...nrh][ncl...nch] */ long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; double **m; void nrerror(char error_text[]); /* allocate pointers to rows */ m=(double **) malloc((size_t)((nrow+1)*sizeof(double*))); if (!m) nrerror("allocation failure 1 in dmatrix()"); m += 1; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(double *) malloc((size_t)((nrow*ncol+1)*sizeof(double))); if (!m[nrl]) nrerror("allocation failure 2 in dmatrix()"); m[nrl] += 1; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } void free_dmatrix(double **m, long nrl,long nrh,long ncl,long nch) { /* free a double matrix allocated by dmatrix() */ free((FREE_ARG) (m[nrl]+ncl-1)); free((FREE_ARG) (m+nrl-1)); } double *dvector(long nl,long nh) { /* allocate a vector with index range v[nl...nh] */ void nrerror(char error_text[]); double *v; v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double))); if (!v) nrerror("allocation failure in dvector()"); return v-nl+1; } void free_dvector(double *v, long nl, long nh) { free((FREE_ARG) (v+nl-1)); } void nrerror(char error_text[]) { /* Numerical Recipes standard error handler */ fprintf(stderr,"Run-time error...\n"); fprintf(stderr,"%s\n",error_text); fprintf(stderr,"...now exiting to system...\n"); getch(); exit(1); } /* random number routines */ double invGamma(double scal,int df, long *setSeed) { /* generates an inverse gamma draw with df degrees of freedom and scale scal */ double x=0.0,result; for(int i=1;i<=df;i++) x += DSQR(gasdev(setSeed)); result=df*scal/x; return result; } double gasdev(long *idum) { /* gets a standard normal deviate */ double ran1(long *idum); static int iset=0; static double gset; double fac,rsq,v1,v2; if(iset==0){ do{ v1=2.0*ran1(idum)-1.0; v2=2.0*ran1(idum)-1.0; rsq=v1*v1+v2*v2; } while (rsq>=1.0 || rsq==0.0); fac=sqrt(-2.0*log(rsq)/rsq); gset=v1*fac; iset=1; return v2*fac; } else{ iset=0; return gset; } } #define IAran1 16807 #define IMran1 2147483647 #define AMran1 (1.0/IMran1) #define IQran1 127773 #define IRran1 2836 #define NDIVran1 (1+(IMran1-1)/32) #define EPSran1 1.2e-7 #define RNMXran1 (1.0-EPSran1) double ran1(long *idum) { /* gets a uniform deviate */ const int NTABran1=32; int j; long k; static long iy=0; static long iv[NTABran1]; double temp; if (*idum <= 0 || !iy) { if (-(*idum)<1) *idum=1; else *idum = -(*idum); for(j=NTABran1+7;j>=0;j--) { k=(*idum)/IQran1; *idum=IAran1*(*idum-k*IQran1)-IRran1*k; if (*idum<0) *idum += IMran1; if (jRNMXran1) return RNMXran1; else return temp; } double gamdev(double ia, long *idum) { /* gets a gamma deviate */ double am,e,s,v1,v2,x,y,p,q; if(ia<=0.0) nrerror("error in gamdev: parameter non-positive"); if(ia==1.0) x=-log(ran1(idum)); else{ if(ia<1.0){ p=exp(1.0)/(ia+exp(1.0)); do{ if(ran1(idum)=q); } else{ do{ do{ do{ v1=ran1(idum); v2=2.0*ran1(idum)-1.0; } while (v1*v1+v2*v2 > 1.0); y=v2/v1; am=ia-1.0; s=sqrt(2.0*am+1.0); x=s*y+am; } while(x <= 0); e=(1.0+y*y)*exp(am*log(x/am)-s*y); } while(ran1(idum) > e); } } return x; } double betadev(double a,double b, long *idum) { /* gets a beta deviate */ double xa,xb; xa=gamdev(a,idum); xb=gamdev(b,idum); return xa/(xa+xb); } /* function for reading in initial values */ void initlz(double **initValu,int noChn) { /* starting values for multiple chains-these initial values can be on the transformed scale (then set trans=1), e.g. log scale if they are var parameters */ int i,j; ifstream filer("c:\\Spluswin\\home\\initValu.dat"); if(!filer){ cerr << "Failure to open filer\n"; getch(); exit(EXIT_FAILURE); } for(i=1;i<=noChn;i++){ for(j=1;j<=noPar;j++){ filer >> initValu[j][i]; } } }