From 36cd7fc4fe847584a959584afdb0383e82e089c4 Mon Sep 17 00:00:00 2001 From: Martin Beroiz Date: Mon, 16 Dec 2019 13:07:49 -0600 Subject: [PATCH] Original AutoMix by David Hastie --- Makefile | 201 +++++ README.txt | 333 +++++++ automix.c | 1599 ++++++++++++++++++++++++++++++++++ ddidata.h | 2366 ++++++++++++++++++++++++++++++++++++++++++++++++++ gammafns.c | 338 ++++++++ register.php | 81 ++ sd.c | 157 ++++ sokal.c | 305 +++++++ usercpt.c | 140 +++ usercptrs.c | 134 +++ userddi.c | 809 +++++++++++++++++ userrb9.c | 206 +++++ usertoy1.c | 126 +++ usertoy2.c | 84 ++ 14 files changed, 6879 insertions(+) create mode 100644 Makefile create mode 100644 README.txt create mode 100644 automix.c create mode 100644 ddidata.h create mode 100644 gammafns.c create mode 100644 register.php create mode 100644 sd.c create mode 100644 sokal.c create mode 100644 usercpt.c create mode 100644 usercptrs.c create mode 100644 userddi.c create mode 100644 userrb9.c create mode 100644 usertoy1.c create mode 100644 usertoy2.c diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..2647875 --- /dev/null +++ b/Makefile @@ -0,0 +1,201 @@ +############################################################################# +# Make file for AutoMix sampler +# Written by: David Hastie, University of Bristol +# Last edited: 14/11/04 +# +# If you use the AutoMix sampler please reference the author and work +# as instructed on the website http://www.davidhastie.me.uk/AutoMix +# Initially this reference 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. +# +# All code has only been tested using the GNU compiler gcc. +############################################################################# + +# The file contains commands for compiling programs using optimisation and +# no debugging. Also included are commands for compiling with debugging. +# The executables and object files for the debugging are identified by +# the final letter "d". For example, amtoy1.exe is the AutoMix sampler for the +# toy problem (see Ph.D. thesis section 5.5.1 for details) and amtoy1d.exe is +# the version that allows debugging using a debugger such as gdb +# Flags and compiler names should be changed as necessary. + +# C compiler +CC=gcc + +# Normal optimised flags +CFLAGS=-O3 -o + +# Debugging flags +CFLAGSD=-g -o + +# Dependency flags +DEPFLAGS=-c + +# Dependency flags with debugging +DEPFLAGSD=-g -c + +# Libraries +LIB=-lm + +# Compiler flags: +# -g -- Enable debugging through gdb +# -O3 -- Optimise progs +# -Wall -- Turn on all warnings +# -lm -- Use the maths library +# -c -- Compile flag for dependencies + + +###### Type "make all" to make all files in the folder ##### + +all: amtoy1 amtoy2 amcpt amcptrs amrb9 amddi amtoy1d amtoy2d amcptd amcptrsd amrb9d amddid + +###### Normal (already debugged) progs ############ + +# Toy example 1 +amtoy1: automix.c usertoy1.o sd.o gammafns.o sokal.o + $(CC) $(CFLAGS) amtoy1 automix.c usertoy1.o sd.o gammafns.o sokal.o $(LIB) + +# Toy example 2 +amtoy2: automix.c usertoy2.o sd.o gammafns.o sokal.o + $(CC) $(CFLAGS) amtoy2 automix.c usertoy2.o sd.o gammafns.o sokal.o $(LIB) + +# Change point problem +amcpt: automix.c usercpt.o gammafns.o sd.o sokal.o + $(CC) $(CFLAGS) amcpt automix.c usercpt.o sd.o gammafns.o sokal.o $(LIB) + +# Rescaled change point problem +amcptrs: automix.c usercptrs.o gammafns.o sd.o sokal.o + $(CC) $(CFLAGS) amcptrs automix.c usercptrs.o sd.o gammafns.o sokal.o $(LIB) + +# Rb9 problem +amrb9: automix.c userrb9.o gammafns.o sd.o sokal.o + $(CC) $(CFLAGS) amrb9 automix.c userrb9.o sd.o gammafns.o sokal.o $(LIB) + +# DDI Clinical trial problem +amddi: automix.c ddidata.h userddi.o gammafns.o sd.o sokal.o + $(CC) $(CFLAGS) amddi automix.c userddi.o sd.o gammafns.o sokal.o $(LIB) + +### AutoMix dependencies (already debugged) + +# Calculates loggamma function +gammafns.o: gammafns.c + $(CC) $(DEPFLAGS) gammafns.c + +# Random number generator +sd.o: sd.c + $(CC) $(DEPFLAGS) sd.c -DDOUB -DRETS + +# Calculates autocorrelation using Sokal's method (Green and Han,92) +sokal.o: sokal.c + $(CC) $(DEPFLAGS) sokal.c + +### User supplied functions (already debugged) + +# Toy example 1 +usertoy1.o: usertoy1.c + $(CC) $(DEPFLAGS) usertoy1.c + +# Toy example 2 +usertoy2.o: usertoy2.c + $(CC) $(DEPFLAGS) usertoy2.c + +# Change point problem +usercpt.o: usercpt.c + $(CC) $(DEPFLAGS) usercpt.c + +# Rescaled change point problem +usercptrs.o: usercptrs.c + $(CC) $(DEPFLAGS) usercptrs.c + +# Rb9 problem +userrb9.o: userrb9.c + $(CC) $(DEPFLAGS) userrb9.c + +# DDI clinical trial problem +userddi.o: userddi.c + $(CC) $(DEPFLAGS) userddi.c + +###### Progs to be debugged ############ + +# Toy example 1 +amtoy1d: automix.c usertoy1d.o sdd.o gammafnsd.o sokald.o + $(CC) $(CFLAGSD) amtoy1d automix.c usertoy1d.o sdd.o gammafnsd.o sokald.o $(LIB) + +# Toy example 2 +amtoy2d: automix.c usertoy2d.o sdd.o gammafnsd.o sokald.o + $(CC) $(CFLAGSD) amtoy2d automix.c usertoy2d.o sdd.o gammafnsd.o sokald.o $(LIB) + +# Change point problem +amcptd: automix.c usercptd.o gammafnsd.o sdd.o sokald.o + $(CC) $(CFLAGSD) amcptd automix.c usercptd.o sdd.o gammafnsd.o sokald.o $(LIB) + +# Rescaled change point problem +amcptrsd: automix.c usercptrsd.o gammafnsd.o sdd.o sokald.o + $(CC) $(CFLAGSD) amcptrsd automix.c usercptrsd.o sdd.o gammafnsd.o sokald.o $(LIB) + +# Rb9 problem +amrb9d: automix.c userrb9d.o gammafnsd.o sdd.o sokald.o + $(CC) $(CFLAGSD) amrb9d automix.c userrb9d.o sdd.o gammafnsd.o sokald.o $(LIB) + +# DDI Clinical trial problem +amddid: automix.c ddidata.h userddid.o gammafnsd.o sdd.o sokald.o + $(CC) $(CFLAGSD) amddid automix.c userddid.o sdd.o gammafnsd.o sokald.o $(LIB) + +# (Old) Toy problems compiled with automix2.c program, implementing +# adaptation through regeneration (automix2.c not included in distribution) + +### AutoMix dependencies (to be debugged) + +# Calculates loggamma function +gammafnsd.o: gammafns.c + cp gammafns.c gammafnsd.c + $(CC) $(DEPFLAGSD) gammafnsd.c + +# Random number generator +sdd.o: sd.c + cp sd.c sdd.c + $(CC) $(DEPFLAGSD) sdd.c -DDOUB -DRETS + +# Calculates autocorrelation using Sokal's method (Green and Han,92) +sokald.o: sokal.c + cp sokal.c sokald.c + $(CC) $(DEPFLAGSD) sokald.c + +### User supplied functions (to be debugged) + +# Toy example 1 +usertoy1d.o: usertoy1.c + cp usertoy1.c usertoy1d.c + $(CC) $(DEPFLAGSD) usertoy1d.c + +# Toy example 2 +usertoy2d.o: usertoy2.c + cp usertoy2.c usertoy2d.c + $(CC) $(DEPFLAGSD) usertoy2d.c + +# Change point problem +usercptd.o: usercpt.c + cp usercpt.c usercptd.c + $(CC) $(DEPFLAGSD) usercptd.c + +# Rescaled change point problem +usercptrsd.o: usercptrs.c + cp usercptrs.c usercptrsd.c + $(CC) $(DEPFLAGSD) usercptrsd.c + +# Rb9 problem +userrb9d.o: userrb9.c + cp userrb9.c userrb9d.c + $(CC) $(DEPFLAGSD) userrb9d.c + +# DDI clinical trial problem +userddid.o: userddi.c + cp userddi.c userddid.c + $(CC) $(DEPFLAGSD) userddid.c + +###### Type "make clean" to remove all executables and object files #### + +clean: + rm *.o + rm am* diff --git a/README.txt b/README.txt new file mode 100644 index 0000000..922592e --- /dev/null +++ b/README.txt @@ -0,0 +1,333 @@ +########################################################################## +# README file - The AutoMix sampler +# +# Last edited 26/09/06. +# Sampler developed by David Hastie, 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 users must +# reference the sampler as instructed below. For commercial +# use permission must be sought from the author. To seek permission +# for such use please send an e-mail to d_hastie@hotmail.com +# 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 +# The sampler may be modified and redistributed as desired but the author +# encourages users to register at the above site so that notice can be +# received of updates of the software. +# +# Before use, please read this README file bundled with this software. +# +# 26/09/06 UPDATE +# Bug fixed in change point user files. Comment had been added +# without closing the comment, meaning function loggamma was not defined. +############################################################################ + +##### INSTALLATION of the AutoMix sampler ##### + +The AutoMix package is a C program for Unix-like systems, implementing +the automatic reversible jump MCMC sampler of the same name +described in Chapters 4, 5, and 6 of David Hastie's Ph.D. thesis + +The package distributed consists of: + main C source file automix.c + auxiliary routines in C gammafns.c sd.c sokal.c + makefile Makefile + example user files usertoy1.c usertoy2.c usercpt.c + usercptrs.c userrb9.c userddi.c + header file ddi example ddidata.h + readme file README.txt + +The package is provided as both a gzipped tarfile and as a compressed +windows archive. For a free Unix-like shell that runs under Windows, +we suggest Cygwin (http:\\www.cygwin.com). + +The program has been run successfully: + using GNU compilers and Cygwin under Windows on Intel PCs + using GNU compilers on Sun Workstations + +We believe that with the appropriate changes to the Makefile the +program will run as required on other platforms. + +##### USING the AutoMix sampler ##### + +The AutoMix package is a C program for Unix-like systems, implementing +the automatic reversible jump MCMC sampler of the same name +described in Chapters 4, 5, and 6 of David Hastie's Ph.D. thesis +These notes assume familiarity with this material. + +In particular, potential users should carefully understand the +limitations of using the AutoMix sampler. The reliability of results +from this sampler depends on many factors, including the scaling of +the parameters and the degree of multimodality of the within-model +conditionals of the target distribution. + +To run the sampler for a particular problem the program must be compiled +with a user-provided file containing four C functions that define the +problem in question. Examples of such files for problems considered +within the aforementioned thesis are bundled with the archived AutoMix +software. More details of the functions are provided below. + +The output of the sampler is in the form of several text files +summarising the MCMC sampler. We discuss these files below. +Subsequent analysis may be performed on the output with the use of +statistical packages. We recommend R (http://www.stats.bris.ac.uk/R/) +as a good free choice of such a package. + +##### WRITING USER FUNCTIONS for the AutoMix sampler ##### + +For each example to which the user wishes to apply the AutoMix sampler +the user must supply a file (which is linked at compile time, see Makefile +for examples) containing four user functions written in C. A familiarity +with the C programming language is assumed. + +These functions must have names and arguments detailed below and return +the following information: + +1. A function to get the number of models: + + void getkmax(int *kmax) + + On exit from the function the (integer) number of models under + consideration should be contained in *kmax. + +2. A function to get the dimension of each model: + + void getnk(int kmax,int *nk){ + + Given the number of models kmax, on exit from the function the + vector nk should contain the dimensions of the models. In particular + the dimension of model k should be in n[k-1], for k=1,...,kmax (Note + that although models run from 1,...,kmax, the vectors run from + the element 0 in C ) + +3. A function to get initial conditions (possibly random) for + the RWM for a given model: + + void getic(int k, int nkk, double *rwm){ + + Given the model index k, and the dimension nkk of that model, + on exit from the function the vector rwm should contain the initial + conditions. In particular, rwm[j-1] should contain + the (possibly random) intial state for component j of the parameter + vector associated with model k, (j=1,...,nkk). + + If random initial states are used, the user must also declare the + random number functions in the file (see the example files). + +4. A function to return the log of the target function pi evaluated at + a given point in the state space, up to an additive constant. + + double lpost(int k,int nkk,double *theta,double *llh1){ + + Given the model index k, and parameter vector theta (of dimension nkk), + the function must return the log of the target function (up to + an additive constant) evaluated at this point. If pi is a posterior + distribution, the double *llh1 should contain the likelihood evaluated + at this point (although this is only necessary for returning the + likelihood to output file, and can contain any other value if preferred). + +Any other functions used by these four functions should also be declared +in the user file. They should also be defined if they are not defined in +other files that are linked at compile time. + +The examples provided, with comments, show typical examples of these user +files for the problems under consideration. + +##### COMPILING the AutoMix sampler ##### + +To compile the AutoMix sampler for the examples included within +the AutoMix package, the user should edit the Makefile (supplied) so that +the appropriate C compiler is used (currently set to be the GNU sampler +gcc). Other aspects of the Makefile may be edited if required (for example +the compile flags, or libraries) + +Typing + make all + +in the shell, in the AutoMix folder where the AutoMix distribution was +unzipped to, will compile all the programs that were distributed with this +package. By default, for each example two executable programs are created +the first using optimisation (for faster run times) and the second +using the debugging flags to enable the user to debug the programs. The +supplied programs have all been debugged but we acknowledge that use +of a debugger can often help to understand how the program works. + +Any of the programs can also be made individually, with the appropriate +make command (for example "make amtoy1" , see the Makefile for further +examples). + +The executables have form + amNAME (for the optimised version) + amNAMEd (for the version that can be debugged) + +where NAME is the name of the example. + +To remove the executables and object (.o) files, type + make clean + +To compile the sampler for a new example, the Makefile can be edited to +contain the appropriate commands (remembering to include command for +compiling any dependencies) for the new program. + +##### RUNNING the AutoMix sampler ##### + +The sampler is run by typing the name of the executable, followed +by run-time flags separated by spaces. The run-time flags control +a number of options within the sampler. If flags are not supplied +default values are used. The flags can be used in any order. + +The flags can be summarised as follows (I is assumed to be a positive +integer): + + -mD controls the mode of the sampler. D=0 is mixture fitting; + D=1 skips stage 1 and 2 if a file containing the mixture + parameters is supplied; D=2 fits AutoMix version of + AutoRJ sampler (see Green, 2003 - full reference in thesis). + (Default uses D=0). + -nI run the sampler for max(I,nkk*10000,100000) iterations + in the stage 1 RWM for each model k. (Default uses I=100000) + -NI run the sampler for I Reversible jump iterations in stage 3. + (Default uses I=100000). + -sI initialises the random number generator with seed I. + (Default uses clock as seed). + -aA controls whether or not adaptation is done in stage 3 RJ. If + A=0 no adaptation is done, if A=1 adaptation is done. (Default + has A=1). + -pP controls whether or not random permutation is done in stage 3 + RJ. If P=0 no permutation is done, if P=1 permutation is done. + (Default has P=0). + -tI Controls whether standard Normal or t distributed variables are + used in RWM and in RJ moves. If I=0 Normal variables are used, + otherwise t-distributed variables with I degrees of freedom are + used. (Default I=0). + -fF Uses the string F as the bases for filenames (e.g. if F=output, + filenames are output_log.data, output_mix.data etc). (Default + is F=output) + +As an example, typing + amtoy1 -m0 -N1000000 -p1 -ftoy1 + +runs the optimised mixture fitting version of the toy1 problem +(see thesis, section 5.5.1) with 1 million RJ sweeps, enabling +permutation and storing the output in files of the type toy1_***.data + +Running the sampler produces a summary of how the run is progressing. + +For each of the models: +In stage 1 a countdown of the number of iterations remaining is +printed to screen; +In stage 2 a summary of the mixture fitting is printed to screen. +This summary consists of a countdown of the number of components in +the current mixture, with the iteration number that the last component +was removed and an indicator n if the component was annihilated +naturally, and f if the annihilation was forced. + +In the RJ stage 3 a countdown of the number of iterations remaining +is printed to screen. + +No summary statistics are printed to screen. Instead all output from +the sampler is written to files. + +##### OUTPUT from the AutoMix sampler ##### + +The following files are returned by the AutoMix sampler (assuming the +filestem is "output") + + output_ac.data + + A file containing the estimated autocorrelation coefficient for the + chain. The coefficients go up to the order that was used to + compute the autocorrelation time using Sokal's method + (see Green and Han, 1992) + + + output_adapt.data + + A file summarising the AAP adaptation process for each model - adapting + the scale parameter of the single component RWM. For a model with nk + components, there are 2*nk columns, with the odd columns showing the + evolution of the scale parameters and the even columns showing the + evolution of the average RWM acceptance rate for the corresponding + scale parameter. + + + output_cf.data + + A file summarising the evolution of the mixture fitting for each + model. In particular, for each model the first column records + the number of mixture components at the current iteration, + the third column summarises the cost function (to be minimised + to find the best mixture) and the final column is an indicator + function which takes the value 1 if a component is annihilated + naturally at that iteration, 2 if a component is removed + by forced annihiliation. (The second column is a term that + contibutes to the cost function, see Figueiredo and Jain, 2002). + + + output_k.data + + A file recording the model index at each sweep of the + RJ sampler (after burn-in). + + + output_log.data + + A log file, containing important run statistics, including + run-time options, run seed, mixture and RWM parameters for each + model, and acceptance rates and autocorrelation statistics for the + RJ stage 3. Model probabilities are also returned in this file. + + + output_lp.data + + A file recording the log posterior (1st column) and log likelihood + (2nd column) at each sweep of the RJ sampler (after burn-in). + + + output_mix.data + + A file containing the fitted mixture parameters for the each model. + The file is for use by the AutoMix program and is not easily readable + by the user. The file contains the number of models and the dimension + of each model. Then for each model in turn, the file records the + adapted random walk scaling parameters, the number of components + of the mixture fitted to that model, and for each component in turn + the weight, the mean and the lower triangle of the B matrix. It is + this file that is required by the AutoMix sampler if it is run in + mode 1. + + + output_pk.data + + A file containing the evolution of the model jumping proposal + parameters for the RJ stage 3. Column k is the probability of + jumping to model k. + + +and for each model k... + + output_thetak.data + + A file recording the parameters conditional on model k at each sweep. + + +##### REFERENCING the AutoMix sampler ##### + +The following reference should be used: + +Hastie, D. I. (2004) Towards Automatic Reversible Jump Markov +Chain Monte Carlo, Ph.D. Thesis, Statistics Group, University +of Bristol. + +It is intended that this reference will change to an academic paper + (currently in preparation) in the not too distant future. diff --git a/automix.c b/automix.c new file mode 100644 index 0000000..b075886 --- /dev/null +++ b/automix.c @@ -0,0 +1,1599 @@ +/* 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 d_hastie@hotmail.com + 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 +#include +#include +#include +#include +#include + +#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[kmax]); + +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=1; + seed=0; + mode=0; + adapt=1; + 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;k12){ + printf("\nInvalid mode entered. Mode must be 0,1,2"); + return -100; + } + else if(mode==1){ + if((check=fscanf(fpmix,"%d",&k1))==EOF){ + printf("\nEnd of file encountered before parameters read:"); + printf("\nContinuing using RWM to estimate parameters"); + mode=0; + goto RWMSTART; + } + if(k1!=kmax){ + printf("\nFile kmax contradicts getkmax function:"); + printf("\nContinuing using RWM to estimate parameters"); + mode=0; + goto RWMSTART; + } + for(k1=0;k11.00001){ + printf("\nComponents weights read do not sum to one for k=%d:",k1); + printf("\nContinuing using RWM to estimate parameters"); + mode=0; + goto RWMSTART; + } + if(sumlambda<1.0||sumlambda>1.0){ + for(l1=0;l1=nburn)&&(fmod((sweep-nburn),((nsweepr-nburn)/10))nburn&&u<0.1){ + if(dof>0){ + rt(Znkk,nkk,dof); + } + else{ + gauss(Znkk,nkk); + } + for(j1=0;j10){ + rt(Z,1,dof); + } + else{ + gauss(Z,1); + } + rwmn[j1]=rwm[j1]+sig[k1][j1]*Z[0]; + lpn=lpost(k1,nkk,rwmn,&llhn); + accept=min(1,exp(max(-30.0,min(0.0,lpn-lp)))); + if(sdrand()0){ + for(l2=0;l20.005){ + /*changed to 0.005 from 0.0 -renormalise else */ + for(j1=0;j10){ + for(l2=0;l21)){ + if(Lkk==1){ + stop=1; + } + else{ + if(fmod(Lkk,5)<0.05){ + printf("\n"); + } + printf("%d(%d-f) ",Lkk,count); + forceann=2; + minlambda=lambda[k1][0]; + ldel=0; + for(l1=1;l1lambda[k1][l1]){ + minlambda=lambda[k1][l1]; + ldel=l1; + } + } + if(ldel<(Lkk-1)){ + for(l1=ldel;l1<(Lkk-1);l1++){ + lambda[k1][l1]=lambda[k1][l1+1]; + for(j1=0;j10){ + for(l2=0;l25000){ + stop=1; + } + costfn=costfnnew; + fprintf(fpcf,"%d %lf %lf %d\n",Lkk,lpn,costfnnew,(natann+forceann)); + fflush(NULL); + } + + for(j1=0;j10){ + rt(Znkk,nkk,dof); + } + else{ + gauss(Znkk,nkk); + } + for(j1=0;j10){ + rt(Z,1,dof); + } + else{ + gauss(Z,1); + } + thetan[j1]=theta[j1]+sig[k][j1]*Z[0]; + lpn=lpost(k,nkk,thetan,&llhn); + if(sdrand()1){ + sum=0.0; + for(l1=0;l10){ + for(l1=0;l10){ + rt(&(work[nkk]),nkkn-nkk,dof); + for(j1=nkk;j10){ + for(j1=nkkn;j11){ + sum=0.0; + for(l1=0;l10){ + for(l1=0;l1nburn){ + for(k1=0;k1nburn){ + (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;k1keep&&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))nburn)&&(fmod(sweep-nburn,(nsweep/10))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[nkk]; + double out; + + for(j1=0;j1 + +#define max(A,B) ((A)>(B)?(A):(B)) +#define min(A,B) ((A)<(B)?(A):(B)) + +extern double sdrand(); + +/* Taken from algama.f (PJG) - converted to C using appropriate machine + constants by DIH 04/11/03 */ + +double loggamma(double X){ + + /* This routine calculates the LOG GAMMA function for a positive real + argument X. Computation is based on an algorithm outlined in + references 1 and 2. The program uses rational functions that + theoretically approximate LOG GAMMA to at least 18 significant + decimal digits. The approximation for X > 12 is from reference + 3, while approximations for X < 12.0 are similar to those in + reference 1, but are unpublished. The accuracy achieved depends + on the arithmetic system, the compiler, the intrinsic functions, + and proper selection of the machine-dependent constants. + + Values taken from float.h and algama.f (for XBIG) + + --------------------------------------------------------------- + + Explanation of machine-dependent constants + + beta - radix for the floating-point representation + maxexp - the smallest positive power of beta that overflows + XBIG - largest argument for which LN(GAMMA(X)) is representable + in the machine, i.e., the solution to the equation + LN(GAMMA(XBIG)) = beta**maxexp + XINF - largest machine representable floating-point number; + approximately beta**maxexp. + EPS - The smallest positive floating-point number such that + 1.0+EPS > 1.0 + FRTBIG - Rough estimate of the fourth root of XBIG + + --------------------------------------------------------------- + + Error returns + + The program returns the value XINF for X <= 0.0 or when + overflow would occur. The computation is believed to + be free of underflow and overflow. + + --------------------------------------------------------------- + References: + + 1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for + the Natural Logarithm of the Gamma Function,' Math. Comp. 21, + 1967, pp. 198-203. + + 2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, + 1969. + + 3) Hart, Et. Al., Computer Approximations, Wiley and sons, New + York, 1968. + + ----------------------------------------------------------------- + Start of code + -----------------------------------------------------------------*/ + + int I; + double CORR,D1,D2,D4,EPS,FRTBIG,FOUR,HALF,ONE,PNT68; + double RES,SQRTPI,THRHAL,TWELVE,TWO,XBIG,XDEN,XINF; + double XM1,XM2,XM4,XNUM,Y,YSQ,ZERO; + double C[7],P1[8],P2[8],P4[8],Q1[8],Q2[8],Q4[8]; + + /*--------------------------------------------------------------------- + Mathematical constants + ---------------------------------------------------------------------*/ + + ONE=1.0E0; + HALF=0.5E0; + TWELVE=12.0E0; + ZERO=0.0E0; + FOUR=4.0E0; + THRHAL=1.5E0; + TWO=2.0E0; + PNT68=0.6796875E0; + SQRTPI=0.9189385332046727417803297E0; /* eh? */ + + /*--------------------------------------------------------------------- + Machine dependent parameters + -------------------------------------------------------------------*/ + + XBIG=2.55E305; + XINF=1.79E308; + EPS=2.22E-16; + FRTBIG=2.25E76; + + /*-------------------------------------------------------------------- + Numerator and denominator coefficients for rational minimax + approximation over (EPS,1.5). + -------------------------------------------------------------------*/ + D1=-5.772156649015328605195174E-1; + P1[0]=4.945235359296727046734888E0; + P1[1]=2.018112620856775083915565E2; + P1[2]=2.290838373831346393026739E3; + P1[3]=1.131967205903380828685045E4; + P1[4]=2.855724635671635335736389E4; + P1[5]=3.848496228443793359990269E4; + P1[6]=2.637748787624195437963534E4; + P1[7]=7.225813979700288197698961E3; + Q1[0]=6.748212550303777196073036E1; + Q1[1]=1.113332393857199323513008E3; + Q1[2]=7.738757056935398733233834E3; + Q1[3]=2.763987074403340708898585E4; + Q1[4]=5.499310206226157329794414E4; + Q1[5]=6.161122180066002127833352E4; + Q1[6]=3.635127591501940507276287E4; + Q1[7]=8.785536302431013170870835E3; + + /*--------------------------------------------------------------------- + Numerator and denominator coefficients for rational minimax + Approximation over (1.5,4.0). + ------------------------------------------------------------------*/ + + D2=4.227843350984671393993777E-1; + P2[0]=4.974607845568932035012064E0; + P2[1]=5.424138599891070494101986E2; + P2[2]=1.550693864978364947665077E4; + P2[3]=1.847932904445632425417223E5; + P2[4]=1.088204769468828767498470E6; + P2[5]=3.338152967987029735917223E6; + P2[6]=5.106661678927352456275255E6; + P2[7]=3.074109054850539556250927E6; + Q2[0]=1.830328399370592604055942E2; + Q2[1]=7.765049321445005871323047E3; + Q2[2]=1.331903827966074194402448E5; + Q2[3]=1.136705821321969608938755E6; + Q2[4]=5.267964117437946917577538E6; + Q2[5]=1.346701454311101692290052E7; + Q2[6]=1.782736530353274213975932E7; + Q2[7]=9.533095591844353613395747E6; + + /*-------------------------------------------------------------------- + Numerator and denominator coefficients for rational minimax + Approximation over (4.0,12.0). + -------------------------------------------------------------------*/ + + D4=1.791759469228055000094023E0; + P4[0]=1.474502166059939948905062E4; + P4[1]=2.426813369486704502836312E6; + P4[2]=1.214755574045093227939592E8; + P4[3]=2.663432449630976949898078E9; + P4[4]=2.940378956634553899906876E10; + P4[5]=1.702665737765398868392998E11; + P4[6]=4.926125793377430887588120E11; + P4[7]=5.606251856223951465078242E11; + Q4[0]=2.690530175870899333379843E3; + Q4[1]=6.393885654300092398984238E5; + Q4[2]=4.135599930241388052042842E7; + Q4[3]=1.120872109616147941376570E9; + Q4[4]=1.488613728678813811542398E10; + Q4[5]=1.016803586272438228077304E11; + Q4[6]=3.417476345507377132798597E11; + Q4[7]=4.463158187419713286462081E11; + + /*--------------------------------------------------------------------- + Coefficients for minimax approximation over (12, INF). + -------------------------------------------------------------------*/ + C[0]=-1.910444077728E-03; + C[1]=8.4171387781295E-04; + C[2]=-5.952379913043012E-04; + C[3]=7.93650793500350248E-04; + C[4]=-2.777777777777681622553E-03; + C[5]=8.333333333333333331554247E-02; + C[6]=5.7083835261E-03; + + /*---------------------------------------------------------------------- + 0 < X <= EPS + --------------------------------------------------------------------*/ + Y=X; + if((Y>0)&&(Y<=XBIG)){ + if(Y<=EPS){ + RES=-log(Y); + } + else if(Y<=THRHAL){ + /*----------------------------------------------------------------------- + EPS < X <= 1.5 + ---------------------------------------------------------------------*/ + if(Y=PNT68)){ + XDEN=ONE; + XNUM=ZERO; + for(I=0;I<8;I++){ + XNUM=XNUM*XM1+P1[I]; + XDEN=XDEN*XM1+Q1[I]; + } + RES=CORR+(XM1*(D1+XM1*(XNUM/XDEN))); + } + else{ + XM2=(Y-HALF)-HALF; /*Is XM2 symbol used to agree with other 2 symbols*/ + XDEN=ONE; + XNUM=ZERO; + for(I=0;I<8;I++){ + XNUM=XNUM*XM2+P2[I]; + XDEN=XDEN*XM2+Q2[I]; + } + RES=CORR+XM2*(D2+XM2*(XNUM/XDEN)); + } + } + else if(Y<=FOUR){ + + /*--------------------------------------------------------------------- + 1.5 < X <= 4.0 + -------------------------------------------------------------------*/ + XM2=Y-TWO; + XDEN=ONE; + XNUM=ZERO; + for(I=0;I<8;I++){ + XNUM=XNUM*XM2+P2[I]; + XDEN=XDEN*XM2+Q2[I]; + } + RES=XM2*(D2+XM2*(XNUM/XDEN)); + } + else if(Y<=TWELVE){ + + /*---------------------------------------------------------------------- + 4.0 < X <= 12.0 + --------------------------------------------------------------------*/ + XM4=Y-FOUR; + XDEN=-ONE; + XNUM=ZERO; + for(I=0;I<8;I++){ + XNUM=XNUM*XM4+P4[I]; + XDEN=XDEN*XM4+Q4[I]; + } + RES=D4+XM4*(XNUM/XDEN); + } + else{ + + /*--------------------------------------------------------------------- + X > 12.0, + -------------------------------------------------------------------*/ + RES=ZERO; + if(Y<=FRTBIG){ + RES=C[6]; + YSQ=Y*Y; + for(I=0;I<6;I++){ + RES=RES/YSQ+C[I]; + } + } + RES=RES/Y; + CORR=log(Y); + RES=RES+SQRTPI-HALF*CORR; + RES=RES+Y*(CORR-ONE); + } + } + else{ + + /*---------------------------------------------------------------------- + Return for bad arguments + --------------------------------------------------------------------*/ + RES=XINF; + } + + /*---------------------------------------------------------------------- + Final return + --------------------------------------------------------------------*/ + return(RES); + +} + + +/* Function for generating a Gamma(s,1) random variable + Converted to C from PJG rgamma FORTRAN function + Note: (1/t)*GAMMA(s,1) is GAMMA(s,t) */ + +double rgamma(double s){ + + double exp1,b,c1,c2,c3,c4,c5,u1,u2,w,bu,out; + + exp1=exp(1.0); + + if(s<1){ + b=(s+exp1)/exp1; + c1 = 1.0/s; + LAB_1: + bu=b*sdrand(); + if(bu<=1.0){ + out=exp(max(-30.0,c1*log(bu))); + if(sdrand()>=exp(-out)){ + goto LAB_1; + } + } + else{ + out=-log((b-bu)/s); + if(sdrand()>=pow(out,(s-1.0))){ + goto LAB_1; + } + } + } + else if(s==1.0){ + out=-log(sdrand()); + } + else{ + c1=s-1.0; + c2=(s-1.0/(6.0*s))/c1; + c3=2.0/c1; + c4=c3+2.0; + c5=1.0/sqrt(s); + LAB_2: + u1=sdrand(); + u2=sdrand(); + if(s>2.5){ + u1=u2+c5*(1.0-1.86*u1); + } + if(u1<=0.0||u1>=1.0){ + goto LAB_2; + } + w=c2*u2/u1; + if((c3*u1+w+1.0/w)<=c4){ + goto LAB_3; + } + if((c3*log(u1)-log(w)+w)>=1.0){ + goto LAB_2; + } + LAB_3: + out=c1*w; + } + + return out; + +} diff --git a/register.php b/register.php new file mode 100644 index 0000000..44ba532 --- /dev/null +++ b/register.php @@ -0,0 +1,81 @@ +\nX-Mailer: PHP/" . phpversion()) or print("Failure E-mail not sent"); + + + } + else { + $register_ok=1; + $userid = mysql_insert_id(); + // Let's mail the user! + $subject = "AutoMix"; + $message = "Dear $first_name $last_name, + +This is an automated response, please do not reply. + +Thank you for downloading the AutoMix software. + +Your name and email address will be stored in our user database to ensure that we can +send you notification of any significant upgrades to the AutoMix software. + +Kind regards +David Hastie"; + + + mail($email_address, $subject, $message, "From: David Hastie\nX-Mailer: PHP/" . phpversion()) or print("User notification e-mail not sent. Please e-mail d_hastie@hotmail.com with this message."); + + $email_address2="d_hastie@hotmail.com"; + $subject1="AutoMix: $first_name $lastname has registered to use the software"; + $message1="New user."; + + mail($email_address2, $subject1, $message1, "From: David Hastie\nX-Mailer: PHP/" . phpversion()) or print("Author notification e-mail not sent. Please e-mail d_hastie@hotmail.com with this message."); + } + } + +if($register_ok==1){ + include '../php/automix_success.php'; +} +else{ + include '../php/automix_failure.php'; +} +?> \ No newline at end of file diff --git a/sd.c b/sd.c new file mode 100644 index 0000000..8725146 --- /dev/null +++ b/sd.c @@ -0,0 +1,157 @@ +/* Combined congruential and Tauseworthe generators from SuperDuper + * package. Should work on machines with unsigned long of at least 32 + * bits. JC and JT must be initialized to values with 0 < JC < 2^32 and + * 0 < JT < 2^32. JC must be odd. + * References: Marsaglia, Ananthanarayanan & Paul, 1973, + * Learmonth & Lewis, 1973, and Dudewicz, 1976) + */ + +/* Compilation flags: + -DSUBR sdrand is used as a subroutine sdrand(u) + instead of a function u = sdrand() + -DSUNF generate Sun Fortran compatible entry names + (with trailing _ symbol) + -DDOUB sdrand or its argument are double precision + -DBSD solves problem with Sun Fortran, pre-Solaris versions + (i.e. BSD), with single precision function option. + Do not use this option for Solaris. + -DRETS returns effective seed in argument to sdrni; + if used, it is essential to use a variable, not + a constant, in calling sdrni. + -DLOG prints seed value in file "rnilog". + + Examples: + cc -c -o sd.o -DSUNF -DRETS sd.c + (single precision Sun Fortran function, returning seed value from sdrni) + + cc -c -o sdc.o sd.c + (single precision C function) +*/ + +#include + +#ifdef SUNF +#ifndef DOUB +#include +#endif +#endif + +static unsigned long JC, JT; +static double Norm=4.656612873E-10; + +#ifdef SUNF +#ifdef SUBR + void sdrand_(u) +#else +#ifdef DOUB + double sdrand_() +#else +#ifdef BSD + FLOATFUNCTIONTYPE sdrand_() +#else + float sdrand_() +#endif +#endif +#endif +#else +#ifdef SUBR + void sdrand(u) +#else +#ifdef DOUB + double sdrand() +#else + float sdrand() +#endif +#endif +#endif + +#ifdef SUBR +#ifdef DOUB + double *u; +#else + float *u; +#endif +#endif + +{ + JC = (JC * 69069) & 037777777777; /* congruential part */ + JT ^= JT >> 15; /* tausworthe part */ + JT ^= (JT << 17) & 037777777777; + +#ifdef SUBR +#ifdef DOUB + *u = ((JT ^ JC) >> 1) * Norm; +#else + *u = (float)(((JT ^ JC) >> 1) * Norm); +#endif +#else +#ifdef DOUB + return(((JT ^ JC) >> 1) * Norm); +#else +#ifdef SUNF +#ifdef BSD + RETURNFLOAT((float)((JT ^ JC) >> 1) * Norm); +#else + return((float)((JT ^ JC) >> 1) * Norm); +#endif +#else + return((float)((JT ^ JC) >> 1) * Norm); +#endif +#endif +#endif + +} + +#ifdef SUNF + void sdrni_(i) +#else + void sdrni(i) +#endif + +unsigned long *i; +{ +#ifdef LOG + FILE *stream; +#endif + unsigned long k=*i; + if(k==0) k=time(0); + JT = k/65536; JC = k-65536*JT; + JT = 65536*JT+1; JC = 32768*JC+1; +#ifdef LOG + stream = fopen("rnilog","a+"); + fprintf(stream,"%12d\n",k); + fclose(stream); +#endif +#ifdef RETS + *i = k; +#endif +} + +#ifdef SUNF + void sdpseed_() +#else + void sdpseed() +#endif +{ +printf("%d %d \n",JC,JT); +} + +#ifdef SUNF + void sdset_(i,j) +#else + void sdset(i,j) +#endif +unsigned long *i,*j; +{ +JC = *i; JT = *j; +} + +#ifdef SUNF + void sdget_(i,j) +#else + void sdget(i,j) +#endif +unsigned long *i,*j; +{ +*i = JC; *j = JT; +} diff --git a/sokal.c b/sokal.c new file mode 100644 index 0000000..468a89d --- /dev/null +++ b/sokal.c @@ -0,0 +1,305 @@ +/* Functions to estimates integrated autocorrelation time using + method of Sokal. Taken from PJG function sokal.f + Note that the definition is the sum from minus infinity to + infinity of the autocorrelation function, hence twice Sokal's + definition. */ + +#include +#include + +#define pi 3.141592653589793 + +void fastfr(int nin, double *xreal, double *ximag); + +void sokal(int n, double *xreal, double *var, double *tau, int *m){ + + int i1; + double ximag[n],c,sum; + + if(n>pow(2.0,20.0)){ + printf("\nAuto-correlation length exceeded"); + return; + } + + for(i1=0;i1 PJG 20 March 1987 -> DIH, March 2004 + + xreal = array which on input contains real part of data + for transformation and on output gives real part of the + result,type real, dimension iabs(isize) or greater. + ximag = array which on input contains imaginary part + for transformation and on output gives imaginary part of + the result, type real, dimension iabs(isize) or greater. + isize = integer variable or constant specifying length and + type of transform. The length is iabs(isize) which must be + a power of 2, minimum length 4, maximum 2**20. Illegal length + leads to warning and return. If isize is positive the forward + transform is calculated. If negative the inverse transform is found. + + The transform is defined by + out(r) = sum(j = 0 to n-1) in(j)*exp(-2*pi*i*r*j/n) + if isize = n > 0, + out(r) = sum(j = 0 to n-1) in(j)*exp(2*pi*i*r*j/n) + if isize = -n < 0, + for r = 0,1,2,...,(n-1), + where i = sqrt(-1), and both in(j) and out(j) + are stored in xreal(j+1)+i*ximag(j+1) */ + + + double z,bcos,bsin,temp,cw1,cw2,cw3,sw1,sw2,sw3; + double xs0,xs1,xs2,xs3,ys0,ys1,ys2,ys3,x1,x2,x3,y1,y2,y3; + + int i0,i1,i2,i3,ul[20],n,counta,countb,time,indic,test; + int j0,j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14,j15,j16,j17,j18,j19; + + n = nin<0 ? -nin : nin; + + if(n<4){ + printf("\nFast Fourier length too short"); + return; + } + + /* if this is to be an inverse transform, conjugate the data */ + + if(nin<0){ + for(i1=0;i11){ + if(test&1){ + printf("\nFast Fourier length must be power of 2"); + return; + } + test>>=1; + } + + counta=n/4; + time=0; + + while(counta>0){ + countb=counta*4; + time+=2; + + /* do the transforms required by this stage */ + + z=pi/countb; + bcos=-2.0*pow(sin(z),2.0); + bsin=sin(2.0*z); + + cw1=1.0; + sw1=0.0; + + /* this is the main calculation of radix 4 transforms */ + + for(j1=0;j11){ + /* set up the transform split for the next stage */ + counta/=4; + indic=1; + } + else{ + counta=0; + } + } + if(indic){ + + /* this is the calculation of a radix two stage */ + for(j1=0;j120){ + printf("\nFast Fourier length too long"); + return; + } + + i1=20-time; + for(j1=0;j1 +#include + +#define tpi 6.283185307179586477 + +/* Function loggamma is in file gammafns.c bundled with the AutoMix software */ + +extern double loggamma(double x); + +/* Hyperparameters */ + +double alpha=1.0,beta=200.0,lambda=3.0; + +/* Function to return number of models */ +void getkmax(int *kmax){ + *kmax=6; + return; +} + +/* Function to return the dimension of each model */ +void getnk(int kmax,int *nk){ + int k; + for(k=0;ktop){ + nj=i-nsofar; + nsofar=i; + llh+=(nj*log(h[j])-h[j]*ds[j]); + j++; + if(j>nsteps){ + return lp; + } + top=s[j+1]; + } + } + nj=191-nsofar; + llh+=nj*log(h[j])-h[j]*ds[j]; + lp+=llh; + *llh1=llh; + return lp; + +} diff --git a/usercptrs.c b/usercptrs.c new file mode 100644 index 0000000..998b193 --- /dev/null +++ b/usercptrs.c @@ -0,0 +1,134 @@ +/* User functions for the rescaled change point example in section 5.5.2 + of Ph.D. thesis. + Original data is for times of coal mining disasters and is taken from + Green, 1995 (see thesis for full reference). */ + +#include +#include + +#define tpi 6.283185307179586477 + +/* Function loggamma is in file gammafns.c bundled with the AutoMix software */ + +extern double loggamma(double x); + +/* Hyperparameters */ +/* beta is now 200/1459 (this means mean is 1459 times orig and var is 1459^2 + times orig - in line with transformation of times*/ +double alpha=1.0,beta=0.137,lambda=3.0; + +/* Function to return number of models */ +void getkmax(int *kmax){ + *kmax=6; + return; +} + +/* Function to return the dimension of each model */ +void getnk(int kmax,int *nk){ + int k; + for(k=0;ktop){ + nj=i-nsofar; + nsofar=i; + llh+=(nj*log(h[j])-h[j]*ds[j]); + j++; + if(j>nsteps){ + return lp; + } + top=s[j+1]; + } + } + nj=191-nsofar; + llh+=nj*log(h[j])-h[j]*ds[j]; + lp+=llh; + *llh1=llh; + return lp; + +} diff --git a/userddi.c b/userddi.c new file mode 100644 index 0000000..07ea334 --- /dev/null +++ b/userddi.c @@ -0,0 +1,809 @@ +/* User functions for the AIDS clinical trial example in section 5.5.4 of + Ph.D. thesis. + Example is taken from Han & Carlin, 2001 (see thesis for full reference). */ + +#include +#include + +/* Data provided in user library - bundled with AutoMix software */ +#include "ddidata.h" + +#define tpi 6.283185307179586477 +#define pi 3.141592653589793238 +#define max(A,B) ((A)>(B)?(A):(B)) + +extern double loggamma(double x); + +extern double sdrand(); + +/* hyperparameters */ +double a=3.0,b=0.005; +int rho=24; +double c0[9]={10.0,0.0,0.0,0.0,0.0,0.0,-3.0,0.0,0.0}; +double c1[6]={10.0,0.0,0.0,0.0,-3.0,0.0}; +double D0min1[9]={0.25,1.0,1.0,100.0,1.0,1.0,1.0,1.0,1.0}; +double D1min1[6]={0.25,1.0,100.0,1.0,1.0,1.0}; +double R0[3][3]={{4.0,0.0,0.0},{0.0,1/16.0,0.0},{0.0,0.0,1.0/16.0}}; +double R1[2][2]={{4.0,0.0},{0.0,1/16.0}}; + +/* Internal functions used in required user functions */ +double boxm2(); + +void chol2(int nkk, double **B,int *posdef); + +double det2(int nkk,double **B); + +double lnormprob2(int nkk, double **B, double *mu, double *datai); + +/* Function to return number of models */ +void getkmax(int *kmax){ + *kmax=2; + return; +} + +/* Function to return the dimension of each model */ +void getnk(int kmax,int *nk){ + + nk[0]=16; + nk[1]=10; + return; + +} + +/* Function to return initial conditions for RWM runs */ + +void getic(int k, int nkk, double *rwm){ + int i,j,j1,j2,posdef; + double u,alpha[9],gamma[6],Vmin1[3][3],Umin1[2][2],V[3][3],U[2][2]; + double sigmasq,tausq; + double **Ctest,**temptest; + + /* Note boxm2() is used to return random initial conditions */ + + BEGIN_IC: + if(k==0){ + + for(j=0;j<9;j++){ + alpha[j]=c0[j]+sqrt(1.0/D0min1[j])*boxm2(); + } + for(j=0;j<3;j++){ + for(j1=0;j1<3;j1++){ + Vmin1[j][j1]=0; + V[j][j1]=0; + } + u=exp(0.01*boxm2()); + Vmin1[j][j]=u; + V[j][j]=1.0/u; + } + sigmasq=100.0*exp(0.1*boxm2()); + for(j=0;j<9;j++){ + rwm[j]=alpha[j]; + } + rwm[9]=Vmin1[0][0]; + rwm[10]=Vmin1[1][0]; + rwm[11]=Vmin1[1][1]; + rwm[12]=Vmin1[2][0]; + rwm[13]=Vmin1[2][1]; + rwm[14]=Vmin1[2][2]; + rwm[15]=sigmasq; + + temptest=(double**)malloc(3*sizeof(double)); + Ctest=(double**)malloc(5*sizeof(double)); + for(j=0;j<5;j++){ + Ctest[j]=(double*)malloc(5*sizeof(double)); + if(j<3){ + temptest[j]=(double*)malloc(5*sizeof(double)); + } + } + + for(i=0;i<467;i++){ + for(j=0;j<3;j++){ + for(j1=0;j1<5;j1++){ + temptest[j][j1]=0.0; + } + for(j1=0;j1 +#include + +#define tpi 6.283185307179586477 + +extern double loggamma(double x); + +extern double sdrand(); + +/* Hyperparameters */ +double alpha1=2.0,alpha2=1.0,beta1=0.1,beta2=2.0; + +/* Internal functions used in required user functions */ +double boxm(); + +/* Function to return number of models */ +void getkmax(int *kmax){ + *kmax=10; + return; +} + +/* Function to return the dimension of each model */ +void getnk(int kmax,int *nk){ + int k; + for(k=0;k<6;k++){ + nk[k]=4; + } + for(k=6;k +#include + +#define tpi 6.283185307179586477 + +extern double sdrand(); + +/* Internal functions used in required user functions */ +double boxm(); + +/* Function to return number of models */ +void getkmax(int *kmax){ + *kmax=2; + return; +} + +/* Function to return the dimension of each model */ +void getnk(int kmax,int *nk){ + int k; + for(k=0;k +#include + +#define tpi 6.283185307179586477 + +/* Function to return number of models */ +void getkmax(int *kmax){ + *kmax=5; + return; +} + +/* Function to return the dimension of each model */ +void getnk(int kmax,int *nk){ + int k; + for(k=0;k