automix.c

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

/* The AutoMix program.
 
   Last edited 25/11/04.
   Developed by David Hastie, Department of Mathematics,
   University of Bristol, UK as a part of a submission for the
   degree of Ph.D. This Ph.D. was supervised by Prof. Peter Green (PJG),
   University of Bristol. Special thanks also to Dr. Christophe Andrieu CA),
   University of Bristol, for advice on adaptive schemes and mixture fitting.
 
   The AutoMix sampler is free for personal and academic use, but must
   reference the sampler as instructed below.  For commercial
   use please permission must be sought from the author. To seek permission
   for such use please send an e-mail to [e-mail]
   outlining the desired usage.
 
   Use of the AutoMix sampler is entirely at the user's own risk. It is the
   responsibility of the user to ensure that any conclusions made through the
   use of the AutoMix sampler are valid. The author accepts no responsibility
   whatsoever for any loss, financial or otherwise, that may arise in
   connection with the use of the sampler.
 
   The AutoMix sampler is available from http://www.davidhastie.me.uk/AutoMix
   Although the sampler may be modified and redistributed, the author
   encourages users to register at the above site so that updates of the
   software can be received.
 
   Before use, please read the README file bundled with this software.
 
   Users should reference the sampler as instructed on the AutoMix website
   (see above). Initially this is likely to be the Ph.D. thesis that
   introduces the AutoMix sampler. However, this will hopefully change to
   be a published paper in the not too distant future.	*/

/* Standard library files */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <fcntl.h>
#include <time.h>

#define MAX(A,B) ((A)>(B)?(A):(B))
#define MIN(A,B) ((A)<(B)?(A):(B))

/* Global constants (please feel free to change as required)
 
   nkmaxmax = maximum dimension of any one model under consideration
   kmaxmax = maximum number of models
   Lkmaxmax = initial number of mixture components fitted in stage 2 of
           AutoMix algorithm */

#define nkmaxmax 20
#define kmaxmax 15
#define Lkmaxmax 30
#define tpi 6.283185307179586477
#define pi 3.141592653589793238
#define logrtpi 0.5*log(tpi)

/* --- Internal functions (described below) ----------------- */

void gauss(double *z,int nkk);

void rt(double *z,int nkk,int dof);

void chol(int nkk,double **B);

void perm(double *work,int nkk);

double ltprob(int dof,double z,double *constt);

double lnormprob(int k,int nkk,int l,double ***mu,
double ****B, double *datai);

double det(int k, int nkk, int l, double ****B);

/* --- External functions -----------------*/

/* Random number functions sdrni and sdrand can be found in sd.c
   bundled with this software. References for these functions can be
   found within sd.c */

extern void sdrni(unsigned long *seed);

extern double sdrand();

/* Functions rgamma and loggamma can be found in gammafns.c bundled with
   with this software. */

extern double rgamma(double s);

extern double loggamma(double s);

/* Function sokal is found in file sokal.c bundled with this software. */

extern void sokal(int n, double *xreal, double *var, double *tau, int *m);

/* --- User supplied functions ------------ */

/* Functions must be supplied in user***.c file (see e.g. usertoy1.c,
   usercpt.c etc bundled with this software).
 
   Descriptions:
   1. lpost(&k,theta,&llh)
   This should be a c function written by the user that evaluates
   log posterior (up to an additive constant) at (k,theta). The function
   can also return the likelihood at this point in llh.
   2. getkmax(&kmax)
   This should be a c function written by the user that returns the
   number of models kmax.
   3. getnk(kmax, nk)
   This should be a c function written by the user that returns the dimensions
   nk for model k=1,...,kmax.
   4. getic(k,nkk,rwm)
   This should be a c function written by the user that returns the
   possibly random starting point for the rwm in stage 1 of the AutoMix
   sampler */

extern double lpost(int k, int nkk, double *theta, double *llh);

extern void getkmax(int *kmax);

extern void getnk(int kmax,int* nk);

extern void getic(int k, int nkk, double *rwm);

/* ---main program-------------------------- */

