automixprime.c

he AutoMix package is a C program for Unix-like systems, implementing the automatic reversible jum

      }
      else{
         l=0;
         palloc[l]=1.0;
      }
      
      
   /* --Section 9.2 - Standardise state variable --------------- */
      
      for(j1=0;j1<nkk;j1++){
         work[j1]=theta[j1]-mu[k][l][j1];
      }
      for(j1=0;j1<nkk;j1++){
         for(j2=0;j2<j1;j2++){
            work[j1]=work[j1]-B[k][l][j1][j2]*work[j2];
         }
         work[j1]=work[j1]/B[k][l][j1][j1];
      }
      
   /* --Section 9.3 - Choose proposed new model and component ----*/
      
      if(kmax==1){
         kn=k;
         logratio=0.0;
      }
      else{
         gamma=pow(1.0/(sweep+1),(2.0/3.0));
         u=sdrand();
         thresh=0.0;
         for(k1=0;k1<kmax;k1++){
            thresh+=pk[k1];
            if(u<thresh){
               kn=k1;
               break;
            }
         }
         logratio=log(pk[k])-log(pk[kn]);
      }
      
      nkkn=nk[kn];
      Lkkn=Lk[kn];
      
      u=sdrand();
      thresh=0.0;
      for(l1=0;l1<Lkkn;l1++){
         thresh+=lambda[kn][l1];
         if(u<thresh){
            ln=l1;
            break;
         }
      }
      
   /* --Section 9.4 Propose new state ----------------*/
      
      if(nkk<nkkn){
         if(dof>0){
            rt(&(work[nkk]),nkkn-nkk,dof);
            for(j1=nkk;j1<nkkn;j1++){
               logratio-=ltprob(dof,work[j1],&constt);
            }
         }
         else{
            gauss(&(work[nkk]),nkkn-nkk);
            for(j1=nkk;j1<nkkn;j1++){
               logratio+=0.5*pow(work[j1],2.0)+logrtpi;
            }
         }
         if(doperm){
            perm(work,nkkn);
         }
      }
      else if(nkk==nkkn){
         if(doperm){
            perm(work,nkk);
         }
      }
      else{
         if(doperm){
            perm(work,nkk);
         }
         if(dof>0){
            for(j1=nkkn;j1<nkk;j1++){
               logratio+=ltprob(dof,work[j1],&constt);
            }
         }
         else{
            for(j1=nkkn;j1<nkk;j1++){
               logratio-=(0.5*pow(work[j1],2.0)+logrtpi);
            }
         }
      }
      
      for(j1=0;j1<nkkn;j1++){
         thetan[j1]=mu[kn][ln][j1];
         for(j2=0;j2<=j1;j2++){
            thetan[j1]+=B[kn][ln][j1][j2]*work[j2];
         }
      }
      
   /* --Section 9.5 - Work out probability of allocating to component
      for acceptance ratio (reverse move) ---------*/
      
      if(Lkkn>1){
         sum=0.0;
         for(l1=0;l1<Lkkn;l1++){
            pallocn[l1]=log(lambda[kn][l1])+lnormprob(kn,nkkn,l1,mu,B,thetan);
            pallocn[l1]=exp(pallocn[l1]);
            sum+=pallocn[l1];
         }
         if(sum>0){
            for(l1=0;l1<Lkkn;l1++){
               pallocn[l1]/=sum;
            }
         }
         else{
            for(l1=0;l1<Lkkn;l1++){
               pallocn[l1]=1.0/Lkkn;
            }
         }
      }
      else{
         pallocn[ln]=1.0;
      }
      
   /* --Section 9.6 - Work out acceptance probability	and new state --*/
      
      lpn=lpost(kn,nkkn,thetan,&llhn);
      
      logratio+=(lpn-lp);
      logratio+=(log(pallocn[ln])-log(palloc[l]));
      logratio+=(log(lambda[k][l])-log(lambda[kn][ln]));
      logratio+=(log(detB[kn][ln])-log(detB[k][l]));
      
      if(sdrand()<exp(max(-30.0,min(0.0,logratio)))){
         for(j1=0;j1<nkkn;j1++){
            theta[j1]=thetan[j1];
         }
         lp=lpn;
         llh=llhn;
         k=kn;
         nkk=nkkn;
         Lkk=Lkkn;
         nacctd++;
      }
      
      if(adapt==1){
         if(sweep>nburn){
            for(k1=0;k1<kmax;k1++){
               if(k1==k){
                  propk[k1]=1.0;
               }
               else{
                  propk[k1]=0.0;
               }
               pk[k1]+=(gamma*(propk[k1]-pk[k1]));
            }
            for(k1=0;k1<kmax;k1++){
               if(pk[k1]<pkllim){
                  reinit=1;
               }
            }
            if(reinit==1){
               reinit=0;
               nreinit++;
               pkllim=1.0/(10.0*nreinit);
               for(k1=0;k1<kmax;k1++){
                  pk[k1]=1.0/kmax;
               }
            }
         }
      }
      
      if(sweep>nburn){
         (ksummary[k])++;
      }
   /* --- Section 10 - Write variables to files --------- */
      
      if(sweep>nburn){
         
         fprintf(fpk,"%d\n",k+1);
         fprintf(fplp,"%lf %lf\n",lp,llh);
         for(k1=0;k1<kmax;k1++){
            fprintf(fpp,"%lf ",pk[k1]);
         }
         fprintf(fpp,"\n");
         for(j1=0;j1<nkk;j1++){
            fprintf(fpt[k],"%lf ",theta[j1]);
         }
         fprintf(fpt[k],"\n");
      }
      
      if(sweep>keep&&fmod(sweep-keep,nsokal)<0.005){
         xr[((sweep-keep)/nsokal)-1]=k;
      }
      
      if(sweep==1){
         printf("\nBurning in");
      }
      if((sweep<=nburn)&&(fmod(sweep,(nburn/10))<tol)){
         printf(" .");
         fflush(NULL);
      }
      
      if(sweep==(nburn+1)){
         printf("\nStart of main sample:");
      }
      
      if((sweep>nburn)&&(fmod(sweep-nburn,(nsweep/10))<tol)){
         printf("\nNo. of iterations remaining: %d",nsweep+nburn-sweep);
      }
      fflush(NULL);
   }
   printf("\n");
   
   
  /* --- Section 11 - Write log file ----------------------*/
   
   sokal(nkeep,xr,&var,&tau,&m);
   fprintf(fpl,"\nAutocorrelation Time:\n");
   fprintf(fpl,"nkeep:%d, nsokal:%d, var:%lf, tau:%lf\n",nkeep,nsokal,
   var,tau);
   for(i1=0;i1<m;i1++){
      fprintf(fpac,"%lf\n",xr[i1]);
   }
   
   fprintf(fpl,"\nPosterior Model Probabilities:\n");
   for(k1=0;k1<kmax;k1++){
      fprintf(fpl,"Model %d: %lf\n",k1+1,(double)ksummary[k1]/(double)nsweep);
   }
   
   fprintf(fpl,"\nAcceptance Rates:\n");
   fprintf(fpl,"Block RWM: %lf\n",(double)naccrwmb/(double)ntryrwmb);
   fprintf(fpl,"Single RWM: %lf\n",(double)naccrwms/(double)ntryrwms);
   fprintf(fpl,"Auto RJ: %lf\n",(double)nacctd/(double)ntrytd);
   
   endtime=clock();
   timesecs=(endtime-starttime)/((double)CLOCKS_PER_SEC);
   fprintf(fpl,"\nRun time:\n");
   fprintf(fpl,"Time: %lf\n",timesecs);
   
   return 0;
   
}


