/************************************************************** * @HemantIshwaran * * Last updated: August 2001 * * ---------------------------- * * Some C-progams interfaced to various S-plus program * * * * Procedure (Splus 3.4 for Unix): * * --------- * * (1) compile in Unix: "cc -c bayesian.c" * * (2) load dynamically in S-PLUS: dyn.load("bayesian.o") * * * * Procedure (Splus 6.0) * * --------- * * For Splus 6.0 for Unix: * * * * (1) compile in Unix: "cc -c bayesian.c" * * (2) load dynamically in S-PLUS: dyn.open("bayesian.o") * * * * For Splus 6.0 for Linux: * * * * (0) Place bayesian.c in an Splus chapter you will use * * (1) Use the unix command "make" on this chapter * * (2) attach chapter * * * * ------------------------------------------------------ * * PLEASE GIVE CREDIT WHERE CREDIT IS DUE * * ------------------------------------------------------ * * Hemant Ishwaran ishwaran@bio.ri.ccf.org * * Dept of Biostatistics/Wb4 phone: 216-444-9932 * * Cleveland Clinic Foundation fax: 216-445-2781 * * 9500 Euclid Avenue * * Cleveland, OH 44195 * * * * www.bio.ri.ccf.org/Resume/Pages/Ishwaran * ***************************************************************/ /************************************************************** Assumptions: ----------- . Variable type int has the same limits as long int (true on most compilers; check /usr/include/sys/limits.h) /************************************************************** * Standard calls * ***************************************************************/ #include #include #include #include #include /*------------------------------------------------------------- . Some global declarations . Define constants for hyperparameters Be !!! CAREFUL !!! not to reuse these variable names later ---------------------------------------------------------------*/ #define FREE_ARG char* #define MAXLINE 10000 #define MCOL 30 #define NR_END 1 #define PI 3.141592654 #define SYSSMALL 1.0e-10 #define SYSLARGE 1.0e+99 static double dminarg1,dminarg2; #define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\ (dminarg1) : (dminarg2)) static double dmaxarg1,dmaxarg2; #define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\ (dmaxarg1) : (dmaxarg2)) static double dsqrarg; #define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg) static int imaxarg1,imaxarg2; #define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\ (imaxarg1) : (imaxarg2)) static int iminarg1,iminarg2; #define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\ (iminarg1) : (iminarg2)) static long lmaxarg1,lmaxarg2; #define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\ (lmaxarg1) : (lmaxarg2)) static long lminarg1,lminarg2; #define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\ (lminarg1) : (lminarg2)) #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) /* nrutil.h calls */ double **convert_dmatrix(double *a,long nrl,long nrh,long ncl,long nch); float **convert_matrix(float *a,long nrl,long nrh,long ncl,long nch); unsigned char *cvector(long nl,long nh); double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh); double **dmatrix(long nrl,long nrh,long ncl,long nch); double **dsubmatrix(double **a,long oldrl,long oldrh,long oldcl,long oldch, long newrl,long newcl); double *dvector(long nl,long nh); float ***f3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh); int **imatrix(long nrl,long nrh,long ncl,long nch); int *ivector(long nl,long nh); unsigned long *lvector(long nl,long nh); float **matrix(long nrl,long nrh,long ncl,long nch); void nrerror(char error_text[]); float **submatrix(float **a,long oldrl,long oldrh,long oldcl,long oldch, long newrl,long newcl); float *vector(long nl,long nh); void free_convert_dmatrix(double **b,long nrl,long nrh,long ncl,long nch); void free_convert_matrix(float **b,long nrl,long nrh,long ncl,long nch); void free_cvector(unsigned char *v,long nl,long nh); void free_d3tensor(double ***t,long nrl,long nrh,long ncl,long nch, long ndl,long ndh); void free_dmatrix(double **m,long nrl,long nrh,long ncl,long nch); void free_dsubmatrix(double **b,long nrl,long nrh,long ncl,long nch); void free_dvector(double *v,long nl,long nh); void free_f3tensor(float ***t,long nrl,long nrh,long ncl,long nch, long ndl,long ndh); void free_imatrix(int **m,long nrl,long nrh,long ncl,long nch); void free_ivector(int *v,long nl,long nh); void free_lvector(unsigned long *v,long nl,long nh); void free_matrix(float **m,long nrl,long nrh,long ncl,long nch); void free_submatrix(float **b,long nrl,long nrh,long ncl,long nch); void free_vector(float *v,long nl,long nh); /************************************************************** * Special Functions ***************************************************************/ /*------------------------------------------------------------- Prints the matrix a[1...r][1..c] in double precision fWmatrix --------------------------------------------------------------*/ void fWmatrix(FILE *fp,double **a,int r,int c) { int i,j; if (r==0 || c==0) return; for (i=1;i<=r;i++){ for (j=1;j<=c;j++){ fprintf(fp,"%15.8f ",a[i][j]); } fprintf(fp,"\n"); } } /*------------------------------------------------------------- Prints the vector a[1...d] in double precision fWvector --------------------------------------------------------------*/ void fWvector(FILE *fp,double *a,int d) { int i; if (d==0) return; for (i=1;i<=d;i++) fprintf(fp,"%15.8f ",a[i]); fprintf(fp,"\n"); } /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software Log of gamma function gammln --------------------------------------------------------------*/ double gammln(double xx) { double x,y,tmp,ser; static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5}; int j; y=x=xx; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.000000000190015; for (j=0;j<=5;j++) ser += cof[j]/++y; return -tmp+log(2.5066282746310005*ser/x); } /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software Series representation for incomplete gamma gser --------------------------------------------------------------*/ #define ITMAX 500 #define EPS 3.0e-7 void gser(double *gamser, double a, double x, double *gln) { double gammln(double xx); void nrerror(char error_text[]); int n; double sum,del,ap; *gln=gammln(a); if (x <= 0.0) { if (x < 0.0) nrerror("x less than 0 in routine gser"); *gamser=0.0; return; } else { ap=a; del=sum=1.0/a; for (n=1;n<=ITMAX;n++) { ++ap; del *= x/ap; sum += del; if (fabs(del) < fabs(sum)*EPS) { *gamser=sum*exp(-x+a*log(x)-(*gln)); return; } } return; } } #undef ITMAX #undef EPS /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software Continued fraction representation for incomplete gamma gcf --------------------------------------------------------------*/ #define ITMAX 500 #define EPS 3.0e-7 #define FPMIN 1.0e-30 void gcf(double *gammcf, double a, double x, double *gln) { double gammln(double xx); void nrerror(char error_text[]); int i; double an,b,c,d,del,h; *gln=gammln(a); b=x+1.0-a; c=1.0/FPMIN; d=1.0/b; h=d; for (i=1;i<=ITMAX;i++) { an = -i*(i-a); b += 2.0; d=an*d+b; if (fabs(d) < FPMIN) d=FPMIN; c=b+an/c; if (fabs(c) < FPMIN) c=FPMIN; d=1.0/d; del=d*c; h *= del; if (fabs(del-1.0) < EPS) break; } *gammcf=exp(-x+a*log(x)-(*gln))*h; } #undef ITMAX #undef EPS #undef FPMIN /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software Returns the incomplete gamma function gammp --------------------------------------------------------------*/ double gammp(double a, double x) { void gcf(double *gammcf, double a, double x, double *gln); void gser(double *gamser, double a, double x, double *gln); void nrerror(char error_text[]); double gamser,gammcf,gln; if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp"); if (x < (a+1.0)) { gser(&gamser,a,x,&gln); return gamser; } else { gcf(&gammcf,a,x,&gln); return 1.0-gammcf; } } /*------------------------------------------------------------- Given a vector arr[1..n] returns: (i) indx[1..n] which ranks arr[1..n] into ascending order (ii) rank[1..n] of the ranks (includes ties). (iii) unique values uniq[1..nclus] (iv) frequency of each unique value nuniq[1..nclus] Only the input quantities n and arr are unchanged. Modified from Numerical Recipes in C, page 338-339 sort --------------------------------------------------------------*/ #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp; #define M 7 void sort(int n,double *arr,int *indx, int *rank, double *uniq,int *nuniq,int *nclus) { int i,indxt,ir=n,temp,j,k,l=1; int jstack=0,*istack; double a; istack=ivector(1,n); for (j=1;j<=n;j++) indx[j]=j; for (;;) { if (ir-l < M) { for (j=l+1;j<=ir;j++) { indxt=indx[j]; a=arr[indxt]; for (i=j-1;i>=l;i--) { if (arr[indx[i]] <= a) break; indx[i+1]=indx[i]; } indx[i+1]=indxt; } if (jstack == 0) break; ir=istack[jstack--]; l=istack[jstack--]; } else { k=(l+ir) >> 1; SWAP(indx[k],indx[l+1]); if (arr[indx[l]] > arr[indx[ir]]) { SWAP(indx[l],indx[ir]) } if (arr[indx[l+1]] > arr[indx[ir]]) { SWAP(indx[l+1],indx[ir]) } if (arr[indx[l]] > arr[indx[l+1]]) { SWAP(indx[l],indx[l+1]) } i=l+1; j=ir; indxt=indx[l+1]; a=arr[indxt]; for (;;) { do i++; while (arr[indx[i]] < a); do j--; while (arr[indx[j]] > a); if (j < i) break; SWAP(indx[i],indx[j]) } indx[l+1]=indx[j]; indx[j]=indxt; jstack += 2; if (ir-i+1 >= j-l) { istack[jstack]=ir; istack[jstack-1]=i; ir=j-1; } else { istack[jstack]=j-1; istack[jstack-1]=l; l=i; } } } /* Compute unique values and their number of occurences */ for (j=1;j<=n;j++){ uniq[j]=0; nuniq[j]=0; } (*nclus)=1; rank[indx[*nclus]]=(*nclus); uniq[*nclus]=arr[indx[*nclus]]; for (j=2;j<=n;j++){ if(arr[indx[j]] > arr[indx[j-1]]){ (*nclus)++; uniq[*nclus]=arr[indx[j]]; } rank[indx[j]]=(*nclus); } for (j=1;j<=n;j++) nuniq[rank[j]]++; free_ivector(istack,1,n); } #undef M #undef SWAP /************************************************************** * Matrix operators ***************************************************************/ /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software Returns the Cholesky decomposition of a square matrix in lower diagonal form - matrix not altered. chold --------------------------------------------------------------*/ void chold(double **a,double **b,int dim) { double **temp,*diag,sum; int i,j,k; temp=dmatrix(1,dim,1,dim); diag=dvector(1,dim); for (i=1;i<=dim;i++){ for (j=i;j<=dim;j++) temp[i][j]=temp[j][i]=a[i][j]; } for (i=1;i<=dim;i++){ for (j=i;j<=dim;j++){ for (sum=temp[i][j],k=i-1;k>=1;k--) sum -= temp[i][k]*temp[j][k]; if (i == j) { if (sum <= 0.0) nrerror("choldc failed"); diag[i]=sqrt(sum); } else temp[j][i]=sum/diag[i]; } } for (i=1;i<=dim;i++){ b[i][i]=diag[i]; for (j=1;j<=i-1;j++){ b[i][j]=temp[i][j]; b[j][i]=0.; } } free_dmatrix(temp,1,dim,1,dim); free_dvector(diag,1,dim); } /*------------------------------------------------------------- convert vector to matrix createMatrix --------------------------------------------------------------*/ void createMatrix(double **newX,double *X,int nrow,int ncol) { int i,j; X=X-1; for (j=1;j<=ncol;j++){ for (i=1;i<=nrow;i++){ newX[i][j]=X[(j-1)*nrow+i]; } } } /*------------------------------------------------------------- convert integer vector to matrix createIntMatrix --------------------------------------------------------------*/ void createIntMatrix(int **newX,int *X,int nrow,int ncol) { int i,j; X=X-1; for (j=1;j<=ncol;j++){ for (i=1;i<=nrow;i++){ newX[i][j]=X[(j-1)*nrow+i]; } } } /*------------------------------------------------------------- convert matrix to vector createVector --------------------------------------------------------------*/ void createVector(double *X,double **newX,int nrow,int ncol) { int i,j; X=X-1; for (j=1;j<=ncol;j++){ for (i=1;i<=nrow;i++){ X[(j-1)*nrow+i]=newX[i][j]; } } } /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software lubskb --------------------------------------------------------------*/ void lubksb(double **a, int n, int *indx, double b[]) { int i,ii=0,ip,j; double sum; for (i=1;i<=n;i++) { ip=indx[i]; sum=b[ip]; b[ip]=b[i]; if (ii) for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; else if (sum) ii=i; b[i]=sum; } for (i=n;i>=1;i--) { sum=b[i]; for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; b[i]=sum/a[i][i]; } } /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software ludcmp --------------------------------------------------------------*/ #define TINY 1.0e-20; void ludcmp(double **a, int n, int *indx, double *d) { int i,imax,j,k; double big,dum,sum,temp; double *vv; vv=dvector(1,n); *d=1.0; for (i=1;i<=n;i++) { big=0.0; for (j=1;j<=n;j++) if ((temp=fabs(a[i][j])) > big) big=temp; if (big == 0.0) nrerror("Singular matrix in routine ludcmp"); vv[i]=1.0/big; } for (j=1;j<=n;j++) { for (i=1;i= big) { big=dum; imax=i; } } if (j != imax) { for (k=1;k<=n;k++) { dum=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=dum; } *d = -(*d); vv[imax]=vv[j]; } indx[j]=imax; if (a[j][j] == 0.0) a[j][j]=TINY; if (j != n) { dum=1.0/(a[j][j]); for (i=j+1;i<=n;i++) a[i][j] *= dum; } } free_dvector(vv,1,n); } #undef TINY /*------------------------------------------------------------- Multiply matrix to vector: c=ab matdot --------------------------------------------------------------*/ int matdot(double **a,double *b,double *c,int r1,int c1,int r2) { int i,k,sdim=c1; for (i=1;i<=r1;i++){ c[i]=0.; for (k=1;k<=sdim;k++){ c[i] += a[i][k] * b[k]; } } return (c1 != r2) ? -1 : 0; } /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software Inverse of matrix a[1..N][1..N] (which is not destroyed). Inverse stored in aInv matinv --------------------------------------------------------------*/ void matinv(double **a,double **aInv,int N) { int i,j,*indx; double **atemp,*col,d; indx=ivector(1,N); atemp=dmatrix(1,N,1,N); col=dvector(1,N); for (j=1;j<=N;j++){ for (i=1;i<=N;i++) atemp[i][j]=a[i][j]; } ludcmp(atemp,N,indx,&d); for (j=1;j<=N;j++){ for (i=1;i<=N;i++) col[i]=0.; col[j]=1.0; lubksb(atemp,N,indx,col); for (i=1;i<=N;i++) aInv[i][j]=col[i]; } free_ivector(indx,1,N); free_dmatrix(atemp,1,N,1,N); free_dvector(col,1,N); } /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software Inverse of matrix a[1..N][1..N] (which is not destroyed). Inverse stored in aInv Also returns determinant of inverse matinvDet --------------------------------------------------------------*/ void matinvDet(double **a,double **aInv,double *det,int N) { int i,j,*indx; double **atemp,*col,d; indx=ivector(1,N); atemp=dmatrix(1,N,1,N); col=dvector(1,N); for (j=1;j<=N;j++){ for (i=1;i<=N;i++) atemp[i][j]=a[i][j]; } ludcmp(atemp,N,indx,&d); (*det)=d; for (j=1;j<=N;j++){ (*det) *= atemp[j][j]; for (i=1;i<=N;i++) col[i]=0.; col[j]=1.0; lubksb(atemp,N,indx,col); for (i=1;i<=N;i++) aInv[i][j]=col[i]; } (*det)=1/(*det); free_ivector(indx,1,N); free_dmatrix(atemp,1,N,1,N); free_dvector(col,1,N); } /*------------------------------------------------------------- Dot product of two vectors: ab vecdot --------------------------------------------------------------*/ double vecdot(double *a,double *b,int d) { int i; double c=0.; for (i=1;i<=d;i++) c += a[i] * b[i]; return c; } /*------------------------------------------------------------- Add two vectors: c=a+b vecadd --------------------------------------------------------------*/ void vecadd(double *a,double *b,double *c,int d) { int i; for (i=1;i<=d;i++) c[i] = a[i] + b[i]; } /************************************************************** * Random Generators ***************************************************************/ /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software dran1 uniform random number generator ---------------------------------------------------------------*/ #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define NTAB 32 #define NDIV (1+(IM-1)/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) double dran1(int *idum) { int j; int k; static int iy=0; static int iv[NTAB]; double temp; if (*idum <= 0 || !iy) { if (-(*idum) < 1) *idum=1; else *idum = -(*idum); for (j=NTAB+7;j>=0;j--) { k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; if (*idum < 0) *idum += IM; if (j < NTAB) iv[j] = *idum; } iy=iv[0]; } k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; if (*idum < 0) *idum += IM; j=iy/NDIV; iy=iv[j]; iv[j] = *idum; if ((temp=AM*iy) > RNMX) return RNMX; else return temp; } #undef IA #undef IM #undef AM #undef IQ #undef IR #undef NTAB #undef NDIV #undef EPS #undef RNMX /*------------------------------------------------------------- Discrete random generator which uses a binary search. Especially useful when np>=30. Modified from Stochastic simulations, page 72 (original algorithm is flawed). rdisc --------------------------------------------------------------*/ double rdisc(double *supp,double *prob,int np,int *idum) { double u=dran1(idum),l=0.,r=np; double cum=0.,*cumprob=dvector(1,np); int i,j,step; for (j=1;j<=np;j++){ cum += prob[j]; /*Form cumulative probabilities*/ cumprob[j]=cum; } for (j=1;j<=np;j++) cumprob[j]=cumprob[j]/cum; for(;;){ i=floor((l+r)/2); if (u>cumprob[i]){ l=i; step=1; } else{ r=i; step=0; } if (l >= r-1){ free_dvector(cumprob,1,np); return supp[i+step]; } } } /*------------------------------------------------------------- Gamma random number generator from Stochastic Simulations pages 88-90 (split cases into alpha<1,alpha=1,alpha>1) rGamma --------------------------------------------------------------*/ double rGamma(double a,int *idum) { double e=exp(1.),b=e/(a+e); double temp,u0,u1,u2,c1=(a-1.),c2=((a-1./(6.*a))/c1); double c3=(2./c1),c4=(2.*a/c1),c5=(1./DSQR(a)),w; if (a == 1){ while((temp=dran1(idum)) != 0.) /* a = 1 */ return -log(temp); } else if (a < 1){ for(;;){ u0=dran1(idum);u1=dran1(idum); /* a < 1 */ if ( u0 > b){ temp=-1.*log((1-u0)/(a*b)); if (u1 <= pow(temp,a-1.)) return temp; } else{ temp=pow(u0/b,1./a); /* a > 1 */ if (u1 <= exp(-1.*temp)) return temp; } } } else{ for(;;){ do{ u1=dran1(idum);u2=dran1(idum); if (a > 2.5) u1=u2+c5*(1.-1.86*u1); } while (u1<0. || u1>1.); w=c2*u2/u1; if ( (c3*u1+w+1/w <= c4) || (c3*log(u1)-log(w)+w < 1) ) return c1*w; } } } /*------------------------------------------------------------- Truncated gamma(-0.5) on [x1,infty) Uses bisection method XMAX corresponds to infinity rGammat --------------------------------------------------------------*/ #define JMAX 150 #define XMAX 10000 #define SPI 1.772454 double rGammat(int *n3,double *x1,double *xacc,double *tgamma,int *seed) { void nrerror(char error_text[]); int i,j; double dx,f,fmid,xmid,rtb,bb,u; x1=x1-1;tgamma=tgamma-1; for (i=1;i<=*n3;i++){ bb=gammp(0.5,x1[i])+exp(-1.*x1[i])/(sqrt(x1[i])*SPI); u=dran1(seed); f=-1.*u; fmid=1.-u; rtb = f < 0.0 ? (dx=XMAX-x1[i],x1[i]) : (dx=x1[i]-XMAX,XMAX); for (j=1;j<=JMAX;j++){ xmid=rtb+(dx *= 0.5); fmid=(gammp(0.5,xmid)+exp(-1.*xmid)/(sqrt(xmid)*SPI)-bb)/(1-bb)-u; if (fmid <= 0.0) rtb=xmid; if (fabs(dx) < (*xacc) || fmid == 0.0) break; } tgamma[i]=rtb; } } #undef JMAX #undef XMAX #undef SPI /*------------------------------------------------------------- Simulate log gamma random variables using rejection sampling for shape parameter 0=(exp(exp(x*a[i])-exp(x))*bound)); plog[i]=x; } } /*------------------------------------------------------------- (C) Copr. 1986-92 Numerical Recipes Software Standard normal random number generator rnormal --------------------------------------------------------------*/ double rnormal(int *idum) { static int iset=0; static double gset; double fac,rsq,v1,v2; if (iset == 0) { do { v1=2.0*dran1(idum)-1.0; v2=2.0*dran1(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; } } /*------------------------------------------------------------- Generate a multivariate normal(mean,sigma) rmnormal --------------------------------------------------------------*/ void rmnormal(double *mean,double **sigma,double *nvec, int dim,int *idum) { double **var,*zz; int i,j; var=dmatrix(1,dim,1,dim); zz=dvector(1,dim); for (i=1;i<=dim;i++){ for (j=i;j<=dim;j++) var[i][j]=var[j][i]=sigma[i][j]; } chold(var,var,dim); for (i=1;i<=dim;i++) zz[i]=rnormal(idum); matdot(var,zz,nvec,dim,dim,dim); vecadd(mean,nvec,nvec,dim); free_dmatrix(var,1,dim,1,dim); free_dvector(zz,1,dim); } /************************************************************** * Bayesian Gibbs/Monte Carlo details ***************************************************************/ /*------------------------------------------------------------- Calculate areas Darea --------------------------------------------------------------*/ void Darea(double *p, double *z,double *yset, double *area,int *nbreaks,int *N) { int i,j; p=p-1;z=z-1;yset=yset-1;area=area-1; for (i=1;i<=((*nbreaks)-1);i++){ for (j=1;j<=*N;j++){ area[i] += p[j]*((yset[i]0){/* check if prior is dpn */ dpn1 = (*alpha)/(*maxtable); dpn2 = 1./(*maxtable); } if (strcmp(*conj,"T") == 0){/* conjugate case ? */ /* first customer details */ for (j=1;j<=*n;j++){ supp[j]=j; xsum[j]=xsum2[j]=gg[j]=0.; } kseq[1]=1; (*nKseq)=1; xsum[1]=x[1]; xsum2[1]=DSQR(x[1]); nc[1]=1; Var=xsum2[1]/(1+(*A)); c1=gammln(*a+0.5)-gammln(*a); c2=(1+*A)*2*(*b); (*lambda)=c1-0.5*log(PI*c2)-(*a+0.5)*log(1+xsum2[1]/c2); /* sequential seating assignment: customers 2...n */ for (i=2;i<=(*n);i++){ for (j=1;j<=(*nKseq);j++) gg[j]=DSQR(x[i])-DSQR(xsum[j]+x[i])*(*A)/(1+*A+nc[j]*(*A)) +DSQR(xsum[j])*(*A)/(1+nc[j]*(*A)); c1=*a+i/2; c2=(1+*A)*(2*(*b)+Var); c3=1+DSQR(x[i])/c2; Const=newtable=(*alpha)*DMAX((1-(*nKseq)*dpn2),0.0)/(sqrt(PI*c2)*pow(c3,c1)); for (j=1;j<=(*nKseq);j++){ c4=2*(*b)+Var; c5=1+nc[j]*(*A); c6=sqrt(c5/(PI*c4*(*A+c5))); c7=1+gg[j]/c4; oldtable[j]=(nc[j]+dpn1)*c6/pow(c7,c1); Const += oldtable[j]; } (*lambda) += log(Const); for (j=1;j<=(*nKseq);j++) prob[j]=oldtable[j]/Const; if ((dpn2==0) || (*nKseq)<(*maxtable)){ prob[(*nKseq)+1]=newtable/Const; rtable=rdisc(supp,prob,(*nKseq+1),seed); if (rtable==(*nKseq+1)) (*nKseq) += 1; } else{ rtable=rdisc(supp,prob,(*nKseq),seed); } kseq[i]=rtable; xsum[rtable] += x[i]; xsum2[rtable] += DSQR(x[i]); Var=0; nc[rtable] += 1; for (j=1;j<=(*nKseq);j++) Var += xsum2[j]-DSQR(xsum[j])/(1/(*A)+nc[j]); } (*theta)=(*b+Var/2)/rGamma((*a+(*n)/2),seed); /* last update for theta */ /* unique y values */ for (j=1;j<=(*nKseq);j++){ yuniq[j] = rnormal(seed)*sqrt(*A*Var/(*A*nc[j]+Var)) + (*A)/(*A*nc[j]+Var)*xsum[j]; } } else if (strcmp(*ngwcr,"F") == 0){/* non-conjugate case */ for (j=1;j<=*n;j++){ /* allows for theta mle */ supp[j]=j; xsum[j]=xsum2[j]=0.; xsum[j]=xsum2[j]=0.; } /* first customer details */ kseq[1]=1; (*nKseq)=1; xsum[1]=x[1]; xsum2[1]=DSQR(x[1]); nc[1]=1; c1=(*alpha)/sqrt(2*PI*(*theta+(*A))); c2=-1/(2*(*theta+(*A))); c3=2*PI*(*theta); c4=-1/(2*(*theta)); (*lambda)=c2*DSQR(x[1])-0.5*log((2*PI*(*theta+(*A)))); /* sequential seating assignment: customers 2...n */ for (i=2;i<=(*n);i++){ if ((i>=(*nmle)) && (strcmp(*thetamle,"T") == 0)){/* estimate theta? */ thetahat=degf=0; for (j=1;j<=(*nKseq);j++){ if (nc[j]>0){ degf += nc[j]-1; thetahat += (xsum2[j]-DSQR(xsum[j])/nc[j]); } } thetahat=thetahat/degf; if (degf>0 && thetahat>TINY) (*theta)=thetahat; c1=(*alpha)/sqrt(2*PI*(*theta+(*A))); c2=-1/(2*(*theta+(*A))); c3=2*PI*(*theta); c4=-1/(2*(*theta)); } Const=newtable=c1*DMAX((1-(*nKseq)*dpn2),0.0)*exp(c2*DSQR(x[i])); for (j=1;j<=(*nKseq);j++){ c5=(*theta)+(*A)*nc[j]; c6=c5+(*A); c7=sqrt(c5/(c3*c6)); oldtable[j]=(nc[j]+dpn1)*c7*exp(c4*(DSQR(x[i])-DSQR(xsum[j]+x[i])*(*A)/c6 +DSQR(xsum[j])*(*A)/c5)); Const += oldtable[j]; } (*lambda) += log(Const); for (j=1;j<=(*nKseq);j++) prob[j]=oldtable[j]/Const; if ((dpn2==0) || (*nKseq)<(*maxtable)){ prob[(*nKseq)+1]=newtable/Const; rtable=rdisc(supp,prob,(*nKseq+1),seed); if (rtable==(*nKseq+1)) (*nKseq) += 1; } else{ rtable=rdisc(supp,prob,(*nKseq),seed); } kseq[i]=rtable; xsum[rtable] += x[i]; xsum2[rtable] += DSQR(x[i]); nc[rtable] += 1; } /* last update for theta */ thetahat=degf=0; for (j=1;j<=(*nKseq);j++){ if (nc[j]>0){ degf += nc[j]-1; thetahat += (xsum2[j]-DSQR(xsum[j])/nc[j]); } } thetahat=thetahat/degf; if (degf>0 && thetahat>TINY) (*theta)=thetahat; /* unique y values */ for (j=1;j<=(*nKseq);j++){ yuniq[j] = rnormal(seed)*sqrt(*A*(*theta)/(*A*nc[j]+(*theta))) + (*A)/(*A*nc[j]+(*theta))*xsum[j]; } } else{/* ngwcr case */ c3=2*PI*(*theta); c4=-1/(2*(*theta)); /* first customer */ for (j=1;j<=(*n);j++){ supp[j]=j; prob[j]=1; xsum[j]=xsum2[j]=0.; } denom[1]=0; for (k=1;k<=(*N);k++){ ix=rdisc(supp,prob,(*n),seed); /* bootstrap reference H_N */ suppz[k]=z[k]=Zngwcr[ix]; ksum[k][1]=DSQR(x[1]-z[k]); denom[1] += exp(c4*ksum[k][1]); } kseq[1]=1; (*nKseq)=1; xsum[1]=x[1]; xsum2[1]=DSQR(x[1]); nc[1]=1; /* sequential seating assignment: customers 2...n */ for (i=2;i<=(*n);i++){ if ((i>=(*nmle)) && (strcmp(*thetamle,"T") == 0)){/* estimate theta? */ thetahat=degf=0; for (j=1;j<=(*nKseq);j++){ if (nc[j]>0){ degf += nc[j]-1; thetahat += (xsum2[j]-DSQR(xsum[j])/nc[j]); } } thetahat=thetahat/degf; if (degf>0 && thetahat>TINY) (*theta)=thetahat; c3=2*PI*(*theta); c4=-1/(2*(*theta)); } newtable=0; for (k=1;k<=(*N);k++){ newtable += exp(c4*DSQR(x[i]-z[k])); } Const=newtable/((*N)*sqrt(c3)); newtable=Const*(*alpha)*DMAX((1-(*nKseq)*dpn2),0.0); for (j=1;j<=(*nKseq);j++){ oldtable[j]=0; for (k=1;k<=(*N);k++) oldtable[j] += exp(c4*(ksum[k][j]+DSQR(x[i]-z[k]))); oldtable[j] = (nc[j]+dpn1)*oldtable[j]/(denom[j]*sqrt(c3)); Const += oldtable[j]; } *lambda += log(Const); for (j=1;j<=(*nKseq);j++) prob[j]=oldtable[j]/Const; if ((dpn2==0) || (*nKseq)<(*maxtable)){ prob[(*nKseq)+1]=newtable/Const; rtable=rdisc(supp,prob,(*nKseq+1),seed); denom[rtable]=0; if (rtable==(*nKseq+1)){ (*nKseq) += 1; for (k=1;k<=(*N);k++){ ksum[k][rtable]=DSQR(x[i]-z[k]); denom[rtable] += exp(c4*ksum[k][rtable]); } } else{ for (k=1;k<=(*N);k++){ ksum[k][rtable] += DSQR(x[i]-z[k]); denom[rtable] += exp(c4*ksum[k][rtable]); } } } else{ rtable=rdisc(supp,prob,(*nKseq),seed); denom[rtable]=0; for (k=1;k<=(*N);k++){ ksum[k][rtable] += DSQR(x[i]-z[k]); denom[rtable] += exp(c4*ksum[k][rtable]); } } kseq[i]=rtable; xsum[rtable] += x[i]; xsum2[rtable] += DSQR(x[i]); nc[rtable] += 1; } /* last update for theta */ thetahat=degf=0; for (j=1;j<=(*nKseq);j++){ if (nc[j]>0){ degf += nc[j]-1; thetahat += (xsum2[j]-DSQR(xsum[j])/nc[j]); } } thetahat=thetahat/degf; if (degf>0 && thetahat>TINY) (*theta)=thetahat; /* unique y values */ for (j=1;j<=(*nKseq);j++){ for (k=1;k<=(*N);k++) probz[k]=exp(c4*ksum[k][j])/(denom[j]*sqrt(c3)); yuniq[j]=rdisc(suppz,probz,(*N),seed); } } free_dmatrix(ksum,1,*N,1,*n); free_dvector(xsum,1,*n); free_dvector(xsum2,1,*n); free_dvector(gg,1,*n); free_dvector(oldtable,1,*n); free_dvector(supp,1,*n); free_dvector(prob,1,*n); free_dvector(z,1,*N); free_dvector(suppz,1,*N); free_dvector(probz,1,*N); free_dvector(denom,1,*n); } /*------------------------------------------------------------- GWCR algorithm for linear mixed models Fits general random effects gwcrlmm --------------------------------------------------------------*/ void gwcrlmm(char **sigmamle, double *Y,double *Z,double *DPalpha,double *sigma, double *alphaUniq,double *alphaMean,double *alphaVar, int *kseq,int *nc,double *delta,double *lambda, int *nKseq,int *nmle,double *A, int *n,int *s,int *maxtable,int *seed) { double ***ZZ; double **YZ,**newZ,**varnew,**varnewInv,**var1,**var2; double **newAlpha,**newAlphamean,**newAlphavar; double *YZnew,*mean1,*mean2,*crossProd,*aa,*Ysum2; double *oldtable,*supp,*prob; double c1,c2,c3,c4,detnew,det1,det2; double Const,newtable,sigmahat,degf; double dpn1=0,dpn2=0,TINY=1e-10,q1,q2; int i,j,k,l,rtable; Y=Y-1;kseq=kseq-1;nc=nc-1;delta=delta-1; ZZ=d3tensor(1,*s,1,*s,1,*n); YZ=dmatrix(1,*s,1,*n); newZ=dmatrix(1,*n,1,*s); varnew=dmatrix(1,*s,1,*s); varnewInv=dmatrix(1,*s,1,*s); var1=dmatrix(1,*s,1,*s); var2=dmatrix(1,*s,1,*s); newAlpha=dmatrix(1,*n,1,*s); newAlphamean=dmatrix(1,*n,1,*s); newAlphavar=dmatrix(1,*n,1,*s); YZnew=dvector(1,*s); mean1=dvector(1,*s); mean2=dvector(1,*s); crossProd=dvector(1,*s); aa=dvector(1,*s); Ysum2=dvector(1,*n); oldtable=dvector(1,*n); supp=dvector(1,*n); prob=dvector(1,*n); createMatrix(newZ,Z,*n,*s); if ((*maxtable)>0){/* check if prior is dpn */ dpn1 = (*DPalpha)/(*maxtable); dpn2 = 1./(*maxtable); } /* START GWCR */ for (i=1;i<=(*n);i++){ supp[i]=i; nc[i]=0; Ysum2[i]=0.; for (k=1;k<=(*s);k++){ YZ[k][i]=0; for (l=k;l<=(*s);l++) ZZ[k][l][i]=0; } } /* first customer details */ kseq[1]=1; (*nKseq)=1; nc[1]=1; Ysum2[1]=DSQR(Y[1]); for (k=1;k<=(*s);k++){ YZnew[k]=YZ[k][1]=Y[1]*newZ[1][k]; for (l=k;l<=(*s);l++){ ZZ[k][l][1]=newZ[1][k]*newZ[1][l]; if (l==k) varnew[k][l]=(*sigma)/(*A)+ZZ[k][l][1]; else varnew[k][l]=varnew[l][k]=ZZ[k][l][1]; } } matinvDet(varnew,varnew,&detnew,*s); matdot(varnew,YZnew,crossProd,*s,*s,*s); c1=pow(*sigma,(*s-1)/2)/sqrt(2*PI*pow(*A,*s)); c2=-1/(2*(*sigma)); c3=2*PI*(*sigma); (*lambda)=log(c1)+0.5*log(detnew)+ c2*((DSQR(Y[1])-vecdot(YZnew,crossProd,*s))); /* sequential seating assignment: customers 2...n */ for (i=2;i<=(*n);i++){ if ((i>=(*nmle)) && (strcmp(*sigmamle,"T") == 0)){/* estimate sigma? */ sigmahat=degf=0; for (j=1;j<=(*nKseq);j++){ if (nc[j]>(*s)){ for (k=1;k<=(*s);k++){ mean1[k]=YZ[k][j]; for (l=k;l<=(*s);l++){ if (l==k) var1[k][l]=delta[k]+ZZ[k][l][j]; else var1[k][l]=var1[l][k]=ZZ[k][l][j]; } } matinv(var1,var1,*s); matdot(var1,mean1,crossProd,*s,*s,*s); sigmahat += Ysum2[j]-vecdot(mean1,crossProd,*s); degf += nc[j]-*s; } } sigmahat = sigmahat/degf; if (degf>0 && sigmahat>TINY) (*sigma)=sigmahat; c1=pow(*sigma,(*s-1)/2)/sqrt(2*PI*pow(*A,*s)); c2=-1/(2*(*sigma)); c3=2*PI*(*sigma); } for (k=1;k<=(*s);k++){ YZnew[k]=Y[i]*newZ[i][k]; for (l=k;l<=(*s);l++){ if (l==k) varnew[k][l]=(*sigma)/(*A)+newZ[i][k]*newZ[i][l]; else varnew[k][l]=varnew[l][k]=newZ[i][k]*newZ[i][l]; } } matinvDet(varnew,varnewInv,&detnew,*s); matdot(varnewInv,YZnew,crossProd,*s,*s,*s); Const=newtable=(*DPalpha)*DMAX((1-(*nKseq)*dpn2),0.0)*c1*sqrt(detnew)*( exp(c2*(DSQR(Y[i])-vecdot(YZnew,crossProd,*s)))); for (j=1;j<=(*nKseq);j++){ for (k=1;k<=(*s);k++){ mean1[k]=YZ[k][j]; mean2[k]=mean1[k]+YZnew[k]; for (l=k;l<=(*s);l++){ if (l==k) var1[k][l]=(*sigma)/(*A)+ZZ[k][l][j]; else var1[k][l]=var1[l][k]=ZZ[k][l][j]; var2[k][l]=var2[l][k]=ZZ[k][l][j]+varnew[k][l]; } } matinvDet(var1,var1,&det1,*s); matinvDet(var2,var2,&det2,*s); matdot(var1,mean1,crossProd,*s,*s,*s); q1=vecdot(mean1,crossProd,*s); matdot(var2,mean2,crossProd,*s,*s,*s); q2=vecdot(mean2,crossProd,*s); oldtable[j]=(nc[j]+dpn1)*sqrt(det2/(c3*det1))*( exp(c2*(DSQR(Y[i])-q2+q1))); Const += oldtable[j]; } (*lambda) += log(Const); for (j=1;j<=(*nKseq);j++) prob[j]=oldtable[j]/Const; if ((dpn2==0) || (*nKseq)<(*maxtable)){ prob[(*nKseq)+1]=newtable/Const; rtable=rdisc(supp,prob,(*nKseq+1),seed); if (rtable==(*nKseq+1)) (*nKseq) += 1; } else{ rtable=rdisc(supp,prob,(*nKseq),seed); } kseq[i]=rtable; Ysum2[rtable] += DSQR(Y[i]); for (k=1;k<=(*s);k++){ YZ[k][rtable] += Y[i]*newZ[i][k]; for (l=k;l<=(*s);l++) ZZ[k][l][rtable] += newZ[i][k]*newZ[i][l]; } nc[rtable] += 1; } /* FINISHED GWCR */ /* now compute unique alpha values */ /* and last value for sigma */ /* careful (!) with small tables */ if (strcmp(*sigmamle,"T") == 0){ sigmahat=degf=0; for (j=1;j<=(*nKseq);j++){ if (nc[j]>(*s)){ for (k=1;k<=(*s);k++){ mean1[k]=YZ[k][j]; for (l=k;l<=(*s);l++){ if (l==k) var1[k][l]=delta[k]+ZZ[k][l][j]; else var1[k][l]=var1[l][k]=ZZ[k][l][j]; } } matinv(var1,var1,*s); matdot(var1,mean1,crossProd,*s,*s,*s); sigmahat += Ysum2[j]-vecdot(mean1,crossProd,*s); degf += nc[j]-*s; } } sigmahat = sigmahat/degf; if (degf>0 && sigmahat>TINY) (*sigma)=sigmahat; } for (j=1;j<=(*nKseq);j++){ for (k=1;k<=(*s);k++){ mean2[k]=YZ[k][j]/(*sigma); for (l=k;l<=(*s);l++){ if (l==k) var2[k][l]=delta[k]+ZZ[k][l][j]/(*sigma); else var2[k][l]=var2[l][k]=ZZ[k][l][j]/(*sigma); } } matinv(var2,var2,*s); matdot(var2,mean2,crossProd,*s,*s,*s); if (nc[j]<(*s)){ for (k=1;k<=(*s);k++){ newAlpha[j][k]=crossProd[k]; newAlphamean[j][k]=crossProd[k]; newAlphavar[j][k]=0; } } else{ rmnormal(crossProd,var2,aa,*s,seed); for (k=1;k<=(*s);k++){ newAlpha[j][k]=aa[k]; newAlphamean[j][k]=crossProd[k]; newAlphavar[j][k]=var2[k][k]; } } } createVector(alphaUniq,newAlpha,*n,*s); createVector(alphaMean,newAlphamean,*n,*s); createVector(alphaVar,newAlphavar,*n,*s); free_d3tensor(ZZ,1,*s,1,*s,1,*n); free_dmatrix(YZ,1,*s,1,*n); free_dmatrix(newZ,1,*n,1,*s); free_dmatrix(varnew,1,*s,1,*s); free_dmatrix(varnewInv,1,*s,1,*s); free_dmatrix(var1,1,*s,1,*s); free_dmatrix(var2,1,*s,1,*s); free_dmatrix(newAlpha,1,*n,1,*s); free_dmatrix(newAlphamean,1,*n,1,*s); free_dmatrix(newAlphavar,1,*n,1,*s); free_dvector(YZnew,1,*s); free_dvector(mean1,1,*s); free_dvector(mean2,1,*s); free_dvector(crossProd,1,*s); free_dvector(aa,1,*s); free_dvector(Ysum2,1,*n); free_dvector(oldtable,1,*n); free_dvector(supp,1,*n); free_dvector(prob,1,*n); } /*------------------------------------------------------------- GWCR algorithm for longitudinal linear mixed models Fits general random effects gwcrlmmlong [sigma update still under construction] --------------------------------------------------------------*/ void gwcrlmmlong(char **sigmamle, int *indx,int *nrec,int *maxrec,int *ntotal, double *Y,double *Z,double *DPalpha,double *sigma, double *alphaUniq,double *alphaMean,double *alphaVar, double *meanprior,double *varprior, int *kseq,int *nc,double *delta,double *lambda, int *nKseq,int *nmle, int *n,int *s,int *maxtable,int *seed) { double ***ZZ; double **YZ,**newZ,**varnew,**varnewInv,**var1,**var2; double **newAlpha,**newAlphamean,**newAlphavar; double *YZnew,*mean1,*mean2,*crossProd,*aa,*Ysum2; double *oldtable,*supp,*prob,*crossprior; double c1,c2,c3,detnew,det1,det2; double Const,newtable,sigmahat,degf,Ycross; double dpn1=0,dpn2=0,TINY=1e-10,q1,q2,qprior=0,dprior=1; int **Indx; int i,j,k,l,m,rtable; nrec=nrec-1;Y=Y-1;kseq=kseq-1;nc=nc-1;delta=delta-1; meanprior=meanprior-1;varprior=varprior-1; ZZ=d3tensor(1,*s,1,*s,1,*n); Indx=imatrix(1,*n,1,*maxrec); YZ=dmatrix(1,*s,1,*n); newZ=dmatrix(1,*ntotal,1,*s); varnew=dmatrix(1,*s,1,*s); varnewInv=dmatrix(1,*s,1,*s); var1=dmatrix(1,*s,1,*s); var2=dmatrix(1,*s,1,*s); newAlpha=dmatrix(1,*n,1,*s); newAlphamean=dmatrix(1,*n,1,*s); newAlphavar=dmatrix(1,*n,1,*s); crossprior=dvector(1,*s); YZnew=dvector(1,*s); mean1=dvector(1,*s); mean2=dvector(1,*s); crossProd=dvector(1,*s); aa=dvector(1,*s); Ysum2=dvector(1,*n); oldtable=dvector(1,*n); supp=dvector(1,*n); prob=dvector(1,*n); createIntMatrix(Indx,indx,*n,*maxrec); createMatrix(newZ,Z,*ntotal,*s); if ((*maxtable)>0){/* check if prior is dpn */ dpn1 = (*DPalpha)/(*maxtable); dpn2 = 1./(*maxtable); } /* START GWCR */ for (i=1;i<=(*n);i++){ supp[i]=i; nc[i]=0; Ysum2[i]=0; for (k=1;k<=(*s);k++){ YZ[k][i]=0; for (l=k;l<=(*s);l++) ZZ[k][l][i]=0; } } /* first customer details */ kseq[1]=1; (*nKseq)=1; nc[1]=1; for (m=1;m<=nrec[1];m++) Ysum2[1] += DSQR(Y[Indx[1][m]]); for (k=1;k<=(*s);k++){ crossprior[k]=meanprior[k]/varprior[k]; qprior += meanprior[k]*crossprior[k]; dprior *= varprior[k]; for (m=1;m<=nrec[1];m++) YZ[k][1] += Y[Indx[1][m]]*newZ[Indx[1][m]][k]; YZnew[k]=(*sigma)*crossprior[k]+YZ[k][1]; for (l=k;l<=(*s);l++){ for (m=1;m<=nrec[1];m++) ZZ[k][l][1] += newZ[Indx[1][m]][k]*newZ[Indx[1][m]][l]; if (l==k) varnew[k][l]=(*sigma)/varprior[k]+ZZ[k][l][1]; else varnew[k][l]=varnew[l][k]=ZZ[k][l][1]; } } matinvDet(varnew,varnew,&detnew,*s); matdot(varnew,YZnew,crossProd,*s,*s,*s); c1=sqrt(pow(*sigma,*s-nrec[1])/(pow(2*PI,nrec[1])*dprior)); c2=-1/(2*(*sigma)); qprior=qprior/2; (*lambda)=log(c1)+0.5*log(detnew)+ c2*(Ysum2[1]-vecdot(YZnew,crossProd,*s))-qprior; /* sequential seating assignment: customers 2...n */ for (i=2;i<=(*n);i++){ if ((i>=(*nmle)) && (strcmp(*sigmamle,"T") == 0)){/* estimate sigma? */ sigmahat=degf=0; for (j=1;j<=(*nKseq);j++){ if (nc[j]>(*s)){ for (k=1;k<=(*s);k++){ mean1[k]=YZ[k][j]; for (l=k;l<=(*s);l++){ if (l==k) var1[k][l]=delta[k]+ZZ[k][l][j]; else var1[k][l]=var1[l][k]=ZZ[k][l][j]; } } matinv(var1,var1,*s); matdot(var1,mean1,crossProd,*s,*s,*s); sigmahat += Ysum2[j]-vecdot(mean1,crossProd,*s); degf += nc[j]-*s; } } sigmahat = sigmahat/degf; c2=-1/(2*(*sigma)); } c1=sqrt(pow(*sigma,*s-nrec[i])/(pow(2*PI,nrec[i])*dprior)); c3=pow(2*PI*(*sigma),nrec[i]); Ycross=0; for (m=1;m<=nrec[i];m++) Ycross += DSQR(Y[Indx[i][m]]); for (k=1;k<=(*s);k++){ YZnew[k]=(*sigma)*crossprior[k]; for (m=1;m<=nrec[i];m++) YZnew[k] += Y[Indx[i][m]]*newZ[Indx[i][m]][k]; for (l=k;l<=(*s);l++){ if (l==k){ varnew[k][l]=(*sigma)/varprior[k]; for (m=1;m<=nrec[i];m++) varnew[k][l] += newZ[Indx[i][m]][k]*newZ[Indx[i][m]][l]; } else{ varnew[k][l]=0; for (m=1;m<=nrec[i];m++) varnew[k][l] += newZ[Indx[i][m]][k]*newZ[Indx[i][m]][l]; varnew[l][k]=varnew[k][l]; } } } matinvDet(varnew,varnewInv,&detnew,*s); matdot(varnewInv,YZnew,crossProd,*s,*s,*s); Const=newtable=(*DPalpha)*DMAX((1-(*nKseq)*dpn2),0.0)*c1*sqrt(detnew)*( exp(c2*(Ycross-vecdot(YZnew,crossProd,*s))-qprior)); for (j=1;j<=(*nKseq);j++){ for (k=1;k<=(*s);k++){ mean1[k]=(*sigma)*crossprior[k]+YZ[k][j]; mean2[k]=YZnew[k]+YZ[k][j]; for (l=k;l<=(*s);l++){ if (l==k) var1[k][l]=(*sigma)/varprior[k]+ZZ[k][l][j]; else var1[k][l]=var1[l][k]=ZZ[k][l][j]; var2[k][l]=var2[l][k]=ZZ[k][l][j]+varnew[k][l]; } } matinvDet(var1,var1,&det1,*s); matinvDet(var2,var2,&det2,*s); matdot(var1,mean1,crossProd,*s,*s,*s); q1=vecdot(mean1,crossProd,*s); matdot(var2,mean2,crossProd,*s,*s,*s); q2=vecdot(mean2,crossProd,*s); oldtable[j]=(nc[j]+dpn1)*sqrt(det2/(c3*det1))*( exp(c2*(Ycross-q2+q1))); Const += oldtable[j]; } (*lambda) += log(Const); for (j=1;j<=(*nKseq);j++) prob[j]=oldtable[j]/Const; if ((dpn2==0) || (*nKseq)<(*maxtable)){ prob[(*nKseq)+1]=newtable/Const; rtable=rdisc(supp,prob,(*nKseq+1),seed); if (rtable==(*nKseq+1)) (*nKseq) += 1; } else{ rtable=rdisc(supp,prob,(*nKseq),seed); } kseq[i]=rtable; for (m=1;m<=nrec[i];m++) Ysum2[rtable] += DSQR(Y[Indx[i][m]]); for (k=1;k<=(*s);k++){ for (m=1;m<=nrec[i];m++) YZ[k][rtable] += Y[Indx[i][m]]*newZ[Indx[i][m]][k]; for (l=k;l<=(*s);l++){ for (m=1;m<=nrec[i];m++) ZZ[k][l][rtable] += newZ[Indx[i][m]][k]*newZ[Indx[i][m]][l]; } } nc[rtable] += 1; } /* FINISHED GWCR */ /* now compute unique alpha values */ /* and last value for sigma */ /* careful (!) with small tables */ if (strcmp(*sigmamle,"T") == 0){ sigmahat=degf=0; for (j=1;j<=(*nKseq);j++){ if (nc[j]>(*s)){ for (k=1;k<=(*s);k++){ mean1[k]=YZ[k][j]; for (l=k;l<=(*s);l++){ if (l==k) var1[k][l]=delta[k]+ZZ[k][l][j]; else var1[k][l]=var1[l][k]=ZZ[k][l][j]; } } matinv(var1,var1,*s); matdot(var1,mean1,crossProd,*s,*s,*s); sigmahat += Ysum2[j]-vecdot(mean1,crossProd,*s); degf += nc[j]-*s; } } sigmahat = sigmahat/degf; if (degf>0 && sigmahat>TINY) (*sigma)=sigmahat; } for (j=1;j<=(*nKseq);j++){ for (k=1;k<=(*s);k++){ mean2[k]=YZ[k][j]/(*sigma)+crossprior[k]; for (l=k;l<=(*s);l++){ if (l==k) var2[k][l]=1/varprior[k]+ZZ[k][l][j]/(*sigma); else var2[k][l]=var2[l][k]=ZZ[k][l][j]/(*sigma); } } matinv(var2,var2,*s); matdot(var2,mean2,crossProd,*s,*s,*s); rmnormal(crossProd,var2,aa,*s,seed); for (k=1;k<=(*s);k++){ newAlpha[j][k]=aa[k]; newAlphamean[j][k]=crossProd[k]; newAlphavar[j][k]=var2[k][k]; } } createVector(alphaUniq,newAlpha,*n,*s); createVector(alphaMean,newAlphamean,*n,*s); createVector(alphaVar,newAlphavar,*n,*s); free_d3tensor(ZZ,1,*s,1,*s,1,*n); free_imatrix(Indx,1,*n,1,*maxrec); free_dmatrix(YZ,1,*s,1,*n); free_dmatrix(newZ,1,*ntotal,1,*s); free_dmatrix(varnew,1,*s,1,*s); free_dmatrix(varnewInv,1,*s,1,*s); free_dmatrix(var1,1,*s,1,*s); free_dmatrix(var2,1,*s,1,*s); free_dmatrix(newAlpha,1,*n,1,*s); free_dmatrix(newAlphamean,1,*n,1,*s); free_dmatrix(newAlphavar,1,*n,1,*s); free_dvector(crossprior,1,*s); free_dvector(YZnew,1,*s); free_dvector(mean1,1,*s); free_dvector(mean2,1,*s); free_dvector(crossProd,1,*s); free_dvector(aa,1,*s); free_dvector(Ysum2,1,*n); free_dvector(oldtable,1,*n); free_dvector(supp,1,*n); free_dvector(prob,1,*n); } /*------------------------------------------------------------- N-GWCR process for poisson process. Can be extended to GWCR case. gwcrPoisson --------------------------------------------------------------*/ void gwcrPoiss(double *x,double *y,double *Zx,double *Zy, double *xUnique, double *yUnique,double *lambda, double *theta1,double *theta2, int *nc,int *ntables,int *n,int *N,int *seed) { double **ksum; double *oldtable,*supp,*prob,*denom,*suppz,*probz,newtable,Const; double c1=2*PI*(*theta2),c2=-1/(2*(*theta2)); double c3=(exp(*theta1)/(1+exp(*theta1))); int i,j,k,rtable; x=x-1;y=y-1;Zx=Zx-1;Zy=Zy-1;xUnique=xUnique-1;yUnique=yUnique-1; nc=nc-1; ksum=dmatrix(1,*N,1,*n); oldtable=dvector(1,*n); supp=dvector(1,*n); prob=dvector(1,*n); denom=dvector(1,*n); suppz=dvector(1,*N); probz=dvector(1,*N); /* first customer */ for (j=1;j<=(*n);j++) supp[j]=j; denom[1]=0; for (k=1;k<=(*N);k++){ suppz[k]=k; ksum[k][1]=DSQR(x[1]-Zx[k])+DSQR(y[1]-Zy[k]); denom[1] += exp(c2*ksum[k][1]); } (*ntables)=1; nc[1]=1; /* sequential seating assignment: customers 2...n */ for (i=2;i<=(*n);i++){ newtable=0; for (k=1;k<=(*N);k++){ newtable += exp(c2*(DSQR(x[i]-Zx[k])+DSQR(y[i]-Zy[k]))); } newtable=c3*newtable/((*N)*c1); Const=newtable; for (j=1;j<=(*ntables);j++){ oldtable[j]=0; for (k=1;k<=(*N);k++) oldtable[j] += exp(c2*(ksum[k][j]+DSQR(x[i]-Zx[k])+DSQR(y[i]-Zy[k]))); oldtable[j] = c3*nc[j]*oldtable[j]/(denom[j]*c1); Const += oldtable[j]; } *lambda += log(Const); for (j=1;j<=(*ntables);j++) prob[j]=oldtable[j]/Const; prob[(*ntables)+1]=newtable/Const; rtable=rdisc(supp,prob,(*ntables+1),seed); denom[rtable]=0; if (rtable==(*ntables+1)){ (*ntables) += 1; for (k=1;k<=(*N);k++){ ksum[k][rtable]=DSQR(x[i]-Zx[k])+DSQR(y[i]-Zy[k]); denom[rtable] += exp(c2*ksum[k][rtable]); } } else{ for (k=1;k<=(*N);k++){ ksum[k][rtable] += DSQR(x[i]-Zx[k])+DSQR(y[i]-Zy[k]); denom[rtable] += exp(c2*ksum[k][rtable]); } } nc[rtable] += 1; } /* unique values */ for (j=1;j<=(*ntables);j++){ for (k=1;k<=(*N);k++) probz[k]=exp(c2*ksum[k][j])/(denom[j]*c3); rtable=rdisc(suppz,probz,(*N),seed); xUnique[j]=Zx[rtable]; yUnique[j]=Zy[rtable]; } free_dmatrix(ksum,1,*N,1,*n); free_dvector(oldtable,1,*n); free_dvector(supp,1,*n); free_dvector(prob,1,*n); free_dvector(denom,1,*n); free_dvector(suppz,1,*N); free_dvector(probz,1,*N); } /*------------------------------------------------------------- Hybrid monte carlo in generalized linear model hybridStep --------------------------------------------------------------*/ void hybridStep(char **glmFamily,double *Y,double *X, double *theta,double *ssize, double *beta,double *b,double *sigma, double *stepH,int *accept, int *ndata,int *xdim,int *ncycle,int *seed) { double **newX; double (*MEAN)(double t),(*MU)(double t); double *zz,*newzz,*newtheta,*newb,*grad,*diff,*newdiff; double total,stepcheck,logRatio=0.,pratio; int i,k,l; Y=Y-1;theta=theta-1;ssize=ssize-1;beta=beta-1;b=b-1;sigma=sigma-1; newX=dmatrix(1,*ndata,1,*xdim); zz=dvector(1,*xdim);newzz=dvector(1,*xdim); newtheta=dvector(1,*ndata); newb=dvector(1,*xdim); grad=dvector(1,*xdim); diff=dvector(1,*xdim);newdiff=dvector(1,*xdim); createMatrix(newX,X,*ndata,*xdim); /* Check glm family type */ if (strcmp(*glmFamily,"bernoulli") == 0 || strcmp(*glmFamily,"binomial") == 0 ){ MEAN=bernMean; MU=bernMU; } else if (strcmp(*glmFamily,"gamma") == 0){ MEAN=gammaMean; MU=gammaMU; } else if (strcmp(*glmFamily,"poisson") == 0){ MEAN=poissMean; MU=poissMU; } else nrerror("Misspecified density type"); /* Simulate standard normals for zz */ for (k=1;k<=*xdim;k++){ newb[k]=beta[k]; newzz[k]=zz[k]=rnormal(seed); grad[k]=0.; } /* First half-step for z */ for (i=1;i<=*ndata;i++){ for (k=1;k<=*xdim;k++){ grad[k]+=(Y[i]-(*MEAN)(theta[i])*ssize[i])*newX[i][k]; } } for (k=1;k<=*xdim;k++){ diff[k]=(beta[k]-b[k])/sqrt(sigma[k]); newzz[k]+=0.5*(*stepH)*(grad[k]-diff[k]/sqrt(sigma[k])); } /* Perform "ncycle" leap-frog discrete steps */ for (l=1;l<=*ncycle;l++){ for (k=1;k<=*xdim;k++){ grad[k]=0.; newb[k]+=(*stepH)*newzz[k]; } for (i=1;i<=*ndata;i++){ total=0.; for (k=1;k<=*xdim;k++) total+=(newb[k]-beta[k])*newX[i][k]; newtheta[i]=theta[i]+total; /*update theta*/ for (k=1;k<=*xdim;k++) grad[k]+=(Y[i]-(*MEAN)(newtheta[i])*ssize[i])*newX[i][k]; } stepcheck = (l<*ncycle) ? 1.0:0.5; for (k=1;k<=*xdim;k++){ newdiff[k]=(newb[k]-b[k])/sqrt(sigma[k]); newzz[k]+=stepcheck*(*stepH)*(grad[k]-newdiff[k]/sqrt(sigma[k])); } } /* Update beta? */ for (i=1;i<=*ndata;i++) logRatio += ( Y[i]*(newtheta[i]-theta[i]) -((*MU)(newtheta[i])-(*MU)(theta[i]))*ssize[i] ); logRatio+=-0.5*( vecdot(newdiff,newdiff,*xdim)-vecdot(diff,diff,*xdim) + vecdot(newzz,newzz,*xdim)-vecdot(zz,zz,*xdim) ); pratio=exp(logRatio); if (dran1(seed) <= DMIN(pratio,1.0)){ for (k=1;k<=*xdim;k++) beta[k]=newb[k]; for (i=1;i<=*ndata;i++) theta[i]=newtheta[i]; /*update theta */ *accept=1; } free_dmatrix(newX,1,*ndata,1,*xdim); free_dvector(zz,1,*xdim);free_dvector(newzz,1,*xdim); free_dvector(newtheta,1,*ndata); free_dvector(newb,1,*xdim); free_dvector(grad,1,*xdim); free_dvector(diff,1,*xdim);free_dvector(newdiff,1,*xdim); } /*------------------------------------------------------------- Sample kseq from generalized linear mixed model kseqglm --------------------------------------------------------------*/ void kseqglm(double *ldFam,double *p,int *kseq,int *Case, int *ndata,int *ncase,int *N,int *seed) { double **newldFam,**prob,*newp,*supp; int i,j; p=p-1;kseq=kseq-1;Case=Case-1; newldFam=dmatrix(1,*ndata,1,*N); prob=dmatrix(1,*ncase,1,*N); newp=dvector(1,*N); supp=dvector(1,*N); createMatrix(newldFam,ldFam,*ndata,*N); /* set up support vector + initialize probabilites */ for (j=1;j<=*N;j++){ supp[j]=j; for (i=1;i<=*ncase;i++) prob[i][j]=0.; } /* work out probabilities */ for (i=1;i<=*ndata;i++){ for (j=1;j<=*N;j++) prob[Case[i]][j] += newldFam[i][j]; } /* sample kseq */ for (i=1;i<=*ncase;i++){ for (j=1;j<=*N;j++) newp[j] = p[j]*exp(prob[i][j]); kseq[i]=rdisc(supp,newp,*N,seed); } free_dmatrix(newldFam,1,*ndata,1,*N); free_dmatrix(prob,1,*ncase,1,*N); free_dvector(newp,1,*N); free_dvector(supp,1,*N); } /*------------------------------------------------------------- Sample kseq from normal mixture model kseqnorm --------------------------------------------------------------*/ void kseqnorm(double *p,double *x, double *m, double *tau, int *kseq,int *n,int *N,int *seed) { double *newp,*supp; int i,k; p=p-1;x=x-1;m=m-1;tau=tau-1;kseq=kseq-1; newp=dvector(1,*N); supp=dvector(1,*N); /* set up support vector */ for (k=1;k<=*N;k++) supp[k]=k; /* sample kseq */ for (i=1;i<=*n;i++){ for (k=1;k<=*N;k++) newp[k]=p[k]*exp(-0.5*DSQR(x[i]-m[k])/tau[k])/sqrt(tau[k]); kseq[i]=rdisc(supp,newp,*N,seed); } free_dvector(newp,1,*N); free_dvector(supp,1,*N); } /*------------------------------------------------------------- Sample kseq from Poisson process model kseqPoiss --------------------------------------------------------------*/ void kseqPoiss(double *W,double *x,double *y, double *Zx,double *Zy, double *theta1,double *theta2,double *taub, int *kseq,int *n,int *N,int *seed) { double *newp,*supp,*beta; int i,k; W=W-1;x=x-1;y=y-1;Zx=Zx-1;Zy=Zy-1;kseq=kseq-1; newp=dvector(1,*N); supp=dvector(1,*N); beta=dvector(1,*N); /* set up support vector + beta-star */ for (k=1;k<=*N;k++){ supp[k]=k; beta[k]=exp(-0.5*(DSQR(Zx[k])+DSQR(Zy[k]))/(*taub))/ (1+exp((*theta1)-0.5*(DSQR(Zx[k])+DSQR(Zy[k]))/(*taub))); } /* sample kseq */ for (i=1;i<=*n;i++){ for (k=1;k<=*N;k++) newp[k]=W[k]*exp(-0.5*(DSQR(x[i]-Zx[k])+DSQR(y[i]-Zy[k]))/(*theta2))*beta[k]; kseq[i]=rdisc(supp,newp,*N,seed); } free_dvector(newp,1,*N); free_dvector(supp,1,*N); free_dvector(beta,1,*N); } /*------------------------------------------------------------- Metropolis for random effects in generalized linear model metglmStep --------------------------------------------------------------*/ void metglmStep(char **glmFamily,double *X,double *V, double *UB,double *ssize, double *Z,double *M,double *tau, double *stepM, int *accept,int *casept, int *ndata,int *nk,int *nran,int *ncycle,int *seed) { double **newV,**newZ,**ZZ; double (*MU)(double t),*logRatio; double oldtheta,newtheta,pratio,oldtotal,newtotal; int i,j,k,l; X=X-1;UB=UB-1;ssize=ssize-1;M=M-1;tau=tau-1;casept=casept-1; newV=dmatrix(1,*ndata,1,*nran); newZ=dmatrix(1,*nk,1,*nran); ZZ=dmatrix(1,*nk,1,*nran); logRatio=dvector(1,*nk); createMatrix(newV,V,*ndata,*nran); createMatrix(newZ,Z,*nk,*nran); /* Check glm family type */ if (strcmp(*glmFamily,"bernoulli") == 0 || strcmp(*glmFamily,"binomial") == 0 ){ MU=bernMU; } else if (strcmp(*glmFamily,"gamma") == 0){ MU=gammaMU; } else if (strcmp(*glmFamily,"poisson") == 0){ MU=poissMU; } else nrerror("Misspecified density type"); for (l=1;l<=*ncycle;l++){ for (j=1;j<=*nk;j++){ logRatio[j]=0; for (k=1;k<=*nran;k++) ZZ[j][k]=(*stepM)*rnormal(seed); } for (i=1;i<=*ndata;i++){ oldtotal=newtotal=0.; for (k=1;k<=*nran;k++){ oldtotal += newZ[casept[i]][k]*newV[i][k]; newtotal += (newZ[casept[i]][k]+ZZ[casept[i]][k])*newV[i][k]; } oldtheta=UB[i]+oldtotal; newtheta=UB[i]+newtotal; logRatio[casept[i]] += ( X[i]*(newtotal-oldtotal) -((*MU)(newtheta)-(*MU)(oldtheta))*ssize[i] ); } for (j=1;j<=*nk;j++){ for (k=1;k<=*nran;k++){ logRatio[j] += -0.5*(DSQR(newZ[j][k]+ZZ[j][k]-M[k])/tau[k]- DSQR(newZ[j][k]-M[k])/tau[k]); } /* Update Z? */ pratio=exp(logRatio[j]); if (dran1(seed) <= DMIN(pratio,1.0)){ for (k=1;k<=*nran;k++) newZ[j][k] += ZZ[j][k]; (*accept) += 1; } } } createVector(Z,newZ,*nk,*nran); free_dmatrix(newV,1,*ndata,1,*nran); free_dmatrix(newZ,1,*nk,1,*nran); free_dmatrix(ZZ,1,*nk,1,*nran); free_dvector(logRatio,1,*nk); } /*------------------------------------------------------------- Exact sampling for normal in random effects generalized linear model normalglm --------------------------------------------------------------*/ void normalglm(double *VX,double *V, double *Z,double *M,double *tau, int *casept,int *kseqn, int *ndata,int *N,int *nran,int *seed) { double ***var; double **newVX,**newV,**newZ,**mean,**tempvar; double *tempmean,*temp,*zz; int i,k,l,m; M=M-1;tau=tau-1;casept=casept-1;kseqn=kseqn-1; var=d3tensor(1,*nran,1,*nran,1,*N); newVX=dmatrix(1,*nran,1,*ndata); newV=dmatrix(1,*nran,1,*ndata); newZ=dmatrix(1,*N,1,*nran); mean=dmatrix(1,*nran,1,*N); tempvar=dmatrix(1,*nran,1,*nran); tempmean=dvector(1,*nran); temp=dvector(1,*nran); zz=dvector(1,*nran); createMatrix(newVX,VX,*nran,*ndata); createMatrix(newV,V,*nran,*ndata); /* Initialize mean + var */ for (k=1;k<=*N;k++){ if (kseqn[k]>0){ for (l=1;l<=*nran;l++){ mean[l][k]=M[l]/tau[l]; for (m=l;m<=*nran;m++){ if (m==l) var[l][m][k]=1./tau[l]; else var[l][m][k]=var[m][l][k]=0.; } } } else{ for (l=1;l<=*nran;l++){ mean[l][k]=M[l]; for (m=l;m<=*nran;m++){ if (m==l) var[l][m][k]=tau[l]; else var[l][m][k]=var[m][l][k]=0.; } } } } /* Calculate mean + var. Simulate mnormal */ for (i=1;i<=*ndata;i++){ for (l=1;l<=*nran;l++){ mean[l][casept[i]] += newVX[l][i]; for (m=1;m<=*nran;m++) var[l][m][casept[i]] += newV[l][i]*newV[m][i]; } } for (k=1;k<=*N;k++){ for (l=1;l<=*nran;l++){ tempmean[l]=mean[l][k]; for (m=l;m<=*nran;m++) tempvar[l][m]=tempvar[m][l]=var[l][m][k]; } if (kseqn[k]>0){ matinv(tempvar,tempvar,*nran); matdot(tempvar,tempmean,temp,*nran,*nran,*nran); rmnormal(temp,tempvar,zz,*nran,seed); } else rmnormal(tempmean,tempvar,zz,*nran,seed); for (l=1;l<=*nran;l++) newZ[k][l]=zz[l]; } createVector(Z,newZ,*N,*nran); free_d3tensor(var,1,*nran,1,*ndata,1,*N); free_dmatrix(newVX,1,*nran,1,*ndata); free_dmatrix(newV,1,*nran,1,*ndata); free_dmatrix(newZ,1,*N,1,*nran); free_dmatrix(mean,1,*nran,1,*N); free_dmatrix(tempvar,1,*nran,1,*N); free_dvector(tempmean,1,*nran); free_dvector(temp,1,*nran); free_dvector(zz,1,*nran); } /*------------------------------------------------------------- Update mean and variance in normal mixture model for Polya urn Gibbs sampler yUpdateNorm --------------------------------------------------------------*/ void yUpdateNorm(char **equalVar,char **accelerate, double *x,double *m,double *tau,double *M, int *n,double *A,double *a,double *b, double *v1,double *v2,double *q0const, int *seed) { double *uniqm,*uniqt,*supp,*prob,*xx,*xx2,*mm,*ss,s; int *indx,*rank,*freq; int i,j,nclus,nct,ix; x=x-1;m=m-1;tau=tau-1; uniqm=dvector(1,*n);uniqt=dvector(1,*n); supp=dvector(1,(*n)+1);prob=dvector(1,(*n)+1); xx=dvector(1,*n);xx2=dvector(1,*n);mm=dvector(1,*n);ss=dvector(1,*n); indx=ivector(1,*n);rank=ivector(1,*n);freq=ivector(1,*n); if (strcmp(*equalVar,"T") != 0){/*unequal variance case*/ for (i=1;i<=((*n)+1);i++) supp[i]=i; for (i=1;i<=(*n);i++){ sort(*n,m,indx,rank,uniqm,freq,&nclus); nct=1; uniqt[nct]=tau[indx[1]]; for (j=2;j<=(*n);j++){/*determine unique tau*/ if (tau[indx[j]] != uniqt[nct]){ nct += 1; uniqt[nct]=tau[indx[j]]; } } prob[1]=((*b)+(*a)*nclus)*(*q0const)/( pow(1+DSQR(x[i]-(*M))/(4*(*v2)),0.5+(*v1))); for (j=1;j<=nclus;j++) prob[1+j]=(freq[j]-(*a))*exp(-0.5*DSQR(x[i]-uniqm[j])/uniqt[j])/( sqrt(2*PI*uniqt[j])); ix=rdisc(supp,prob,(nclus+1),seed); if (ix==1){ tau[i]=((*v2)+DSQR(x[i]-(*M))/4)/rGamma((*v1)+0.5,seed); m[i]=rnormal(seed)*sqrt(tau[i]/2)+((*M)+x[i])/2; } else{ m[i]=uniqm[ix-1]; tau[i]=uniqt[ix-1]; } } if (strcmp(*accelerate,"T") == 0){/*accelerate?*/ sort(*n,m,indx,rank,uniqm,freq,&nclus); for (i=1;i<=(*n);i++) xx[i]=xx2[i]=0; for (i=1;i<=(*n);i++){ xx[rank[i]] += x[i]; xx2[rank[i]] += DSQR(x[i]); } for (j=1;j<=nclus;j++){ ss[j]=((*v2)+DSQR(*M)/2+xx2[j]/2 -DSQR((*M)+xx[j])/(2*(freq[j]+1)))/rGamma((*v1)+freq[j]/2,seed); mm[j]=rnormal(seed)*sqrt(ss[j]/(freq[j]+1)) + ((*M)+xx[j])/(freq[j]+1); } for (i=1;i<=(*n);i++){ m[i]=mm[rank[i]]; tau[i]=ss[rank[i]]; } } } else{/*equal variance case*/ for (i=1;i<=((*n)+1);i++) supp[i]=i; s=tau[1]; for (i=1;i<=(*n);i++){ sort(*n,m,indx,rank,uniqm,freq,&nclus); prob[1]=((*b)+(*a)*nclus)*exp(0.5*DSQR(x[i]/s+(*M)/(*A))/(1/s+1/(*A)) -0.5*DSQR(x[i])/s-0.5*DSQR(*M)/(*A))/sqrt(2*PI*((*A)+s)); for (j=1;j<=nclus;j++) prob[1+j]=(freq[j]-(*a))*exp(-0.5*DSQR(x[i]-uniqm[j])/s)/( sqrt(2*PI*s)); ix=rdisc(supp,prob,(nclus+1),seed); if (ix==1) m[i]=rnormal(seed)/sqrt(1/s+1/(*A))+(x[i]/s+(*M)/(*A))/ (1/s+1/(*A)); else m[i]=uniqm[ix-1]; } if (strcmp(*accelerate,"T") == 0){/*accelerate?*/ sort(*n,m,indx,rank,uniqm,freq,&nclus); for (i=1;i<=(*n);i++) xx[i]=0; for (i=1;i<=(*n);i++) xx[rank[i]] += x[i]; for (j=1;j<=nclus;j++) mm[j]=rnormal(seed)/sqrt(freq[j]/s+1/(*A)) + (xx[j]/s+(*M)/(*A))/(freq[j]/s+1/(*A)); for (i=1;i<=(*n);i++) m[i]=mm[rank[i]]; } s=(*v2); for (i=1;i<=*n;i++) s += DSQR(x[i]-m[i])/2.; s=s/rGamma((*v1)+(*n)/2,seed); for (i=1;i<=*n;i++) tau[i]=s; } free_dvector(uniqm,1,*n);free_dvector(uniqt,1,*n); free_dvector(supp,1,(*n)+1);free_dvector(prob,1,(*n)+1); free_dvector(xx,1,*n);free_dvector(xx2,1,*n); free_dvector(mm,1,*n);free_dvector(ss,1,*n); free_ivector(indx,1,*n);free_ivector(rank,1,*n); free_ivector(freq,1,*n); } /*------------------------------------------------------------- Update mean and variance in normal mixture model for blocked Gibbs sampler zUpdateNorm --------------------------------------------------------------*/ void zUpdateNorm(char **equalVar,char **tauTrunc, double *x,double *m,double *tau,double *M, int *kseq,int *kseqn, int *n,int *N,double *A,double *a,double *b, int *seed) { double *mean,*var,*scale,*shape,sc,sh,tt; int i,k; x=x-1;m=m-1;tau=tau-1;kseq=kseq-1;kseqn=kseqn-1; mean=dvector(1,*N);var=dvector(1,*N); scale=dvector(1,*N);shape=dvector(1,*N); /* initialize mean and variances */ for (k=1;k<=*N;k++){ mean[k]=0; var[k]=1./(kseqn[k]/tau[k]+1./(*A)); } /* update m */ for (i=1;i<=*n;i++) mean[kseq[i]] += x[i]; for (k=1;k<=*N;k++) m[k] = rnormal(seed)*sqrt(var[k])+ var[k]*(mean[k]/tau[k]+(*M)/(*A)); /* initialize scale and shape */ if (strcmp(*equalVar,"T") == 0){ sc=(*b); sh=(*a)+(*n)/2.; for (i=1;i<=*n;i++) sc += DSQR(x[i]-m[kseq[i]])/2.; tt = sc/rGamma(sh,seed); for (k=1;k<=*N;k++) tau[k] = tt; } else if (strcmp(*tauTrunc,"T") == 0){ for (k=1;k<=*N;k++) tau[k]=0.; /* update truncation level (store in tau for now) */ for (i=1;i<=*n;i++) tau[kseq[i]] += DSQR(x[i]-m[kseq[i]])/2.; } else{ for (k=1;k<=*N;k++){ scale[k]=(*b); shape[k]=(*a)+kseqn[k]/2.; } /* update tau */ for (i=1;i<=*n;i++) scale[kseq[i]] += DSQR(x[i]-m[kseq[i]])/2.; for (k=1;k<=*N;k++) tau[k] = scale[k]/rGamma(shape[k],seed); } free_dvector(mean,1,*N);free_dvector(var,1,*N); free_dvector(scale,1,*N);free_dvector(shape,1,*N); } /*------------------------------------------------------------- Update Z and theta in Poisson spatial process mixture model for blocked Gibbs sampler zUpdatePoiss --------------------------------------------------------------*/ void zUpdatePoiss(double *x,double *y,double *Zx,double *Zy, double *theta1,double *theta2, double *t1V,double *tstep,double *accept, double *taua,double *taub, int *kseq,int *kseqn,int *n,int *N,int *seed) { double *meanx,*meany,*var,*Zxnew,*Zynew; double zold,znew,zsqr,metr,bold,bnew,t1new,t2new; int i,k; x=x-1;y=y-1;Zx=Zx-1;Zy=Zy-1;kseq=kseq-1;kseqn=kseqn-1; meanx=dvector(1,*N);meany=dvector(1,*N);var=dvector(1,*N); Zxnew=dvector(1,*N);Zynew=dvector(1,*N); /* initialize mean and variances */ for (k=1;k<=*N;k++){ meanx[k]=meany[k]=0; var[k]=((*taua)*(*theta2))/(kseqn[k]*(*taua)+(*theta2)); } /* update Zx, Zy using independence Metropolis */ for (i=1;i<=*n;i++){ meanx[kseq[i]] += x[i]; meany[kseq[i]] += y[i]; } for (k=1;k<=*N;k++){ Zxnew[k] = rnormal(seed)*sqrt(var[k])+meanx[k]*var[k]/(*theta2); Zynew[k] = rnormal(seed)*sqrt(var[k])+meany[k]*var[k]/(*theta2); zold=DSQR(Zx[k])+DSQR(Zy[k]); znew=DSQR(Zxnew[k])+DSQR(Zynew[k]); if (kseqn[k] == 0){ metr = exp(-log(1+exp((*theta1)-0.5*znew/(*taub)))/(*N) +log(1+exp((*theta1)-0.5*zold/(*taub)))/(*N)); if (dran1(seed) <= DMIN(metr,1.0)){ Zx[k]=Zxnew[k];Zy[k]=Zynew[k]; } } else{ bold = -0.5*kseqn[k]*zold/(*taub) -(kseqn[k]+1/(*N))*log(1+exp((*theta1)-0.5*zold/(*taub))); bnew = -0.5*kseqn[k]*znew/(*taub) -(kseqn[k]+1/(*N))*log(1+exp((*theta1)-0.5*znew/(*taub))); metr = exp(bnew-bold); if (dran1(seed) <= DMIN(metr,1.0)){ Zx[k]=Zxnew[k];Zy[k]=Zynew[k]; } } } /* update theta1 using RWMH. */ t1new = (*theta1)+rnormal(seed)*(*tstep); metr = (*n)*(t1new-(*theta1))+0.5*(DSQR(*theta1)-DSQR(t1new))/(*t1V); for (k=1;k<=*N;k++){ zsqr=DSQR(Zx[k])+DSQR(Zy[k]); metr += (kseqn[k]+1/(*N))*log(1+exp((*theta1)-0.5*zsqr/(*taub))) -(kseqn[k]+1/(*N))*log(1+exp(t1new-0.5*zsqr/(*taub))); } metr = exp(metr); if (dran1(seed) <= DMIN(metr,1.0)){ (*theta1)=t1new; (*accept) += 1; } free_dvector(meanx,1,*N);free_dvector(meany,1,*N);free_dvector(var,1,*N); free_dvector(Zxnew,1,*N);free_dvector(Zynew,1,*N); } /************************************************************* * nrutil.h routines * ***************************************************************/ /*------------------------------------------------------------- allocate a double matrix a[nrl..nrh][ncl..nch] that points to the matrix declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 and ncol=nch-ncl+1. The routine should be called with the address &a[0][0] as the first argument. convert_dmatrix --------------------------------------------------------------*/ double **convert_dmatrix(double *a,long nrl,long nrh,long ncl,long nch) { long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; double **m; /* allocate pointers to rows */ m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*))); if (!m) nrerror("allocation failure in convert_dmatrix()"); m += NR_END; m -= nrl; /* set pointers to rows */ m[nrl]=a-ncl; for(i=1,j=nrl+1;i