int main(int argc,char *argv[]){
   
  /*---Section 1 Declare Variables -------------------------*/
   
  /* ---clock variables ---------------------- */
   clock_t starttime,endtime;
   double timesecs;
   
  /* ---indexing variables ------------------- */
   int t1,t2,i1,i2,j1,j2,k1,l1,l2,sweep,remain;
   
  /* ---counting variables ------------------- */
   int nsweep,count,nsweep2,naccrwmb,naccrwms,nacctd,ntryrwmb,ntryrwms,ntrytd;
   int nburn,nsokal,nkeep,keep,nsweepr,*nacc,*ntry,*ksummary;
   
  /* ---command line reading parameters ------ */
   int numargs,sametest;
   char word[20],selector[3],iparam[18];
   
  /* ---filename variables ------------------- */
   int check;
   FILE *fpk,*fpl,*fpt[kmaxmax],*fpcf,*fpmix,*fplp,*fpp,*fpac,*fpad;
   char fname[18],fname1[18],kno[6];
   
  /* ---random no. variables ----------------- */
   unsigned long seed;
   double u,constt;
   int dof;
   
  /* ---logical variables -------------------- */
   int doperm;
   
  /* ---State parameters and variables ------- */
   int k,kn,kmax,nkk,nkkn,nkmax,lendata;
   int *nk;
   double *theta,*thetan;
   double **data,*propk,*pk;
   
  /* ---Mixture parameters --------------------*/
   int l,ln,Lkk,Lkkn,Lkkmin,nparams,ldel,Lkmax;
   int *Lk;
   double ***mu,****B,****BBT,**detB;
   double ***mumin,****Bmin,****BBTmin;
   double **lambda,**lambdamin,**w,*logw,**lpdatagivenl,*palloc,*pallocn;
   double tol=0.00001,costfn,costfnnew,costfnmin,minlambda;
   
  /* ---RWM parameters ------------------------*/
   double *rwm,*rwmn,Z[1],*Znkk;
   double **sig,gamma,accept,alphastar=0.25;
   
  /* ---Probabilities ------------------------ */
   double lp,lpn,logratio,llh,llhn;
   
  /* ---working arrays and variables --------- */
   int indic,stop,natann,forceann,mode;
   int *init;
   double sum,sigma,wnew,wnewl1,*sumw,sumwnew,sumlambda,*work,thresh;
   double *datamean,**M1;
   
  /* ---autocorrelation variables ------------ */
   double *xr,var,tau;
   int m;
   
  /* ---adaptation parameters ---------------- */
   int adapt,reinit,nreinit;
   double pkllim;
   
   
  /* --- Section 2 - Read in Comand Line Variables ----------------- */
   
   starttime=clock();
   
  /* Definition of command line variables and explanation
   
    Prog variable ~ Command line variable ~ Explanation
   
    mode ~ m ~ 0 if mixture fitting, 1 if user supplied mixture params,
            2 if AutoRJ
    nsweep ~ N ~ no. of reversible jump sweeps in stage 3
    nsweep2 ~ n ~ max(n,10000*nk,100000) sweeps in within-model RWM in stage 1
    adapt ~ a ~ 1 if RJ adaptation done in stage 3, 0 otherwise
    doperm ~ p ~ 1 if random permutation done in stage 3 RJ move, 0 otherwise
    seed ~ s ~ random no. seed, 0 uses clock seed
    dof ~ t ~ 0 if Normal random variables used in RWM and RJ moves, otherwise
            specify integer degrees of freedom of student t variables
    fname ~ f ~ filename base */
   
   
  /* Default values */
   
   nsweep=100000;
   nsweep2=100000;
   numargs=argc-1;
   strcpy(fname,"output");
   doperm=0;
   seed=0;
   mode=0;
   adapt=0;
   dof=0;
   
  /* Override defaults if user supplies command line options */
   
   if(numargs>0){
      for(t1=1;t1<=numargs;t1++){
         
         strcpy(word,argv[t1]);
         for(t2=0;t2<2;t2++){
            selector[t2]=word[t2];
         }
         selector[2]='\0';
         for(t2=0;t2<17;t2++){
            iparam[t2]=word[t2+2];
         }
         iparam[17]='\0';
         
         
         sametest=strcmp(selector,"-f");
         if(sametest==0){
            strcpy(fname,iparam);
            continue;
         }
         sametest=strcmp(selector,"-N");
         if(sametest==0){
            nsweep=atoi(iparam);
            continue;
         }
         
         sametest=strcmp(selector,"-n");
         if(sametest==0){
            nsweep2=max(atoi(iparam),100000);
            continue;
         }
         
         sametest=strcmp(selector,"-s");
         if(sametest==0){
            seed=atoi(iparam);
            continue;
         }
         sametest=strcmp(selector,"-p");
         if(sametest==0){
            doperm=atoi(iparam);
            continue;
         }
         sametest=strcmp(selector,"-m");
         if(sametest==0){
            mode=atoi(iparam);
            continue;
         }
         
         sametest=strcmp(selector,"-a");
         if(sametest==0){
            adapt=atoi(iparam);
            continue;
         }
         
         sametest=strcmp(selector,"-t");
         if(sametest==0){
            dof=atoi(iparam);
            continue;
         }
         
      }
   }
   
   sdrni(&seed);
   
   
  /* --- Section 3 - Initial File handling ---------------------  */
   
   sprintf(fname1,fname);
   strcat(fname1,"_log.data");
   fpl = fopen(fname1,"w");
   sprintf(fname1,fname);
   strcat(fname1,"_pk.data");
   fpp = fopen(fname1,"w");
   sprintf(fname1,fname);
   strcat(fname1,"_ac.data");
   fpac = fopen(fname1,"w");
   sprintf(fname1,fname);
   strcat(fname1,"_adapt.data");
   fpad = fopen(fname1,"w");
   sprintf(fname1,fname);
   strcat(fname1,"_cf.data");
   fpcf = fopen(fname1,"w");
   
  /* Print user options to log file */
   
   fprintf(fpl,"seed: %ld\n",seed);
   fprintf(fpl,"m: %d\n",mode);
   fprintf(fpl,"a: %d\n",adapt);
   fprintf(fpl,"p: %d\n",doperm);
   fprintf(fpl,"n: %d\n",nsweep2);
   fprintf(fpl,"N: %d\n",nsweep);
   
  /* Check user has supplied mixture parameters if trying to use mode 1.
    If not default back to mode 0 */
   
   if(mode==1){
      sprintf(fname1,fname);
      strcat(fname1,"_mix.data");
      if((fpmix = fopen(fname1,"r"))==NULL){
         printf("\nMixture file doesn't exist:");
         printf("\nContinuing using RWM to estimate parameters");
         mode=0;
      }
   }
   
  /* --- Section 4.0 - Read in key variables from user functions -*/
   
   getkmax(&kmax);
   if(kmax>kmaxmax){
      printf("\nError:kmax too large \n");
      return 0;
   }
   else if(kmax<0){
      printf("\nError:negative kmax \n");
      return 0;
   }
   
   nk=(int*)malloc(kmax*sizeof(int));
   Lk=(int*)malloc(kmax*sizeof(int));
   ksummary=(int*)malloc(kmax*sizeof(int));
   
   getnk(kmax,nk);
   nkmax=nk[0];
   ksummary[0]=0;
   for(k1=1;k1<kmax;k1++){
      nkmax=max(nk[k1],nkmax);
      ksummary[k1]=0;
   }
   
   lambda=(double**)malloc(kmax*sizeof(double));
   lambdamin=(double**)malloc(kmax*sizeof(double));
   mu=(double***)malloc(kmax*sizeof(double));
   mumin=(double***)malloc(kmax*sizeof(double));
   BBT=(double****)malloc(kmax*sizeof(double));
   BBTmin=(double****)malloc(kmax*sizeof(double));
   B=(double****)malloc(kmax*sizeof(double));
   Bmin=(double****)malloc(kmax*sizeof(double));
   detB=(double**)malloc(kmax*sizeof(double));
   sig=(double**)malloc(kmax*sizeof(double));
   for(k1=0;k1<kmax;k1++){
      nkk=nk[k1];
      lambda[k1]=(double*)malloc(Lkmaxmax*sizeof(double));
      lambdamin[k1]=(double*)malloc(Lkmaxmax*sizeof(double));
      mu[k1]=(double**)malloc(Lkmaxmax*sizeof(double));
      mumin[k1]=(double**)malloc(Lkmaxmax*sizeof(double));
      BBT[k1]=(double***)malloc(Lkmaxmax*sizeof(double));
      BBTmin[k1]=(double***)malloc(Lkmaxmax*sizeof(double));
      B[k1]=(double***)malloc(Lkmaxmax*sizeof(double));
      Bmin[k1]=(double***)malloc(Lkmaxmax*sizeof(double));
      detB[k1]=(double*)malloc(Lkmaxmax*sizeof(double));
      sig[k1]=(double*)malloc(nkk*sizeof(double));
      for(l1=0;l1<Lkmaxmax;l1++){
         mu[k1][l1]=(double*)malloc(nkk*sizeof(double));
         mumin[k1][l1]=(double*)malloc(nkk*sizeof(double));
         BBT[k1][l1]=(double**)malloc(nkk*sizeof(double));
         BBTmin[k1][l1]=(double**)malloc(nkk*sizeof(double));
         B[k1][l1]=(double**)malloc(nkk*sizeof(double));
         Bmin[k1][l1]=(double**)malloc(nkk*sizeof(double));
         for(j1=0;j1<nkk;j1++){
            BBT[k1][l1][j1]=(double*)malloc(nkk*sizeof(double));
            BBTmin[k1][l1][j1]=(double*)malloc(nkk*sizeof(double));
            B[k1][l1][j1]=(double*)malloc(nkk*sizeof(double));
            Bmin[k1][l1][j1]=(double*)malloc(nkk*sizeof(double));
         }
      }
   }
   
   
  /* --- Section 5.1 - Read in mixture parameters if mode 1 (m=1) --- */
   
  /* These parameters are used if mode 1 (m=1) of the AutoMix sampler is