void gauss(double *z,int n){
   
  /* Uses Box mueller method to simulate n N(0,1) variables and stores them
    in z */
   
   int i,n1;
   double u,v;
   
   n1=n-1;
   for(i=0;i<n1;i+=2){
      u=sqrt(-2.0*log(sdrand()));
      v=tpi*sdrand();
      z[i]=u*sin(v);
      z[i+1]=u*cos(v);
   }
   if(fmod(n,2)<0.5){
      return;
   }
   else{
      u=sqrt(-2.0*log(sdrand()));
      v=tpi*sdrand();
      z[n-1]=u*sin(v);
   }
   return;
}

void rt(double *z,int n,int dof){
   
  /* Simulates n random t variable with dof degrees of freedom
    by simulating standard normals and chi-squared random variables.
    Chi-squared rvs simulated by rgamma function that simulates random gamma
    variables (see gammafns file for details)*/
   
   int j1;
   double s;
   
   gauss(z,n);
   s=0.5*dof;
   for(j1=0;j1<n;j1++){
      z[j1]/=sqrt(rgamma(s)/s);
   }
   return;
}

void chol(int nkk, double **A){
   
  /* Performs cholesky decompositon of A and returns result in the
    same matrix - adapted from PJG Fortran function*/
   
   int j1,j2,j3;
   double sum;
   
   for(j1=0;j1<nkk;j1++){
      sum=A[j1][j1];
      for(j2=0;j2<j1;j2++){
         sum-=pow(A[j1][j2],2);
      }
      A[j1][j1]=sqrt(sum);
      
      for(j2=j1+1;j2<nkk;j2++){
         sum=A[j2][j1];
         for(j3=0;j3<j1;j3++){
            sum-=A[j2][j3]*A[j1][j3];
         }
         A[j2][j1]=sum/A[j1][j1];
      }
   }
}

void perm(double *work, int nkk){
   
  /* Randomly permutes the nkk-vector work */
   int j1,j2;
   double temp;
   
   for(j1=0;j1<(nkk-1);j1++){
      j2=j1+(int)((nkk-j1)*sdrand());
      if(j2!=j1){
         temp=work[j2];
         work[j2]=work[j1];
         work[j1]=temp;
      }
   }
   return;
}

double ltprob(int dof,double z,double *constt){
   
  /* Evaluates the log of p.d.f. of a t variable with dof degrees of freedom
    at point z */
   
   double out;
   
  /* only calculate const of proportionality once */
   if((*constt)>10000.0){
      *constt=loggamma(0.5*(dof+1))-loggamma(0.5*dof)-0.5*log(dof*pi);
   }
   out=(*constt)-0.5*(dof+1)*log(1.0+pow(z,2.0)/dof);
   return out;
}


double lnormprob(int k,int nkk,int l,double ***mu,
double ****B, double *datai){
   
  /* Evaluates log of p.d.f. for a multivariate normal for model
    k, of dimension nkk, component l. The summary of means and
    sqrt of cov matrices (for all models and all component)
    are supplied in mu and B */
   
   int j1,j2;
   double* work;
   double out;
   
   work = (double*)malloc(sizeof(double)*nkk);
   
   for(j1=0;j1<nkk;j1++){
      work[j1]=datai[j1]-mu[k][l][j1];
   }
   for(j1=0;j1<nkk;j1++){
      for(j2=0;j2<j1;j2++){
         (work[j1])-=B[k][l][j1][j2]*work[j2];
      }
      (work[j1])/=B[k][l][j1][j1];
   }
   out=0.0;
   for(j1=0;j1<nkk;j1++){
      out+=(work[j1]*work[j1]);
   }
   out=-0.5*out-(nkk/2.0)*log(tpi)-log(det(k,nkk,l,B));
   
   free(work);
   return out;
   
}

double det(int k, int nkk, int l, double ****B){
   
  /* Evaluates the determinant of a matrix in B corresponding to model k,
    component l. */
   int j1;
   double out;
   out=1.0;
   for(j1=0;j1<nkk;j1++){
      out*=B[k][l][j1][j1];
   }
   return out;
}