# include # include # include # include # include # include "dream1.h" # include "rnglib.h" /******************************************************************************/ void chain_init ( int chain_num, double fit[], int gen_num, int mnor_ind[], int mnor_num, int par_num, double parm[], int prior_den_par_num[], double prior_den_par_val[], int prior_den_type[], double sample_likelihood ( int par_num, double zp[] ), int uni_ind[], int uni_num, double z[] ) /******************************************************************************/ /* Purpose: CHAIN_INIT starts Markov chains from a prior distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 May 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int CHAIN_NUM, the total number of chains. Output, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of each sample. Input, int GEN_NUM, the total number of generations. Input, int MNOR_IND[MNOR_NUM], indices of parameters associated with multivariate normal distributions. Input, int MNOR_NUM, the number of parameters with a multivariate normal distribution. Input, int PAR_NUM, the total number of parameters. Input, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters for a multivariate normal distribution. Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of variables associated with the prior density for each parameter. Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values associated with prior densities. Input, int PRIOR_DEN_TYPE[PAR_NUM], the type of the prior density. Input, double SAMPLE_LIKELIHOOD ( int PAR_NUM, double ZP[] ), a user-defined function to evaluate the log likelihood function of a sample vector ZP of dimension PAR_NUM. Input, int UNI_IND[UNI_NUM], indices of parameters associated with univariate distributions. Input, int UNI_NUM, the number of parameters with a univariate distribution. Output, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. */ { int c; int p; double *zp; for ( c = 0; c < chain_num; c++ ) { zp = prior_sample ( mnor_ind, mnor_num, par_num, parm, prior_den_par_num, prior_den_par_val, prior_den_type, uni_ind, uni_num ); for ( p = 0; p < par_num; p++ ) { z[p+c*par_num+0*par_num*chain_num] = zp[p]; } fit[c+0*chain_num] = sample_likelihood ( par_num, zp ); free ( zp ); } return; } /******************************************************************************/ void chain_init_print ( int chain_num, double fit[], int gen_num, int par_num, int restart_flag, char *restart_read_filename, double z[] ) /******************************************************************************/ /* Purpose: CHAIN_INIT_PRINT prints the initial values for Markov chains. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 April 2013 Author: John Burkardt Parameters: Input, int CHAIN_NUM, the total number of chains. Input, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of each sample. Input, int GEN_NUM, the total number of generations. Input, int PAR_NUM, the total number of parameters. Input, int RESTART_FLAG, is TRUE if the chains are to be initialized from a restart file. Input, char *RESTART_READ_FILENAME, the name of the restart file. Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. */ { int i; int j; printf ( "\n" ); printf ( "CHAIN_INIT_PRINT\n" ); printf ( " Display initial values of Markov chains.\n" ); if ( restart_flag ) { printf ( " Initialization from restart file \"%s\"\n", restart_read_filename ); } else { printf ( " Initialization by sampling prior density.\n" ); } for ( j = 0; j < chain_num; j++ ) { printf ( "\n" ); printf ( " Chain %d\n", j ); printf ( " Fitness %g\n", fit[j+0*chain_num] ); for ( i = 0; i < par_num; i++ ) { printf ( " %g", z[i+j*par_num+0*par_num*chain_num] ); if ( ( i + 1 ) % 5 == 0 || i == par_num - 1 ) { printf ( "\n" ); } } } return; } /******************************************************************************/ void chain_outliers ( int chain_num, int gen_index, int gen_num, int par_num, double fit[], double z[] ) /******************************************************************************/ /* Purpose: CHAIN_OUTLIERS identifies and modifies outlier chains during burn-in. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Reference: Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, Dave Higdon, Accelerating Markov Chain Monte Carlo Simulation by Differential Evolution with Self-Adaptive Randomized Subspace Sampling, International Journal of Nonlinear Sciences and Numerical Simulation, Volume 10, Number 3, March 2009, pages 271-288. Parameters: Input, int CHAIN_NUM, the total number of chains. Input, int GEN_INDEX, the index of the current generation. 2 <= GEN_INDEX <= GEN_NUM. Input, int GEN_NUM, the total number of generations. Input, int PAR_NUM, the total number of parameters. Input/output, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of each sample. Input/output, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. */ { double *avg; double avg_max; double *avg_sorted; int best; int i; int ind1; int ind3; int j; int klo; int knum; int k; int outlier_num; double q1; double q3; double qr; double t; klo = ( ( gen_index + 1 ) / 2 ) - 1; knum = gen_index + 1 - klo; avg = ( double * ) malloc ( chain_num * sizeof ( double ) ); for ( j = 0; j < chain_num; j++ ) { t = 0.0; for ( k = klo; k <= gen_index; k++ ) { t = t + fit[j+k*chain_num]; } avg[j] = t / ( double ) ( knum ); } /* Set BEST to be the index of the chain with maximum average. */ best = 0; avg_max = avg[0]; for ( j = 1; j < chain_num; j++ ) { if ( avg_max < avg[j] ) { best = j; avg_max = avg[j]; } } /* Determine the indices of the chains having averages 1/4 "above" and "below" the average. */ avg_sorted = r8vec_copy_new ( chain_num, avg ); r8vec_sort_heap_a ( chain_num, avg_sorted ); ind1 = r8_round_i4 ( 0.25 * ( double ) ( chain_num ) ); ind3 = r8_round_i4 ( 0.75 * ( double ) ( chain_num ) ); q1 = avg_sorted[ind1]; q3 = avg_sorted[ind3]; qr = q3 - q1; free ( avg_sorted ); /* Identify outlier chains, and replace their later samples with values from the "best" chain. */ outlier_num = 0; for ( j = 0; j < chain_num; j++ ) { if ( avg[j] < q1 - 2.0 * qr ) { outlier_num = outlier_num + 1; for ( i = 0; i < par_num; i++ ) { z[i+j*par_num+gen_index*par_num*chain_num] = z[i+best*par_num+gen_index*par_num*chain_num]; } for ( k = klo; k <= gen_index; k++ ) { fit[j+k*chain_num] = fit[best+k*chain_num]; } } } /* List the outlier chains. */ if ( 0 < outlier_num ) { printf ( "\n" ); printf ( "CHAIN_OUTLIERS:\n" ); printf ( " At iteration %d found %d outlier chains,\n", gen_index, outlier_num ); printf ( " whose indices appear below, and for which samples\n" ); printf ( " from the chain with the largest log likelihood function,\n" ); printf ( " index number %d will be substituted.\n", best ); for ( j = 0; j < chain_num; j++ ) { if ( avg[j] < q1 - 2.0 * qr ) { printf ( " %d\n", j ); } } } free ( avg ); return; } /******************************************************************************/ void chain_write ( char *chain_filename, int chain_num, double fit[], int gen_num, int par_num, double z[] ) /******************************************************************************/ /* Purpose: CHAIN_WRITE writes samples of each chain to separate files. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version John Burkardt. Parameters: Input, char *CHAIN_FILENAME, the base name for the files. It should include a string of 0's that will be replaced by the indices of each chain. Input, int CHAIN_NUM, the total number of chains. Input, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of each sample. Input, int GEN_NUM, the total number of generations. Input, int PAR_NUM, the total number of parameters. Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. */ { FILE *chain_unit; char chain_filename2[255]; int i; int j; int k; /* Make a temporary copy of the filename template, which we can alter. */ strcpy ( chain_filename2, chain_filename ); /* Write parameter samples of all chains. */ printf ( "\n" ); printf ( "CHAIN_WRITE:\n" ); for ( j = 0; j < chain_num; j++ ) { chain_unit = fopen ( chain_filename2, "wt" ); if ( ! chain_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CHAIN_WRITE - Fatal error!\n" ); fprintf ( stderr, " Could not open file \"%s\".\n", chain_filename2 ); exit ( 1 ); } fprintf ( chain_unit, "DREAM.C:Parameters_and_log_likelihood_for_chain_#%d\n", j ); for ( k = 0; k < gen_num; k++ ) { fprintf ( chain_unit, " %d %g", k, fit[j+k*chain_num] ); for ( i = 0; i < par_num; i++ ) { fprintf ( chain_unit, " %g", z[i+j*par_num+k*par_num*chain_num] ); } fprintf ( chain_unit, "\n" ); } fclose ( chain_unit ); printf ( " Created file \"%s\".\n", chain_filename2 ); filename_inc ( chain_filename2 ); } return; } /******************************************************************************/ void cr_dis_update ( int chain_index, int chain_num, double cr_dis[], int cr_index, int cr_num, int cr_ups[], int gen_index, int gen_num, int par_num, double z[] ) /******************************************************************************/ /* Purpose: CR_DIS_UPDATE updates the CR distance. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2018 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int CHAIN_INDEX, the index of the chain. 0 <= CHAIN_INDEX < CHAIN_NUM. Input, int CHAIN_NUM, the total number of chains. Input/output, double CR_DIS[CR_NUM], the CR distances. Input, int CR_INDEX, the index of the CR. 0 <= CR_INDEX < CR_NUM. Input, int CR_NUM, the total number of CR values. Input/output, int CR_UPS[CR_NUM], the number of updates for each CR. Input, int GEN_INDEX, the index of the generation. 0 <= GEN_INDEX < GEN_NUM. Input, int GEN_NUM, the total number of generations. Input, int PAR_NUM, the total number of parameters. Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. */ { int i; int i1; int i2; double *std; /* Compute the standard deviations. */ std = std_compute ( chain_num, gen_index, gen_num, par_num, z ); /* Increment the update count. */ cr_ups[cr_index] = cr_ups[cr_index] + 1; /* Update the CR distance. */ for ( i = 0; i < par_num; i++ ) { i1 = i + chain_index * par_num + gen_index * par_num * chain_num; i2 = i + chain_index * par_num + ( gen_index - 1 ) * par_num * chain_num; cr_dis[cr_index] = cr_dis[cr_index] + pow ( ( z[i2] - z[i1] ) / std[i], 2 ); } free ( std ); return; } /******************************************************************************/ int cr_index_choose ( int cr_num, double cr_prob[] ) /******************************************************************************/ /* Purpose: CR_INDEX_CHOOSE chooses a CR index. Discussion: Index I is chosen with probability CR_PROB(I). Licensing: This code is distributed under the GNU LGPL license. Modified: 14 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int CR_NUM, the total number of CR values. Input, double CR_PROB[CR_NUM], the probability of each CR. Output, int CR_INDEX_CHOOSE, the index of the CR. 0 <= CR_INDEX_CHOOSE < CR_NUM. */ { int cr_index; int i; int n; int *tmp_index; if ( cr_num == 1 ) { cr_index = 0; } else { n = 1; tmp_index = i4vec_multinomial_sample ( n, cr_prob, cr_num ); for ( i = 0; i < cr_num; i++ ) { if ( tmp_index[i] == 1 ) { cr_index = i; break; } } free ( tmp_index ); } return cr_index; } /******************************************************************************/ void cr_init ( double cr[], double cr_dis[], int cr_num, double cr_prob[], int cr_ups[] ) /******************************************************************************/ /* Purpose: CR_INIT initializes the crossover probability values. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Output, double CR[CR_NUM], the CR values. Output, double CR_DIS[CR_NUM], the CR distances. Input, int CR_NUM, the total number of CR values. Output, double CR_PROB[CR_NUM], the probability of each CR. Output, int CR_UPS[CR_NUM], the number of updates for each CR. */ { int i; for ( i = 0; i < cr_num; i++ ) { cr[i] = ( double ) ( i + 1 ) / ( double ) ( cr_num ); cr_dis[i] = 1.0; cr_prob[i] = 1.0 / ( double ) ( cr_num ); cr_ups[i] = 1; } return; } /******************************************************************************/ void cr_prob_update ( double cr_dis[], int cr_num, double cr_prob[], int cr_ups[] ) /******************************************************************************/ /* Purpose: CR_PROB_UPDATE updates the CR probabilities. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double CR_DIS[CR_NUM], the CR distances. Input, int CR_NUM, the total number of CR values. Output, double CR_PROB[CR_NUM], the updated CR probabilities. Input, int CR_UPS[CR_NUM], the number of updates for each CR. */ { double cr_prob_sum; int i; for ( i = 0; i < cr_num - 1; i++ ) { cr_prob[i] = cr_dis[i] / ( double ) cr_ups[i]; } cr_prob_sum = r8vec_sum ( cr_num, cr_prob ); for ( i = 0; i < cr_num - 1; i++ ) { cr_prob[i] = cr_prob[i] / cr_prob_sum; } return; } /******************************************************************************/ double *diff_compute ( int chain_num, int gen_index, int gen_num, int jump_dim[], int jump_num, int pair_num, int par_num, int r[], double z[] ) /******************************************************************************/ /* Purpose: DIFF_COMPUTE computes the differential evolution. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Reference: Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, Dave Higdon, Accelerating Markov Chain Monte Carlo Simulation by Differential Evolution with Self-Adaptive Randomized Subspace Sampling, International Journal of Nonlinear Sciences and Numerical Simulation, Volume 10, Number 3, March 2009, pages 271-288. Parameters: Input, int CHAIN_NUM, the total number of chains. Input, int GEN_INDEX, the index of the current generation. 1 <= GEN_INDEX <= GEN_NUM. Input, int GEN_NUM, the total number of generations. Input, int JUMP_DIM[JUMP_NUM], the dimensions in which a jump is to be made. Input, int JUMP_NUM, the number of dimensions in which a jump will be made. 0 <= JUMP_NUM <= PAR_NUM. Input, int PAIR_NUM, the number of pairs of crossover chains. Input, int PAR_NUM, the total number of parameters. Input, int R[2*PAIR_NUM], pairs of chains used to compute differences. Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. Output, double DIFF_COMPUTE[PAR_NUM], the vector of pair differences. */ { double *diff; int i1; int i2; int j; int k; int pair; int r1; int r2; /* Produce the difference of the pairs used for population evolution. */ diff = r8vec_zero_new ( par_num ); for ( pair = 0; pair < pair_num; pair++ ) { r1 = r[0+pair*2]; r2 = r[1+pair*2]; for ( j = 0; j < jump_num; j++ ) { k = jump_dim[j]; i1 = k+r1*par_num+(gen_index-1)*par_num*chain_num; i2 = k+r2*par_num+(gen_index-1)*par_num*chain_num; diff[k] = diff[k] + ( z[i1] - z[i2] ); } } return diff; } /******************************************************************************/ void dream_algm ( int chain_num, double cov_mat[], int cr_num, double fit[], int gen_num, double gr[], int *gr_conv, int *gr_count, int gr_num, double gr_threshold, double jumprate_table[], int jumpstep, double limits[], double logdet, int mnor_ind[], double mnor_mean[], int mnor_num, int pair_num, int par_num, int printstep, int prior_den_par_num[], double prior_den_par_val[], int prior_den_type[], double sample_likelihood ( int par_num, double zp[] ), int uni_ind[], int uni_num, double z[] ) /******************************************************************************/ /* Purpose: DREAM_ALGM gets a candidate parameter sample. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 May 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Reference: Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, Dave Higdon, Accelerating Markov Chain Monte Carlo Simulation by Differential Evolution with Self-Adaptive Randomized Subspace Sampling, International Journal of Nonlinear Sciences and Numerical Simulation, Volume 10, Number 3, March 2009, pages 271-288. Parameters: Input, int CHAIN_NUM, the total number of chains. Input, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance matrix for a multivariate normal density. Input, int CR_NUM, the total number of CR values. Input, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of each sample. Input, int GEN_NUM, the total number of generations. Input, double GR[PAR_NUM*GR_NUM], the Gelman-Rubin R statistic. Input/output, int *GR_CONV, the Gelman-Rubin convergence flag. Input/output, int *GR_COUNT, counts the number of generations at which the Gelman-Rubin statistic has been computed. Input, int GR_NUM, the number of times the Gelman-Rubin statistic may be computed. Input, double GR_THRESHOLD, the convergence tolerance for the Gelman-Rubin statistic. Input, double JUMPRATE_TABLE[PAR_NUM], the jumprate table. Input, int JUMPSTEP, forces a "long jump" every JUMPSTEP generations. Input, double LIMITS[2*PAR_NUM], lower and upper bounds for each parameter. Input, double LOGDET, the logarithm of the determinant of the covariance matrix. Input, int MNOR_IND[MNOR_NUM], the indices of parameters with multivariate density. Input, int MNOR_MEAN[MNOR_NUM], the means associated with the multivariate density. Input, int MNOR_NUM, the number of parameters with a multivariate normal distribution. Input, int PAIR_NUM, the number of pairs of crossover chains. Input, int PAR_NUM, the total number of parameters. Input, int PRINTSTEP, the interval between generations on which the Gelman-Rubin statistic will be computed and written to a file. Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of variables associated with the prior density for each parameter. Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values associated with prior densities. Input, int PRIOR_DEN_TYPE[PAR_NUM], the type of the prior density. Input, double SAMPLE_LIKELIHOOD ( int PAR_NUM, double ZP[] ), a user-defined function to evaluate the log likelihood function of a sample vector ZP of dimension PAR_NUM. Input, int UNI_IND[UNI_NUM], indices of parameters associated with univariate distributions. Input, int UNI_NUM, the number of parameters with a univariate distribution. Output, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. Local parameters: Local, int CHAIN_INDEX, the index of the current chain. 1 <= CHAIN_INDEX <= CHAIN_NUM. Local, double CR[CR_NUM], the CR values. Local, double CR_DIS[CR_NUM], the CR distances. Local, int CR_INDEX, the index of the selected CR value. 1 <= CR_INDEX <= CR_NUM. Local, double CR_PROB[CR_NUM], the probability of each CR. Local, double CR_UPS[CR_NUM], the number of updates for each CR. Local, int GEN_INDEX, the index of the current generation. 1 <= GEN_INDEX <= GEN_NUM. Local, double ZP[PAR_NUM], a candidate sample. Local, int ZP_ACCEPT, the number of candidates accepted. Local, double ZP_ACCEPT_RATE, the rate at which generated candidates were accepted. Local, int ZP_COUNT, the number of candidates generated. Local, double ZP_RATIO, the Metropolis ratio for a candidate. */ { int chain_index; double *cr; double *cr_dis; int cr_index; double *cr_prob; int *cr_ups; int gen_index; int i; double pd1; double pd2; double r; double *zp; int zp_accept; double zp_accept_rate; int zp_count; double zp_fit; double *zp_old; double zp_old_fit; double zp_ratio; zp_old = ( double * ) malloc ( par_num * sizeof ( double ) ); zp_count = 0; zp_accept = 0; /* Initialize the CR values. */ cr = ( double * ) malloc ( cr_num * sizeof ( double ) ); cr_dis = ( double * ) malloc ( cr_num * sizeof ( double ) ); cr_prob = ( double * ) malloc ( cr_num * sizeof ( double ) ); cr_ups = ( int * ) malloc ( cr_num * sizeof ( int ) ); cr_init ( cr, cr_dis, cr_num, cr_prob, cr_ups ); for ( gen_index = 1; gen_index < gen_num; gen_index++ ) { for ( chain_index = 0; chain_index < chain_num; chain_index++ ) { /* Choose CR_INDEX, the index of a CR. */ cr_index = cr_index_choose ( cr_num, cr_prob ); /* Generate a sample candidate ZP. */ zp = sample_candidate ( chain_index, chain_num, cr, cr_index, cr_num, gen_index, gen_num, jumprate_table, jumpstep, limits, pair_num, par_num, z ); zp_count = zp_count + 1; /* Compute the log likelihood function for ZP. */ zp_fit = sample_likelihood ( par_num, zp ); for ( i = 0; i < par_num; i++ ) { zp_old[i] = z[i+chain_index*par_num+(gen_index-1)*par_num*chain_num]; } zp_old_fit = fit[chain_index+(gen_index-1)*chain_num]; /* Compute the Metropolis ratio for ZP. */ pd1 = prior_density ( cov_mat, logdet, mnor_ind, mnor_mean, mnor_num, par_num, prior_den_par_num, prior_den_par_val, prior_den_type, uni_ind, uni_num, zp ); pd2 = prior_density ( cov_mat, logdet, mnor_ind, mnor_mean, mnor_num, par_num, prior_den_par_num, prior_den_par_val, prior_den_type, uni_ind, uni_num, z+0+chain_index*par_num+(gen_index-1)*par_num*chain_num ); zp_ratio = exp ( ( zp_fit + log ( pd1 ) ) - ( zp_old_fit + log ( pd2 ) ) ); zp_ratio = r8_min ( zp_ratio, 1.0 ); /* Accept the candidate, or copy the value from the previous generation. */ r = r8_uniform_01_sample ( ); if ( r <= zp_ratio ) { for ( i = 0; i < par_num; i++ ) { z[i+chain_index*par_num+gen_index*par_num*chain_num] = zp[i]; } zp_accept = zp_accept + 1; fit[chain_index+gen_index*chain_num] = zp_fit; } else { for ( i = 0; i < par_num; i++ ) { z[i+chain_index*par_num+gen_index*par_num*chain_num] = zp_old[i]; } fit[chain_index+gen_index*chain_num] = zp_old_fit; } /* Update the CR distance. */ if ( ! gr_conv ) { if ( 1 < cr_num ) { cr_dis_update ( chain_index, chain_num, cr_dis, cr_index, cr_num, cr_ups, gen_index, gen_num, par_num, z ); } } free ( zp ); } /* Update the multinomial distribution of CR. */ if ( ! gr_conv ) { if ( 1 < cr_num ) { if ( ( gen_index + 1 ) % 10 == 0 ) { cr_prob_update ( cr_dis, cr_num, cr_prob, cr_ups ); } } } /* Every PRINTSTEP interval, * compute the Gelman Rubin R statistic for this generation, and determine if convergence has occurred. */ if ( ( gen_index + 1 ) % printstep == 0 ) { gr_compute ( chain_num, gen_index, gen_num, gr, gr_conv, gr_count, gr_num, gr_threshold, par_num, z ); } /* Check for outlier chains. */ if ( ! gr_conv ) { if ( ( gen_index + 1 ) % 10 == 0 ) { chain_outliers ( chain_num, gen_index, gen_num, par_num, fit, z ); } } } /* Compute the acceptance rate. */ zp_accept_rate = ( double ) ( zp_accept ) / ( double ) ( zp_count ); printf ( "\n" ); printf ( " The acceptance rate is %g\n", zp_accept_rate ); free ( cr ); free ( cr_dis ); free ( cr_prob ); free ( cr_ups ); free ( zp_old ); return; } /******************************************************************************/ void filename_inc ( char *filename ) /******************************************************************************/ /* Purpose: FILENAME_INC increments a partially numeric file name. Discussion: It is assumed that the digits in the name, whether scattered or connected, represent a number that is to be increased by 1 on each call. If this number is all 9's on input, the output number is all 0's. Non-numeric letters of the name are unaffected. If the name is empty, then the routine stops. If the name contains no digits, the empty string is returned. Example: Input Output ----- ------ "a7to11.txt" "a7to12.txt" (typical case. Last digit incremented) "a7to99.txt" "a8to00.txt" (last digit incremented, with carry.) "a9to99.txt" "a0to00.txt" (wrap around) "cat.txt" " " (no digits to increment) " " STOP! (error) Licensing: This code is distributed under the GNU LGPL license. Modified: 22 November 2011 Author: John Burkardt Parameters: Input/output, char *FILENAME, the filename to be incremented. */ { int change; int n; char *t; n = strlen ( filename ); if ( n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILENAME_INC - Fatal error!\n" ); fprintf ( stderr, " The input string is empty.\n" ); exit ( 1 ); } change = 0; t = filename + n - 1; while ( 0 < n ) { if ( '0' <= *t && *t <= '9' ) { change = change + 1; if ( *t == '9' ) { *t = '0'; } else { *t = *t + 1; return; } } t--; n--; } /* No digits were found. Return blank. */ if ( change == 0 ) { n = strlen ( filename ); t = filename + n - 1; while ( 0 < n ) { *t = ' '; t--; n--; } } return; } /******************************************************************************/ void gr_compute ( int chain_num, int gen_index, int gen_num, double gr[], int *gr_conv, int *gr_count, int gr_num, double gr_threshold, int par_num, double z[] ) /******************************************************************************/ /* Purpose: GR_COMPUTE computes the Gelman Rubin statistics R used to check convergence. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Reference: Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, Dave Higdon, Accelerating Markov Chain Monte Carlo Simulation by Differential Evolution with Self-Adaptive Randomized Subspace Sampling, International Journal of Nonlinear Sciences and Numerical Simulation, Volume 10, Number 3, March 2009, pages 271-288. Parameters: Input, int CHAIN_NUM, the total number of chains. Input, int GEN_INDEX, the index of the current generation. 0 < GEN_INDEX < GEN_NUM. Input, int GEN_NUM, the total number of generations. Output, double GR[PAR_NUM*GR_NUM], the Gelman-Rubin R statistic. Output, int *GR_CONV, the Gelman-Rubin convergence flag. Input/output, int *GR_COUNT, counts the number of generations at which the Gelman-Rubin statistic has been computed. Input, int GR_NUM, the number of times the Gelman-Rubin statistic may be computed. Input, double GR_THRESHOLD, the convergence tolerance for the Gelman-Rubin statistic. Input, int PAR_NUM, the total number of parameters. Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. */ { double b_var; int chain_index; int ind0; int k; double mean_all; double *mean_chain; int par_index; double rnd0; double s; double s_sum; double var; double w_var; ind0 = ( ( gen_index + 1 ) / 2 ) - 1; rnd0 = ( double ) ( ind0 + 1 ); mean_chain = ( double * ) malloc ( chain_num * sizeof ( double ) ); for ( par_index = 0; par_index < par_num; par_index++ ) { for ( chain_index = 0; chain_index < chain_num; chain_index++ ) { mean_chain[chain_index] = 0.0; for ( k = ind0; k <= gen_index; k++ ) { mean_chain[chain_index] = mean_chain[chain_index] + z[par_index+chain_index*par_num+k*par_num*chain_num]; } mean_chain[chain_index] = mean_chain[chain_index] / rnd0; } mean_all = r8vec_sum ( chain_num, mean_chain ) / ( double ) chain_num; b_var = 0.0; for ( chain_index = 0; chain_index < chain_num; chain_index++ ) { b_var = b_var + pow ( mean_chain[chain_index] - mean_all, 2 ); } b_var = rnd0 * b_var / ( double ) ( chain_num - 1 ); s_sum = 0.0; for ( chain_index = 0; chain_index < chain_num; chain_index++ ) { s = 0.0; for ( k = ind0; k <= gen_index; k++ ) { s = s + pow ( z[par_index+chain_index*par_num+k*par_num*chain_num] - mean_chain[chain_index], 2 ); } s_sum = s_sum + s; } s_sum = s_sum / ( rnd0 - 1.0 ); w_var = s_sum / ( double ) ( chain_num ); var = ( ( rnd0 - 1.0 ) * w_var + b_var ) / rnd0; gr[par_index+(*gr_count)*par_num] = sqrt ( var / w_var ); } /* Set the convergence flag. */ *gr_conv = 1; for ( par_index = 0; par_index < par_num; par_index++ ) { if ( gr_threshold < gr[par_index+(*gr_count)*par_num] ) { *gr_conv = 0; break; } } if ( *gr_conv ) { printf ( "\n" ); printf ( "GR_COMPUTE:\n" ); printf ( " GR convergence at iteration: %d\n", gen_index ); } free ( mean_chain ); *gr_count = *gr_count + 1; return; } /******************************************************************************/ void gr_init ( double gr[], int *gr_conv, int *gr_count, int gr_num, int par_num ) /******************************************************************************/ /* Purpose: GR_INIT initializes Gelman-Rubin variables. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Output, double GR[PAR_NUM*GR_NUM], the Gelman-Rubin statistic. Output, int *GR_CONV, the convergence flag. Output, int *GR_COUNT, counts the number of generations at which the Gelman-Rubin statistic has been computed. Input, int GR_NUM, the number of times the Gelman-Rubin statistic may be computed. Input, int PAR_NUM, the number of parameters. */ { int i; int j; for ( j = 0; j < gr_num; j++ ) { for ( i = 0; i < par_num; i++ ) { gr[i+j*par_num] = 0.0; } } *gr_conv = 0; *gr_count = 0; return; } /******************************************************************************/ void gr_write ( double gr[], char *gr_filename, int gr_num, int par_num, int printstep ) /******************************************************************************/ /* Purpose: GR_WRITE writes Gelman-Rubin R statistics into a file. Licensing: This code is distributed under the GNU LGPL license. Modified: 30 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double GR[PAR_NUM*GR_NUM], the Gelman-Rubin R statistic. Input, char *GR_FILENAME, the Gelman-Rubin filename. Input, int GR_NUM, the number of times the Gelman-Rubin statistic may be computed. Input, int PAR_NUM, the total number of parameters. Input, int PRINTSTEP, the interval between generations on which the Gelman-Rubin statistic will be computed and written to a file. */ { FILE *gr_unit; int i; int j; gr_unit = fopen ( gr_filename, "wt" ); if ( ! gr_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GR_WRITE - Fatal error!\n" ); fprintf ( stderr, " Could not open the file \"%s\"\n", gr_filename ); exit ( 1 ); } fprintf ( gr_unit, "DREAM.C:Monitored_parameter_interchains_Gelman_Rubin_R_statistic\n" ); for ( j = 0; j < gr_num; j++ ) { fprintf ( gr_unit, "%d", printstep * ( j + 1 ) - 1 ); for ( i = 0; i < par_num; i++ ) { fprintf ( gr_unit, " %f", gr[i+j*par_num] ); } fprintf ( gr_unit, "\n" ); } fclose ( gr_unit ); printf ( "\n" ); printf ( "GR_WRITE:\n" ); printf ( " Created the file \"%s\".\n", gr_filename ); return; } /******************************************************************************/ int i4_binomial_sample ( int n, double pp ) /******************************************************************************/ /* Purpose: I4_BINOMIAL_SAMPLE generates a binomial random deviate. Discussion: This procedure generates a single random deviate from a binomial distribution whose number of trials is N and whose probability of an event in each trial is P. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Reference: Voratas Kachitvichyanukul, Bruce Schmeiser, Binomial Random Variate Generation, Communications of the ACM, Volume 31, Number 2, February 1988, pages 216-222. Parameters: Input, int N, the number of binomial trials, from which a random deviate will be generated. 0 < N. Input, double PP, the probability of an event in each trial of the binomial distribution from which a random deviate is to be generated. 0.0 < PP < 1.0. Output, int I4_BINOMIAL_SAMPLE, a random deviate from the distribution. */ { double al; double alv; double amaxp; double c; double f; double f1; double f2; double ffm; double fm; double g; int i; int ix; int ix1; int k; int m; int mp; double p; double p1; double p2; double p3; double p4; double q; double qn; double r; double t; double u; double v; int value; double w; double w2; double x; double x1; double x2; double xl; double xll; double xlr; double xm; double xnp; double xnpq; double xr; double ynorm; double z; double z2; if ( pp <= 0.0 || 1.0 <= pp ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_BINOMIAL_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " PP is out of range.\n" ); exit ( 1 ); } p = r8_min ( pp, 1.0 - pp ); q = 1.0 - p; xnp = ( double ) ( n ) * p; if ( xnp < 30.0 ) { qn = pow ( q, n ); r = p / q; g = r * ( double ) ( n + 1 ); for ( ; ; ) { ix = 0; f = qn; u = r8_uniform_01_sample ( ); for ( ; ; ) { if ( u < f ) { if ( 0.5 < pp ) { ix = n - ix; } value = ix; return value; } if ( 110 < ix ) { break; } u = u - f; ix = ix + 1; f = f * ( g / ( double ) ( ix ) - r ); } } } ffm = xnp + p; m = ffm; fm = m; xnpq = xnp * q; p1 = ( int ) ( 2.195 * sqrt ( xnpq ) - 4.6 * q ) + 0.5; xm = fm + 0.5; xl = xm - p1; xr = xm + p1; c = 0.134 + 20.5 / ( 15.3 + fm ); al = ( ffm - xl ) / ( ffm - xl * p ); xll = al * ( 1.0 + 0.5 * al ); al = ( xr - ffm ) / ( xr * q ); xlr = al * ( 1.0 + 0.5 * al ); p2 = p1 * ( 1.0 + c + c ); p3 = p2 + c / xll; p4 = p3 + c / xlr; /* Generate a variate. */ for ( ; ; ) { u = r8_uniform_01_sample ( ) * p4; v = r8_uniform_01_sample ( ); /* Triangle */ if ( u < p1 ) { ix = xm - p1 * v + u; if ( 0.5 < pp ) { ix = n - ix; } value = ix; return value; } /* Parallelogram */ if ( u <= p2 ) { x = xl + ( u - p1 ) / c; v = v * c + 1.0 - fabs ( xm - x ) / p1; if ( v <= 0.0 || 1.0 < v ) { continue; } ix = x; } else if ( u <= p3 ) { ix = xl + log ( v ) / xll; if ( ix < 0 ) { continue; } v = v * ( u - p2 ) * xll; } else { ix = xr - log ( v ) / xlr; if ( n < ix ) { continue; } v = v * ( u - p3 ) * xlr; } k = abs ( ix - m ); if ( k <= 20 || xnpq / 2.0 - 1.0 <= k ) { f = 1.0; r = p / q; g = ( n + 1 ) * r; if ( m < ix ) { mp = m + 1; for ( i = mp; i <= ix; i++ ) { f = f * ( g / i - r ); } } else if ( ix < m ) { ix1 = ix + 1; for ( i = ix1; i <= m; i++ ) { f = f / ( g / i - r ); } } if ( v <= f ) { if ( 0.5 < pp ) { ix = n - ix; } value = ix; return value; } } else { amaxp = ( k / xnpq ) * ( ( k * ( k / 3.0 + 0.625 ) + 0.1666666666666 ) / xnpq + 0.5 ); ynorm = - ( double ) ( k * k ) / ( 2.0 * xnpq ); alv = log ( v ); if ( alv < ynorm - amaxp ) { if ( 0.5 < pp ) { ix = n - ix; } value = ix; return value; } if ( ynorm + amaxp < alv ) { continue; } x1 = ( double ) ( ix + 1 ); f1 = fm + 1.0; z = ( double ) ( n + 1 ) - fm; w = ( double ) ( n - ix + 1 ); z2 = z * z; x2 = x1 * x1; f2 = f1 * f1; w2 = w * w; t = xm * log ( f1 / x1 ) + ( n - m + 0.5 ) * log ( z / w ) + ( double ) ( ix - m ) * log ( w * p / ( x1 * q )) + ( 13860.0 - ( 462.0 - ( 132.0 - ( 99.0 - 140.0 / f2 ) / f2 ) / f2 ) / f2 ) / f1 / 166320.0 + ( 13860.0 - ( 462.0 - ( 132.0 - ( 99.0 - 140.0 / z2 ) / z2 ) / z2 ) / z2 ) / z / 166320.0 + ( 13860.0 - ( 462.0 - ( 132.0 - ( 99.0 - 140.0 / x2 ) / x2 ) / x2 ) / x2 ) / x1 / 166320.0 + ( 13860.0 - ( 462.0 - ( 132.0 - ( 99.0 - 140.0 / w2 ) / w2 ) / w2 ) / w2 ) / w / 166320.0; if ( alv <= t ) { if ( 0.5 < pp ) { ix = n - ix; } value = ix; return value; } } } return value; } /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ void i4mat_print ( int m, int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4MAT_PRINT prints an I4MAT. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, int A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { i4mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void i4mat_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: I4MAT_PRINT_SOME prints some of an I4MAT. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the GNU LGPL license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, int A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 10 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } /* Print the columns of the matrix, in strips of INCX. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; j2hi = i4_min ( j2hi, n ); j2hi = i4_min ( j2hi, jhi ); fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col:" ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %6d", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ i2lo = i4_max ( ilo, 1 ); i2hi = i4_min ( ihi, m ); for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to INCX) entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %6d", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ int *i4vec_multinomial_sample ( int n, double p[], int ncat ) /******************************************************************************/ /* Purpose: I4VEC_MULTINOMIAL_SAMPLE generates a multinomial random deviate. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Reference: Luc Devroye, Non-Uniform Random Variate Generation, Springer, 1986, ISBN: 0387963057, LC: QA274.D48. Parameters: Input, int N, the number of events, which will be classified into one of the NCAT categories. Input, double P[NCAT-1]. P(I) is the probability that an event will be classified into category I. Thus, each P(I) must be between 0.0 and 1.0. Only the first NCAT-1 values of P must be defined since P(NCAT) would be 1.0 minus the sum of the first NCAT-1 P's. Input, int NCAT, the number of categories. Output, int I4VEC_MULTINOMIAL_SAMPLE[NCAT], a random observation from the multinomial distribution. All IX(i) will be nonnegative and their sum will be N. */ { int i; int icat; int *ix; int ntot; double prob; double ptot; if ( n < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " N < 0\n" ); exit ( 1 ); } if ( ncat <= 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " NCAT <= 1\n" ); exit ( 1 ); } for ( i = 0; i < ncat - 1; i++ ) { if ( p[i] < 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " Some P(i) < 0.\n" ); exit ( 1 ); } if ( 1.0 < p[i] ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " Some 1 < P(i).\n" ); exit ( 1 ); } } ptot = 0.0; for ( i = 0; i < ncat - 1; i++ ) { ptot = ptot + p[i]; } if ( 0.99999 < ptot ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_MULTINOMIAL_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " 1.0 < Sum of P().\n" ); exit ( 1 ); } /* Initialize variables. */ ntot = n; ptot = 1.0; ix = ( int * ) malloc ( ncat * sizeof ( int ) ); for ( i = 0; i < ncat; i++ ) { ix[i] = 0; } /* Generate the observation. */ for ( icat = 0; icat < ncat - 1; icat++ ) { prob = p[icat] / ptot; ix[icat] = i4_binomial_sample ( ntot, prob ); ntot = ntot - ix[icat]; if ( ntot <= 0 ) { return ix; } ptot = ptot - p[icat]; } ix[ncat-1] = ntot; return ix; } /******************************************************************************/ void i4vec_transpose_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4VEC_TRANSPOSE_PRINT prints an I4VEC "transposed". Discussion: An I4VEC is a vector of I4's. Example: A = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 } TITLE = "My vector: " My vector: 1 2 3 4 5 6 7 8 9 10 11 Licensing: This code is distributed under the GNU LGPL license. Modified: 13 December 2012 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, int A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; int ihi; int ilo; int title_len; title_len = strlen ( title ); for ( ilo = 1; ilo <= n; ilo = ilo + 5 ) { ihi = i4_min ( ilo + 5 - 1, n ); if ( ilo == 1 ) { printf ( "%s", title ); } else { for ( i = 1; i <= title_len; i++ ) { printf ( " " ); } } for ( i = ilo; i <= ihi; i++ ) { printf ( "%12d", a[i-1] ); } printf ( "\n" ); } return; } /******************************************************************************/ int *i4vec_zero_new ( int n ) /******************************************************************************/ /* Purpose: I4VEC_ZERO_NEW creates and zeroes an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, int I4VEC_ZERO_NEW[N], a vector of zeroes. */ { int *a; int i; a = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { a[i] = 0; } return a; } /******************************************************************************/ void input_print ( int chain_num, int cr_num, char *input_filename, double gr_threshold, int jumpstep, double limits[], int gen_num, char *ob_filename, int ob_num, int pair_num, int par_num, int printstep, char *restart_string ) /******************************************************************************/ /* Purpose: INPUT_PRINT prints the data from the input file. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 April 2013 Author: John Burkardt Parameters: Input, int CHAIN_NUM, the total number of chains. Input, int CR_NUM, the total number of CR values. Input, char *INPUT_FILENAME, the name of the input file. Input, double GR_THRESHOLD, the convergence tolerance for the Gelman-Rubin statistic. Input, int JUMPSTEP, forces a "long jump" every JUMPSTEP generations. Input, double LIMITS[2*PAR_NUM], lower and upper limits for each parameter. Input, int GEN_NUM, the total number of generations. Input, char *OB_FILENAME, the name of the observation file. Input, int OB_NUM, the total number of observations. Input, int PAIR_NUM, the number of pairs of crossover chains. Input, int PAR_NUM, the total number of parameters. Input, int PRINTSTEP, the interval between generations on which the Gelman-Rubin statistic will be computed and written to a file. Input, char *RESTART_STRING, indicates whether a restart data file should be used to initialize the chains. */ { int j; printf ( "\n" ); printf ( "INPUT_PRINT:\n" ); printf ( " Input data read from \"%s\".\n", input_filename ); printf ( "\n" ); printf ( " Number of parameters\n" ); printf ( " PAR_NUM = %d\n", par_num ); printf ( "\n" ); printf ( " LIMITS: Lower and upper limits for each parameter:\n" ); printf ( "\n" ); printf ( " Index Lower Upper\n" ); printf ( "\n" ); for ( j = 0; j < par_num; j++ ) { printf ( " %5d %14.6g %14.6g\n", j, limits[0+j*2], limits[1+j*2] ); } printf ( "\n" ); printf ( " Number of generations:\n" ); printf ( " GEN_NUM = %d\n", gen_num ); printf ( "\n" ); printf ( " Number of simultaneous chains:\n" ); printf ( " CHAIN_NUM = %d\n", chain_num ); printf ( "\n" ); printf ( " Observation filename:\n" ); printf ( " OB_FILENAME = \"%s\".\n", ob_filename ); printf ( "\n" ); printf ( " Number of observations:\n" ); printf ( " OB_NUM = %d\n", ob_num ); printf ( "\n" ); printf ( " Number of pairs of chains for crossover:\n" ); printf ( " PAIR_NUM = %d\n", pair_num ); printf ( "\n" ); printf ( " Number of crossover values:\n" ); printf ( " CR_NUM = %d\n", cr_num ); printf ( "\n" ); printf ( " Number of steps til a long jump:\n" ); printf ( " JUMPSTEP = %d\n", jumpstep ); printf ( "\n" ); printf ( " Interval between Gelman-Rubin computations:\n" ); printf ( " PRINTSTEP = %d\n", printstep ); printf ( "\n" ); printf ( " Gelman-Rubin convergence tolerance:\n" ); printf ( " GR_THRESHOLD = %g\n", gr_threshold ); printf ( "\n" ); printf ( " Restart string (Y/N):\n" ); printf ( " RESTART_STRING = \"%s\".\n", restart_string ); return; } /******************************************************************************/ void input_read ( int *chain_num, int *cr_num, char *input_filename, double *gr_threshold, int *jumpstep, double limits[], int *gen_num, char *ob_filename, int *ob_num, int *pair_num, int par_num, int *printstep, char *restart_string ) /******************************************************************************/ /* Purpose: INPUT_READ reads the data from the input file. Discussion: This function gave me an extraordinary amount of trouble. Reading irregular text files with C is a real nightmare. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Output, int *CHAIN_NUM, the total number of chains. Output, int *CR_NUM, the total number of CR values. Input, char *INPUT_FILENAME, the name of the input file. Output, double *GR_THRESHOLD, the convergence tolerance for the Gelman-Rubin statistic. Output, int *JUMPSTEP, forces a "long jump" every JUMPSTEP generations. Output, double LIMITS[2*PAR_NUM], lower and upper limits for each parameter. Output, int *GEN_NUM, the total number of generations. Output, char *OB_FILENAME, the name of the observation file. Output, int *OB_NUM, the total number of observations. Output, int *PAIR_NUM, the number of pairs of crossover chains. Input, int PAR_NUM, the total number of parameters. Output, int *PRINTSTEP, the interval between generations on which the Gelman-Rubin statistic will be computed and written to a file. Output, char *RESTART_STRING, indicates whether a restart data file should be used to initialize the chains. */ { FILE *input_unit; int j; double l; char *line; size_t n; double u; n = 255; line = ( char * ) malloc ( n + 1 ); input_unit = fopen ( input_filename, "rt" ); if ( !input_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "INPUT_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open the file \"%s\".\n", input_filename ); exit ( 1 ); } /* Read and ignore line 1. */ fgets ( line, 255, input_unit ); /* Read and ignore line 2. */ fgets ( line, 255, input_unit ); /* Read and ignore line 3. */ fgets ( line, 255, input_unit ); /* Line 4 is actually PAR_NUM lines, each a pair of values. */ for ( j = 0; j < par_num; j++ ) { fscanf ( input_unit, "%lf %lf", &l, &u ); limits[0+j*2] = l; limits[1+j*2] = u; } fgets ( line, 255, input_unit ); /* Read and ignore line 5. */ fgets ( line, 255, input_unit ); /* Line 6 is GEN_NUM. */ fscanf ( input_unit, "%d", gen_num ); fgets ( line, 255, input_unit ); /* Read and ignore line 7. */ fgets ( line, 255, input_unit ); /* Line 8 is CHAIN_NUM. */ fscanf ( input_unit, "%d", chain_num ); fgets ( line, 255, input_unit ); /* Read and ignore line 9. */ fgets ( line, 255, input_unit ); /* Line 10 is OB_FILENAME. */ fscanf ( input_unit, "%s", ob_filename ); fgets ( line, 255, input_unit ); /* Read and ignore line 11. */ fgets ( line, 255, input_unit ); /* Line 12 is OB_NUM. */ fscanf ( input_unit, "%d", ob_num ); fgets ( line, 255, input_unit ); /* Read and ignore line 13. */ fgets ( line, 255, input_unit ); /* Line 14 is PAIR_NUM. */ fscanf ( input_unit, "%d", pair_num ); fgets ( line, 255, input_unit ); /* Read and ignore line 15. */ fgets ( line, 255, input_unit ); /* Line 16 is CR_NUM. */ fscanf ( input_unit, "%d", cr_num ); fgets ( line, 255, input_unit ); /* Read and ignore line 17. */ fgets ( line, 255, input_unit ); /* Line 18 is JUMPSTEP. */ fscanf ( input_unit, "%d", jumpstep ); fgets ( line, 255, input_unit ); /* Read and ignore line 19. */ fgets ( line, 255, input_unit ); /* Line 20 is PRINTSTEP. */ fscanf ( input_unit, "%d", printstep ); fgets ( line, 255, input_unit ); /* Read and ignore line 21. */ fgets ( line, 255, input_unit ); /* Line 22 is GR_THRESHOLD. */ fscanf ( input_unit, "%lf", gr_threshold ); fgets ( line, 255, input_unit ); /* Read and ignore line 23. */ fgets ( line, 255, input_unit ); /* Line 24 is RESTART_STRING. */ fscanf ( input_unit, "%s", restart_string ); fgets ( line, 255, input_unit ); fclose ( input_unit ); free ( line ); return; } /******************************************************************************/ void input_size ( char *input_filename, int *par_num ) /******************************************************************************/ /* Purpose: INPUT_SIZE reads the parameter size variable from the input file. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 April 2013 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Output, int *PAR_NUM, the total number of parameters. */ { char *flag; FILE *input_unit; char line[255]; input_unit = fopen ( input_filename, "rt" ); if ( !input_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "INPUT_SIZE - Fatal error!\n" ); fprintf ( stderr, " Could not open the file \"%s\".\n", input_filename ); exit ( 1 ); } flag = fgets ( line, 255, input_unit ); if ( !flag ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "INPUT_SIZE - Fatal error!\n" ); fprintf ( stderr, " Could not read line 1 of \"%s\".\n", input_filename ); exit ( 1 ); } flag = fgets ( line, 80, input_unit ); if ( !flag ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "INPUT_SIZE - Fatal error!\n" ); fprintf ( stderr, " Could not read line 2 of \"%s\".\n", input_filename ); exit ( 1 ); } sscanf ( line, "%d", par_num ); fclose ( input_unit ); return; } /******************************************************************************/ void jumprate_choose ( double cr[], int cr_index, int cr_num, int gen_index, int jump_dim[], int *jump_num, double *jumprate, double jumprate_table[], int jumpstep, int par_num ) /******************************************************************************/ /* Purpose: JUMPRATE_CHOOSE chooses a jump rate from the jump rate table. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Reference: Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, Dave Higdon, Accelerating Markov Chain Monte Carlo Simulation by Differential Evolution with Self-Adaptive Randomized Subspace Sampling, International Journal of Nonlinear Sciences and Numerical Simulation, Volume 10, Number 3, March 2009, pages 271-288. Parameters: Input, double CR[CR_NUM], the CR values. Input, int CR_INDEX, the index of the CR. 1 <= CR_INDEX <= CR_NUM. Input, int CR_NUM, the total number of CR values. Input, int GEN_INDEX, the current generation. 1 <= GEN_INDEX <= GEN_NUM. Output, int JUMP_DIM[PAR_NUM], the indexes of the parameters to be updated. Output, int *JUMP_NUM, the number of dimensions in which a jump will be made. 0 <= JUMP_NUM <= PAR_NUM. Output, double *JUMPRATE, the jump rate. Input, double JUMPRATE_TABLE[PAR_NUM], the jump rate table. Input, int JUMPSTEP, forces a "long jump" every JUMPSTEP generations. Input, int PAR_NUM, the total number of parameters. */ { int i; double r; /* Determine the dimensions that will be updated. */ *jump_num = 0; for ( i = 0; i < par_num; i++ ) { jump_dim[i] = 0; } for ( i = 0; i < par_num; i++ ) { r = r8_uniform_01_sample ( ); if ( 1.0 - cr[cr_index] < r ) { jump_dim[*jump_num] = i; *jump_num = *jump_num + 1; } } /* Calculate the general jump rate. */ if ( *jump_num == 0 ) { *jumprate = 0.0; } else { *jumprate = jumprate_table[*jump_num-1]; } /* If parameter dimension is 1, 2, or 3, fix the jump rate to 0.6. */ if ( par_num <= 3 ) { *jumprate = 0.6; } /* Determine if a long jump is forced. */ if ( ( gen_index % jumpstep ) == 0 ) { *jumprate = 0.98; } return; } /******************************************************************************/ double *jumprate_table_init ( int pair_num, int par_num ) /******************************************************************************/ /* Purpose: JUMPRATE_TABLE_INIT initializes the jump rate table. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int PAIR_NUM, the number of pairs of crossover chains. Input, int PAR_NUM, the total number of parameters. Output, double JUMPRATE_TABLE_INIT[PAR_NUM], the jumprate table. */ { double c; int i; double *jumprate_table; jumprate_table = ( double * ) malloc ( par_num * sizeof ( double ) ); c = 2.38 / sqrt ( ( double ) ( 2 * pair_num ) ); for ( i = 0; i < par_num; i++ ) { jumprate_table[i] = c / sqrt ( ( double ) ( i + 1 ) ); } return jumprate_table; } /******************************************************************************/ void jumprate_table_print ( double jumprate_table[], int pair_num, int par_num ) /******************************************************************************/ /* Purpose: JUMPRATE_TABLE_PRINT prints the jump rate table. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 April 2013 Author: John Burkardt Parameters: Input, double JUMPRATE_TABLE[PAR_NUM], the jumprate table. Input, int PAIR_NUM, the number of pairs of crossover chains. Input, int PAR_NUM, the total number of parameters. */ { int i; printf ( "\n" ); printf ( "JUMPRATE_TABLE_PRINT\n" ); printf ( "\n" ); printf ( " I Jumprate\n" ); printf ( "\n" ); for ( i = 0; i < par_num; i++ ) { printf ( " %2d %14.6g\n", i, jumprate_table[i] ); } return; } /******************************************************************************/ void mnor_init ( double cov_mat[], double *logdet, int mnor_ind[], double mnor_mean[], int mnor_num, int par_num, double parm[], double prior_den_par_val[] ) /******************************************************************************/ /* Purpose: MNOR_INIT initializes information for multivariate normal prior densities. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Output, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance matrix for a multivariate normal density. Output, double *LOGDET, the logarithm of the determinant of the covariance matrix. Input, int MNOR_IND[MNOR_NUM], the indices of parameters with multivariate density. Output, int MNOR_MEAN[MNOR_NUM], the means associated with the multivariate density. Input, int MNOR_NUM, the number of parameters with a multivariate normal distribution. Input, int PAR_NUM, the total number of parameters. Output, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters for a multivariate normal distribution. Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values associated with prior densities. */ { double a; int i; int j; /* Collect the mean values of the multinormal density. */ for ( i = 0; i < mnor_num; i++ ) { j = mnor_ind[i]; a = prior_den_par_val[0+j*2]; mnor_mean[i] = a; } /* Pack data into PARM. */ mnor_process ( cov_mat, logdet, mnor_mean, mnor_num, parm ); return; } /******************************************************************************/ void mnor_process ( double cov_mat[], double *logdet, double mnor_mean[], int mnor_num, double parm[] ) /******************************************************************************/ /* Purpose: MNOR_PROCESS sets data for the generation of multivariate normal deviates. Discussion: The final form and purpose of this routine is not decided yet. Licensing: This code is distributed under the GNU LGPL license. Modified: 26 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Parameters: Input/output, double COVM[MNOR_NUM*MNOR_NUM]. On input, the covariance matrix of the multivariate distribution. On output, the information has been overwritten. Output, double *LOGDET, the logarithm of the determinant of the covariance matrix. Input, double MNOR_MEAN[MNOR_NUM], the means of the multivariate normal distribution. Input, int MNOR_NUM, the number of parameters with a multivariate normal distribution. Output, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters needed to generate multivariate normal deviates. */ { int i; int icount; int info; int j; /* Store MNOR_NUM. */ parm[0] = mnor_num; /* Store MEANV. */ for ( i = 1; i < mnor_num; i++ ) { parm[i] = mnor_mean[i-1]; } /* Compute the Cholesky decomposition. */ info = r8mat_pofa ( cov_mat, mnor_num, mnor_num ); if ( info != 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "MNOR_PROCESS - Fatal error!\n" ); fprintf ( stderr, " Matrix not positive definite.\n" ); exit ( 1 ); } /* Store the upper half of the Cholesky factor. */ icount = mnor_num + 1; for ( i = 0; i < mnor_num; i++ ) { for ( j = i; j < mnor_num; j++ ) { parm[icount] = cov_mat[i+j*mnor_num]; icount = icount + 1; } } /* Compute the determinant...NOT DONE! */ *logdet = 1.0; return; } /******************************************************************************/ void mnor_read ( double cov_mat[], char *mnor_filename, int mnor_num ) /******************************************************************************/ /* Purpose: MNOR_READ reads a file for multivariate normal prior densities. Licensing: This code is distributed under the GNU LGPL license. Modified: 25 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Output, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance matrix for a multivariate normal density. Input, char *MNOR_FILENAME, the name of the file containing the covariance matrix information. Input, int MNOR_NUM, the number of parameters with a multivariate normal distribution. */ { int cov_dim; FILE *cov_unit; int i; int j; char *line; size_t n; cov_unit = fopen ( mnor_filename, "rt" ); if ( ! cov_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "MNOR_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open the file \"%s\".\n", mnor_filename ); exit ( 1 ); } n = 255; line = ( char * ) malloc ( n + 1 ); /* Read and ignore the first line. */ fgets ( line, 255, cov_unit ); /* Get the dimension from the second line, then discard the line. */ fscanf ( cov_unit, "%d", &cov_dim ); if ( cov_dim != mnor_num ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "MNOR_READ - Fatal error!\n" ); fprintf ( stderr, " Dimension mismatch.\n" ); fprintf ( stderr, " Expecting matrix dimension MNOR_NUM = %d\n", mnor_num ); fprintf ( stderr, " Found matrix dimension COV_DIM = %d\n", cov_dim ); exit ( 1 ); } fgets ( line, 255, cov_unit ); /* Read the matrix. */ for ( i = 0; i < mnor_num; i++ ) { for ( j = 0; j < mnor_num; j++ ) { fscanf ( cov_unit, "%lf", cov_mat+i+j*mnor_num ); } } fclose ( cov_unit ); free ( line ); return; } /******************************************************************************/ void ob_read ( double ob_data[], char *ob_filename, int ob_num ) /******************************************************************************/ /* Purpose: OB_READ reads observational data from a file. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Output, double OB_DATA[OB_NUM], the data. Input, char *OB_FILENAME, the name of the file. Input, int OB_NUM, the total number of observations. */ { int i; char *line; size_t n; FILE *ob_unit; n = 255; line = ( char * ) malloc ( n + 1 ); ob_unit = fopen ( ob_filename, "rt" ); if ( ! ob_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "OB_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open file \"%s\".\n", ob_filename ); exit ( 1 ); } /* Read and ignore the header line. */ fgets ( line, 255, ob_unit ); /* Read OB_NUM lines of data. */ for ( i = 0; i < ob_num; i++ ) { fscanf ( ob_unit, "%lf", ob_data + i ); } fclose ( ob_unit ); return; } /******************************************************************************/ double prior_density ( double cov_mat[], double logdet, int mnor_ind[], double mnor_mean[], int mnor_num, int par_num, int prior_den_par_num[], double prior_den_par_val[], int prior_den_type[], int uni_ind[], int uni_num, double zp[] ) /******************************************************************************/ /* Purpose: PRIOR_DENSITY evaluates the prior density function. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double COV_MAT[MNOR_NUM*MNOR_NUM], the covariance matrix for a multivariate normal density. Input, double LOGDET, the logarithm of the determinant of the multivariate normal covariance matrix. Input, int MNOR_IND[MNOR_NUM], the indices of parameters with a multivariate normal distribution. Input, double MNOR_MEAN[PAR_NUM], the mean of parameters with a multivariate normal distribution. Input, int MNOR_NUM, the number of parameters with a multivariate normal distribution. Input, int PAR_NUM, the total number of parameters. Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of variables associated with the prior density for each parameter. Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values associated with prior densities. Input, int PRIOR_DEN_TYPE[PAR_NUM], the type of the prior density. Input, int UNI_IND[UNI_NUM], indices of parameters associated with univariate distributions. Input, int UNI_NUM, the number of parameters with a univariate distribution. Input, double ZP[PAR_NUM], the argument of the density function. Output, real PRIOR_DENSITY, the value of the prior density function. */ { double a; double b; int i; int j; double *tmp; int tp; double value; value = 1.0; /* Univariate density: */ for ( i = 0; i < uni_num; i++ ) { j = uni_ind[i]; tp = prior_den_type[j]; if ( 1 <= prior_den_par_num[j] ) { a = prior_den_par_val[0+j*2]; } else { a = 0.0; } if ( 2 <= prior_den_par_num[j] ) { b = prior_den_par_val[1+j*2]; } else { b = 0.0; } if ( tp == 1 ) { value = value * r8_beta_pdf ( a, b, zp[j] ); } else if ( tp == 2 ) { value = value * r8_chi_pdf ( a, zp[j] ); } else if ( tp == 3 ) { value = value * r8_exponential_pdf ( a, zp[j] ); } else if ( tp == 4 ) { value = value * r8_gamma_pdf ( a, b, zp[j] ); } else if ( tp == 5 ) { value = value * r8_invchi_pdf ( a, zp[j] ); } else if ( tp == 6 ) { value = value * r8_invgam_pdf ( a, b, zp[j] ); } else if ( tp == 7 ) { value = value * r8_normal_pdf ( a, b, zp[j] ); } else if ( tp == 8 ) { value = value * r8_scinvchi_pdf ( a, b, zp[j] ); } else if ( tp == 9 ) { value = value * r8_uniform_pdf ( a, b, zp[j] ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRIOR_DENSITY - Fatal error!\n" ); fprintf ( stderr, " Unrecognized type for prior density.\n" ); exit ( 1 ); } } /* Calculate density values for multinormal density */ if ( 1 < mnor_num ) { tmp = ( double * ) malloc ( mnor_num * sizeof ( double ) ); for ( i = 0; i < mnor_num; i++ ) { tmp[i] = zp[mnor_ind[i]]; } /* ? I would assume R8_MNOR_PDF would be accessing PARM, not these pieces. */ value = r8_mnor_pdf ( mnor_num, mnor_mean, cov_mat, logdet, tmp ) ; free ( tmp ); } return value; } /******************************************************************************/ void prior_print ( int mnor_ind[], int mnor_num, int par_num, int prior_den_par_num[], double prior_den_par_val[], int prior_den_type[], char *prior_filename, int uni_ind[], int uni_num ) /******************************************************************************/ /* Purpose: PRIOR_PRINT prints information about the prior density. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int MNOR_IND[MNOR_NUM], the indices of parameters with a multivariate normal distribution. Input, int MNOR_NUM, the number of parameters with a multivariate normal distribution. Input, int PAR_NUM, the total number of parameters. Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of variables associated with the prior density for each parameter. Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values associated with prior densities. Input, int PRIOR_DEN_TYPE[PAR_NUM], the type of the prior density. Input, char *PRIOR_FILENAME, the name of the file containing information about the prior density. Input, int UNI_IND[UNI_NUM], the indices of variables with univariate density. Input, int UNI_NUM, the number of parameters with a univariate distribution. */ { int i; int j; int k; int pn; int pt; char type_string[40]; printf ( "\n" ); printf ( "PRIOR_PRINT:\n" ); printf ( " Prior density information read from \"%s\".\n", prior_filename ); printf ( "\n" ); printf ( " Number of univariate parameters = %d\n", uni_num ); printf ( "\n" ); printf ( " Full Uni\n" ); printf ( " Index Index Type #" ); printf ( " Arg1 Arg2\n" ); printf ( "\n" ); for ( i = 0; i < uni_num; i++ ) { j = uni_ind[i]; pt = prior_den_type[j]; if ( pt == 1 ) { strcpy ( type_string, "beta" ); } else if ( pt == 2 ) { strcpy ( type_string, "chi-square" ); } else if ( pt == 3 ) { strcpy ( type_string, "exponential" ); } else if ( pt == 4 ) { strcpy ( type_string, "gamma" ); } else if ( pt == 5 ) { strcpy ( type_string, "inv-chi-square" ); } else if ( pt == 6 ) { strcpy ( type_string, "inv-gamma" ); } else if ( pt == 7 ) { strcpy ( type_string, "normal" ); } else if ( pt == 8 ) { strcpy ( type_string, "scaled-inv-chi-square" ); } else if ( pt == 9 ) { strcpy ( type_string, "uniform" ); } else if ( pt == 10 ) { strcpy ( type_string, "multinormal" ); } else { strcpy ( type_string, "?" ); } pn = prior_den_par_num[j]; printf ( " %5d %5d %20s %1d", i, j, type_string, pn ); for ( k = 0; k < pn; k++ ) { printf ( " %14.6g", prior_den_par_val[k+j*2] ); } printf ( "\n" ); } printf ( "\n" ); printf ( " Number of multivariate parameters = %d\n", mnor_num ); if ( 0 < mnor_num ) { printf ( "\n" ); printf ( " Index Parameter index\n" ); printf ( "\n" ); for ( i = 0; i < mnor_num; i++ ) { j = mnor_ind[i]; printf ( " %5d %2d\n", i, j ); } } return; } /******************************************************************************/ void prior_read ( int mnor_ind[], int *mnor_num, int par_num, int prior_den_par_num[], double prior_den_par_val[], int prior_den_type[], char *prior_filename, int uni_ind[], int *uni_num ) /******************************************************************************/ /* Purpose: PRIOR_READ reads information about the prior density. Licensing: This code is distributed under the GNU LGPL license. Modified: 19 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Output, int MNOR_IND[MNOR_NUM], the indices of parameters with a multivariate normal distribution. Output, int *MNOR_NUM, the number of parameters with a multivariate normal distribution. Input, int PAR_NUM, the total number of parameters. Output, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of variables associated with the prior density for each parameter. Output, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values associated with prior densities. Output, int PRIOR_DEN_TYPE[PAR_NUM], the type of the prior density. Input, char *PRIOR_FILENAME, the name of the file containing information about the prior density. Output, int UNI_IND[UNI_NUM], the indices of variables with univariate density. Output, int *UNI_NUM, the number of parameters with a univariate distribution. */ { int k; char *line; size_t n; int par_index; int pn; double pv; FILE *prior_unit; char type_string[40]; n = 255; line = ( char * ) malloc ( n + 1 ); *mnor_num = 0; *uni_num = 0; prior_unit = fopen ( prior_filename, "rt" ); if ( !prior_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRIOR_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open the file \"%s\".\n", prior_filename ); exit ( 1 ); } /* Read and ignore line 1. */ fgets ( line, 255, prior_unit ); /* Read prior information for each parameter dimension. */ for ( par_index = 0; par_index < par_num; par_index++ ) { fscanf ( prior_unit, "%s", type_string ); fgets ( line, 255, prior_unit ); if ( !strcasecmp ( type_string, "beta" ) ) { k = 1; } else if ( !strcasecmp ( type_string, "chi-square" ) ) { k = 2; } else if ( !strcasecmp ( type_string, "exponential" ) ) { k = 3; } else if ( !strcasecmp ( type_string, "gamma" ) ) { k = 4; } else if ( !strcasecmp ( type_string, "inv-chi-square" ) ) { k = 5; } else if ( !strcasecmp ( type_string, "inv-gamma" ) ) { k = 6; } else if ( !strcasecmp ( type_string, "normal" ) ) { k = 7; } else if ( !strcasecmp ( type_string, "scaled-inv-chi-square" ) ) { k = 8; } else if ( !strcasecmp ( type_string, "uniform" ) ) { k = 9; } else if ( !strcasecmp ( type_string, "multinormal" ) ) { k = 10; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRIOR_READ - Fatal error!\n" ); fprintf ( stderr, " Unrecognized input value of type string.\n" ); exit ( 1 ); } prior_den_type[par_index] = k; /* Read parameters for the prior density function. */ fscanf ( prior_unit, "%d", &pn ); prior_den_par_num[par_index] = pn; if ( 1 <= pn ) { fscanf ( prior_unit, "%lf", &pv ); prior_den_par_val[0+par_index*2] = pv; } if ( 2 <= pn ) { fscanf ( prior_unit, "%lf", &pv ); prior_den_par_val[1+par_index*2] = pv; } fgets ( line, 255, prior_unit ); if ( k == 10 ) { mnor_ind[*mnor_num] = par_index; *mnor_num = *mnor_num + 1; } else { uni_ind[*uni_num] = par_index; *uni_num = *uni_num + 1; } } fclose ( prior_unit ); free ( line ); return; } /******************************************************************************/ double *prior_sample ( int mnor_ind[], int mnor_num, int par_num, double parm[], int prior_den_par_num[], double prior_den_par_val[], int prior_den_type[], int uni_ind[], int uni_num ) /******************************************************************************/ /* Purpose: PRIOR_SAMPLE samples from the prior distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 18 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int MNOR_IND[MNOR_NUM], indices of parameters associated with multivariate normal distributions. Input, int MNOR_NUM, the number of parameters with a multivariate normal distribution. Input, int PAR_NUM, the total number of parameters. Input, double PARM[MNOR_NUM*(MNOR_NUM+3)/2+1], parameters for a multivariate normal distribution. Input, int PRIOR_DEN_PAR_NUM[PAR_NUM], the number of variables associated with the prior density for each parameter. Input, double PRIOR_DEN_PAR_VAL[2*PAR_NUM], parameter values associated with prior densities. Input, int PRIOR_DEN_TYPE[PAR_NUM], the type of the prior density. Input, int UNI_IND[UNI_NUM], indices of parameters associated with univariate distributions. Input, int UNI_NUM, the number of parameters with a univariate distribution. Output, double PRIOR_SAMPLE[PAR_NUM], the sample from the distribution. */ { double a; double a1; double b; double b1; int i; int j; int pt; double *tmp; double *zp; /* Generate samples following the univariate density. */ zp = ( double * ) malloc ( par_num * sizeof ( double ) ); for ( i = 0; i < uni_num; i++ ) { j = uni_ind[i]; pt = prior_den_type[j]; if ( 1 <= prior_den_par_num[j] ) { a = prior_den_par_val[0+j*2]; } else { a = 0.0; } if ( 2 <= prior_den_par_num[j] ) { b = prior_den_par_val[1+j*2]; } else { b = 0.0; } if ( pt == 1 ) { zp[j] = r8_beta_sample ( a, b ); } else if ( pt == 2 ) { zp[j] = r8_chi_sample ( a ); } else if ( pt == 3 ) { zp[j] = r8_exponential_sample ( a ); } else if ( pt == 4 ) { zp[j] = r8_gamma_sample ( a, b ); } else if ( pt == 5 ) { a1 = 0.5; b1 = 0.5 * a; zp[j] = r8_gamma_sample ( a1, b1 ); if ( zp[j] != 0.0 ) { zp[j] = 1.0 / zp[j]; } } else if ( pt == 6 ) { zp[j] = r8_gamma_sample ( a, b ); if ( zp[j] != 0.0 ) { zp[j] = 1.0 / zp[j]; } } else if ( pt == 7 ) { zp[j] = r8_normal_sample ( a, b ); } else if ( pt == 8 ) { a1 = 0.5 * a * b; b1 = 0.5 * a; zp[j] = r8_gamma_sample ( a1, b1 ); if ( zp[j] != 0.0 ) { zp[j] = 1.0 / zp[j]; } } else if ( pt == 9 ) { zp[j] = r8_uniform_sample ( a, b ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRIOR_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " Unrecognized type for prior density.\n" ); exit ( 1 ); } } /* Generate samples following a multinormal density. */ if ( 1 < mnor_num ) { tmp = r8vec_mnor_sample ( parm ); for ( i = 0; i < mnor_num; i++ ) { j = mnor_ind[i]; zp[j] = tmp[i]; } free ( tmp ); } return zp; } /******************************************************************************/ double r8_beta_pdf ( double alpha, double beta, double rval ) /******************************************************************************/ /* Purpose: R8_BETA_PDF evaluates the PDF of a beta distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double ALPHA, BETA, shape parameters. 0.0 < ALPHA, BETA. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_BETA_PDF, the value of the PDF at RVAL. */ { double temp; double value; if ( alpha <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_BETA_PDF - Fatal error!\n" ); fprintf ( stderr, " Parameter ALPHA is not positive.\n" ); exit ( 1 ); } if ( beta <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_BETA_PDF- Fatal error!\n" ); fprintf ( stderr, " Parameter BETA is not positive.\n" ); exit ( 1 ); } if ( rval < 0.0 || 1.0 < rval ) { value = 0.0; } else { temp = r8_gamma_log ( alpha + beta ) - r8_gamma_log ( alpha ) - r8_gamma_log ( beta ); value = exp ( temp ) * pow ( rval, alpha - 1.0 ) * pow ( 1.0 - rval, beta - 1.0 ); } return value; } /******************************************************************************/ double r8_beta_sample ( double aa, double bb ) /******************************************************************************/ /* Purpose: R8_BETA_SAMPLE generates a beta random deviate. Discussion: This procedure returns a single random deviate from the beta distribution with parameters A and B. The density is x^(a-1) * (1-x)^(b-1) / Beta(a,b) for 0 < x < 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Reference: Russell Cheng, Generating Beta Variates with Nonintegral Shape Parameters, Communications of the ACM, Volume 21, Number 4, April 1978, pages 317-322. Parameters: Input, double AA, the first parameter of the beta distribution. 0.0 < AA. Input, double BB, the second parameter of the beta distribution. 0.0 < BB. Output, double R8_BETA_SAMPLE, a beta random variate. */ { double a; double alpha; double b; double beta; double delta; double gamma; double k1; double k2; static double log4 = 1.3862943611198906188; static double log5 = 1.6094379124341003746; double r; double s; double t; double u1; double u2; double v; double value; double w; double y; double z; if ( aa <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_BETA_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " AA <= 0.0\n" ); exit ( 1 ); } if ( bb <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_BETA_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " BB <= 0.0\n" ); exit ( 1 ); } /* Algorithm BB */ if ( 1.0 < aa && 1.0 < bb ) { a = r8_min ( aa, bb ); b = r8_max ( aa, bb ); alpha = a + b; beta = sqrt ( ( alpha - 2.0 ) / ( 2.0 * a * b - alpha ) ); gamma = a + 1.0 / beta; for ( ; ; ) { u1 = r8_uniform_01_sample ( ); u2 = r8_uniform_01_sample ( ); v = beta * log ( u1 / ( 1.0 - u1 ) ); w = a * exp ( v ); z = u1 * u1 * u2; r = gamma * v - log4; s = a + r - w; if ( 5.0 * z <= s + 1.0 + log5 ) { break; } t = log ( z ); if ( t <= s ) { break; } if ( t <= ( r + alpha * log ( alpha / ( b + w ) ) ) ) { break; } } } /* Algorithm BC */ else { a = r8_max ( aa, bb ); b = r8_min ( aa, bb ); alpha = a + b; beta = 1.0 / b; delta = 1.0 + a - b; k1 = delta * ( 1.0 / 72.0 + b / 24.0 ) / ( a / b - 7.0 / 9.0 ); k2 = 0.25 + ( 0.5 + 0.25 / delta ) * b; for ( ; ; ) { u1 = r8_uniform_01_sample ( ); u2 = r8_uniform_01_sample ( ); if ( u1 < 0.5 ) { y = u1 * u2; z = u1 * y; if ( k1 <= 0.25 * u2 + z - y ) { continue; } } else { z = u1 * u1 * u2; if ( z <= 0.25 ) { v = beta * log ( u1 / ( 1.0 - u1 ) ); w = a * exp ( v ); if ( aa == a ) { value = w / ( b + w ); } else { value = b / ( b + w ); } return value; } if ( k2 < z ) { continue; } } v = beta * log ( u1 / ( 1.0 - u1 ) ); w = a * exp ( v ); if ( log ( z ) <= alpha * ( log ( alpha / ( b + w ) ) + v ) - log4 ) { break; } } } if ( aa == a ) { value = w / ( b + w ); } else { value = b / ( b + w ); } return value; } /******************************************************************************/ double r8_chi_pdf ( double df, double rval ) /******************************************************************************/ /* Purpose: R8_CHI_PDF evaluates the PDF of a chi-squared distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C by John Burkardt. Parameters: Input, double DF, the degrees of freedom. 0.0 < DF. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_CHI_PDF, the value of the PDF at RVAL. */ { double temp1; double temp2; double value; if ( df <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_CHI_PDF - Fatal error!\n" ); fprintf ( stderr, " Degrees of freedom must be positive.\n" ); exit ( 1 ); } if ( rval <= 0.0 ) { value = 0.0; } else { temp2 = df * 0.5; temp1 = ( temp2 - 1.0 ) * log ( rval ) - 0.5 * rval - temp2 * log ( 2.0 ) - r8_gamma_log ( temp2 ); value = exp ( temp1 ); } return value; } /******************************************************************************/ double r8_chi_sample ( double df ) /******************************************************************************/ /* Purpose: R8_CHI_SAMPLE generates a Chi-Square random deviate. Discussion: This procedure generates a random deviate from the chi square distribution with DF degrees of freedom random variable. The algorithm exploits the relation between chisquare and gamma. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Parameters: Input, double DF, the degrees of freedom. 0.0 < DF. Output, double R8_CHI_SAMPLE, a random deviate from the distribution. */ { double arg1; double arg2; double value; if ( df <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_CHI_SAMPLE - Fatal error!\n" ); fprintf ( stderr, " DF <= 0.\n" ); fprintf ( stderr, " Value of DF: %g\n", df ); exit ( 1 ); } arg1 = 1.0; arg2 = df / 2.0; value = 2.0 * r8_gamma_sample ( arg1, arg2 ); return value; } /******************************************************************************/ double r8_epsilon ( void ) /******************************************************************************/ /* Purpose: R8_EPSILON returns the R8 round off unit. Discussion: R8_EPSILON is a number R which is a power of 2 with the property that, to the precision of the computer's arithmetic, 1 < 1 + R but 1 = ( 1 + R / 2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 01 September 2012 Author: John Burkardt Parameters: Output, double R8_EPSILON, the R8 round-off unit. */ { static double value = 2.220446049250313E-016; return value; } /******************************************************************************/ double r8_exponential_pdf ( double beta, double rval ) /******************************************************************************/ /* Purpose: R8_EXPONENTIAL_PDF evaluates the PDF of an exponential distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double BETA, the scale value. 0.0 < BETA. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_EXPONENTIAL_PDF, the value of the PDF at RVAL. */ { double value; if ( beta <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_EXPONENTIAL_PDF - Fatal error!\n" ); fprintf ( stderr, " BETA parameter must be positive.\n" ); exit ( 1 ); } if ( rval < 0.0 ) { value = 0.0; } else { value = exp ( - rval / beta ) / beta; } return value; } /******************************************************************************/ double r8_exponential_sample ( double lambda ) /******************************************************************************/ /* Purpose: R8_EXPONENTIAL_SAMPLE samples the exponential PDF. Discussion: Note that the parameter LAMBDA is a multiplier. In some formulations, it is used as a divisor instead. Licensing: This code is distributed under the GNU LGPL license. Modified: 20 April 2013 Author: John Burkardt Parameters: Input, double LAMBDA, the parameter of the PDF. Output, double R8_EXPONENTIAL_SAMPLE, a sample of the PDF. */ { double r; double value; r = r8_uniform_01_sample ( ); value = - log ( r ) * lambda; return value; } /******************************************************************************/ double r8_exponential_01_sample ( ) /******************************************************************************/ /* Purpose: R8_EXPONENTIAL_01_SAMPLE samples the standard exponential PDF. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 April 2013 Author: John Burkardt Parameters: Output, double R8_EXPONENTIAL_01_SAMPLE, a sample of the PDF. */ { double r; double value; r = r8_uniform_01_sample ( ); value = - log ( r ); return value; } /******************************************************************************/ double r8_gamma_log ( double x ) /******************************************************************************/ /* Purpose: R8_GAMMA_LOG evaluates the logarithm of the gamma function. Discussion: 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. Licensing: This code is distributed under the GNU LGPL license. Modified: 19 April 2013 Author: Original FORTRAN77 version by William Cody, Laura Stoltz. C version by John Burkardt. Reference: William Cody, Kenneth Hillstrom, Chebyshev Approximations for the Natural Logarithm of the Gamma Function, Mathematics of Computation, Volume 21, Number 98, April 1967, pages 198-203. Kenneth Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May 1969. John Hart, Ward Cheney, Charles Lawson, Hans Maehly, Charles Mesztenyi, John Rice, Henry Thatcher, Christoph Witzgall, Computer Approximations, Wiley, 1968, LC: QA297.C64. Parameters: Input, double X, the argument of the function. Output, double R8_GAMMA_LOG, the value of the function. */ { double c[7] = { -1.910444077728E-03, 8.4171387781295E-04, -5.952379913043012E-04, 7.93650793500350248E-04, -2.777777777777681622553E-03, 8.333333333333333331554247E-02, 5.7083835261E-03 }; double corr; static double d1 = -5.772156649015328605195174E-01; static double d2 = 4.227843350984671393993777E-01; static double d4 = 1.791759469228055000094023; static double frtbig = 2.25E+76; int i; double p1[8] = { 4.945235359296727046734888, 2.018112620856775083915565E+02, 2.290838373831346393026739E+03, 1.131967205903380828685045E+04, 2.855724635671635335736389E+04, 3.848496228443793359990269E+04, 2.637748787624195437963534E+04, 7.225813979700288197698961E+03 }; double p2[8] = { 4.974607845568932035012064, 5.424138599891070494101986E+02, 1.550693864978364947665077E+04, 1.847932904445632425417223E+05, 1.088204769468828767498470E+06, 3.338152967987029735917223E+06, 5.106661678927352456275255E+06, 3.074109054850539556250927E+06 }; double p4[8] = { 1.474502166059939948905062E+04, 2.426813369486704502836312E+06, 1.214755574045093227939592E+08, 2.663432449630976949898078E+09, 2.940378956634553899906876E+10, 1.702665737765398868392998E+11, 4.926125793377430887588120E+11, 5.606251856223951465078242E+11 }; double q1[8] = { 6.748212550303777196073036E+01, 1.113332393857199323513008E+03, 7.738757056935398733233834E+03, 2.763987074403340708898585E+04, 5.499310206226157329794414E+04, 6.161122180066002127833352E+04, 3.635127591501940507276287E+04, 8.785536302431013170870835E+03 }; double q2[8] = { 1.830328399370592604055942E+02, 7.765049321445005871323047E+03, 1.331903827966074194402448E+05, 1.136705821321969608938755E+06, 5.267964117437946917577538E+06, 1.346701454311101692290052E+07, 1.782736530353274213975932E+07, 9.533095591844353613395747E+06 }; double q4[8] = { 2.690530175870899333379843E+03, 6.393885654300092398984238E+05, 4.135599930241388052042842E+07, 1.120872109616147941376570E+09, 1.488613728678813811542398E+10, 1.016803586272438228077304E+11, 3.417476345507377132798597E+11, 4.463158187419713286462081E+11 }; double res; static double sqrtpi = 0.9189385332046727417803297; static double xbig = 2.55E+305; double xden; static double xinf = 1.79E+308; double xm1; double xm2; double xm4; double xnum; double y; double ysq; y = x; if ( 0.0 < y && y <= xbig ) { if ( y <= r8_epsilon ( ) ) { res = - log ( y ); } /* EPS < X <= 1.5. */ else if ( y <= 1.5 ) { if ( y < 0.6796875 ) { corr = -log ( y ); xm1 = y; } else { corr = 0.0; xm1 = ( y - 0.5 ) - 0.5; } if ( y <= 0.5 || 0.6796875 <= y ) { xden = 1.0; xnum = 0.0; 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 - 0.5 ) - 0.5; xden = 1.0; xnum = 0.0; for ( i = 0; i < 8; i++ ) { xnum = xnum * xm2 + p2[i]; xden = xden * xm2 + q2[i]; } res = corr + xm2 * ( d2 + xm2 * ( xnum / xden ) ); } } /* 1.5 < X <= 4.0. */ else if ( y <= 4.0 ) { xm2 = y - 2.0; xden = 1.0; xnum = 0.0; for ( i = 0; i < 8; i++ ) { xnum = xnum * xm2 + p2[i]; xden = xden * xm2 + q2[i]; } res = xm2 * ( d2 + xm2 * ( xnum / xden ) ); } /* 4.0 < X <= 12.0. */ else if ( y <= 12.0 ) { xm4 = y - 4.0; xden = -1.0; xnum = 0.0; for ( i = 0; i < 8; i++ ) { xnum = xnum * xm4 + p4[i]; xden = xden * xm4 + q4[i]; } res = d4 + xm4 * ( xnum / xden ); } /* Evaluate for 12 <= argument. */ else { res = 0.0; 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 - 0.5 * corr; res = res + y * ( corr - 1.0 ); } } /* Return for bad arguments. */ else { res = xinf; } /* Final adjustments and return. */ return res; } /******************************************************************************/ double r8_gamma_pdf ( double beta, double alpha, double rval ) /******************************************************************************/ /* Purpose: R8_GAMMA_PDF evaluates the PDF of a gamma distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double BETA, the rate parameter. 0.0 < BETA. Input, double ALPHA, the shape parameter. 0.0 < ALPHA. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_GAMMA_PDF, the value of the PDF at RVAL. */ { double temp; double value; if ( alpha <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_GAMMA_PDF - Fatal error!\n" ); fprintf ( stderr, " Parameter ALPHA is not positive.\n" ); exit ( 1 ); } if ( beta <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_GAMMA_PDF - Fatal error!\n" ); fprintf ( stderr, " Parameter BETA is not positive.\n" ); exit ( 1 ); } if ( rval <= 0.0 ) { value = 0.0; } else { temp = alpha * log ( beta ) + ( alpha - 1.0 ) * log ( rval ) - beta * rval - r8_gamma_log ( alpha ); value = exp ( temp ); } return value; } /******************************************************************************/ double r8_gamma_sample ( double a, double r ) /******************************************************************************/ /* Purpose: R8_GAMMA_SAMPLE generates a Gamma random deviate. Discussion: This procedure generates random deviates from the gamma distribution whose density is (A^R)/Gamma(R) * X^(R-1) * Exp(-A*X) Licensing: This code is distributed under the GNU LGPL license. Modified: 22 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Reference: Joachim Ahrens, Ulrich Dieter, Generating Gamma Variates by a Modified Rejection Technique, Communications of the ACM, Volume 25, Number 1, January 1982, pages 47-54. Joachim Ahrens, Ulrich Dieter, Computer Methods for Sampling from Gamma, Beta, Poisson and Binomial Distributions, Computing, Volume 12, Number 3, September 1974, pages 223-246. Parameters: Input, double A, the location parameter. Input, double R, the shape parameter. Output, double R8_GAMMA_SAMPLE, a random deviate from the distribution. */ { double value; value = r8_gamma_01_sample ( r ) / a; return value; } /******************************************************************************/ double r8_gamma_01_sample ( double a ) /******************************************************************************/ /* Purpose: R8_GAMMA_01_SAMPLE samples the standard Gamma distribution. Discussion: This procedure corresponds to algorithm GD in the reference. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Reference: Joachim Ahrens, Ulrich Dieter, Generating Gamma Variates by a Modified Rejection Technique, Communications of the ACM, Volume 25, Number 1, January 1982, pages 47-54. Parameters: Input, double A, the parameter of the standard gamma distribution. 0.0 < A < 1.0. Output, double R8_GAMMA_01_SAMPLE, a random deviate from the distribution. */ { double a1 = 0.3333333; double a2 = -0.2500030; double a3 = 0.2000062; double a4 = -0.1662921; double a5 = 0.1423657; double a6 = -0.1367177; double a7 = 0.1233795; double b; double c; double d; double e; double e1 = 1.0; double e2 = 0.4999897; double e3 = 0.1668290; double e4 = 0.0407753; double e5 = 0.0102930; double p; double q; double q0; double q1 = 0.04166669; double q2 = 0.02083148; double q3 = 0.00801191; double q4 = 0.00144121; double q5 = -0.00007388; double q6 = 0.00024511; double q7 = 0.00024240; double r; double s; double s2; double si; double sqrt32 = 5.6568542494923801952; double t; double u; double v; double value; double w; double x; if ( 1.0 <= a ) { s2 = a - 0.5; s = sqrt ( s2 ); d = sqrt32 - 12.0 * s; // // Immediate acceptance. // t = r8_normal_01_sample ( ); x = s + 0.5 * t; value = x * x; if ( 0.0 <= t ) { return value; } // // Squeeze acceptance. // u = r8_uniform_01_sample ( ); if ( d * u <= t * t * t ) { return value; } r = 1.0 / a; q0 = (((((( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r + q2 ) * r + q1 ) * r; // // Approximation depending on size of parameter A. // if ( 13.022 < a ) { b = 1.77; si = 0.75; c = 0.1515 / s; } else if ( 3.686 < a ) { b = 1.654 + 0.0076 * s2; si = 1.68 / s + 0.275; c = 0.062 / s + 0.024; } else { b = 0.463 + s + 0.178 * s2; si = 1.235; c = 0.195 / s - 0.079 + 0.16 * s; } // // Quotient test. // if ( 0.0 < x ) { v = 0.5 * t / s; if ( 0.25 < fabs ( v ) ) { q = q0 - s * t + 0.25 * t * t + 2.0 * s2 * log ( 1.0 + v ); } else { q = q0 + 0.5 * t * t * (((((( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v + a2 ) * v + a1 ) * v; } if ( log ( 1.0 - u ) <= q ) { return value; } } for ( ; ; ) { e = r8_exponential_01_sample ( ); u = 2.0 * r8_uniform_01_sample ( ) - 1.0; if ( 0.0 <= u ) { t = b + fabs ( si * e ); } else { t = b - fabs ( si * e ); } // // Possible rejection. // if ( t < -0.7187449 ) { continue; } // // Calculate V and quotient Q. // v = 0.5 * t / s; if ( 0.25 < fabs ( v ) ) { q = q0 - s * t + 0.25 * t * t + 2.0 * s2 * log ( 1.0 + v ); } else { q = q0 + 0.5 * t * t * (((((( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v + a2 ) * v + a1 ) * v; } // // Hat acceptance. // if ( q <= 0.0 ) { continue; } if ( 0.5 < q ) { w = exp ( q ) - 1.0; } else { w = (((( e5 * q + e4 ) * q + e3 ) * q + e2 ) * q + e1 ) * q; } // // May have to sample again. // if ( c * fabs ( u ) <= w * exp ( e - 0.5 * t * t ) ) { break; } } x = s + 0.5 * t; value = x * x; } // // Method for A < 1. // else if ( a < 1.0 ) { b = 1.0 + 0.3678794 * a; for ( ; ; ) { p = b * r8_uniform_01_sample ( ); if ( p < 1.0 ) { value = exp ( log ( p ) / a ); if ( value <= r8_exponential_01_sample ( ) ) { break; } } else { value = - log ( ( b - p ) / a ); if ( ( 1.0 - a ) * log ( value ) <= r8_exponential_01_sample ( ) ) { break; } } } } return value; } /******************************************************************************/ double r8_invchi_pdf ( double df, double rval ) /******************************************************************************/ /* Purpose: R8_INVCHI_PDF evaluates the PDF of an inverse chi-squared distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double DF, the degrees of freedom. 0.0 < DF. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_INVCHI_PDF, the value of the PDF at RVAL. */ { double temp1; double temp2; double value; if ( df <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_INVCHI_PDF - Fatal error!\n" ); fprintf ( stderr, " Degrees of freedom must be positive.\n" ); exit ( 1 ); } if ( rval <= 0.0 ) { value = 0.0; } else { temp2 = df * 0.5; temp1 = - temp2 * log ( 2.0 ) - ( temp2 + 1.0 ) * log ( rval ) - 0.5 / rval - r8_gamma_log ( temp2 ); value = exp ( temp1 ); } return value; } /******************************************************************************/ double r8_invgam_pdf ( double beta, double alpha, double rval ) /******************************************************************************/ /* Purpose: R8_INVGAM_PDF evaluates the PDF of an inverse gamma distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double BETA, the rate parameter. 0.0 < BETA. Input, double ALPHA, the shape parameter. 0.0 < ALPHA. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_INVGAM_PDF, the value of the PDF at RVAL. */ { double temp; double value; if ( alpha <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_INVGAM_PDF - Fatal error!\n" ); fprintf ( stderr, " Parameter ALPHA is not positive.\n" ); exit ( 1 ); } if ( beta <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_INVGAM_PDF - Fatal error!\n" ); fprintf ( stderr, " Parameter BETA is not positive.\n" ); exit ( 1 ); } if ( rval <= 0.0 ) { value = 0.0; } else { temp = alpha * log ( beta ) - ( alpha + 1.0 ) * log ( rval ) - beta / rval - r8_gamma_log ( alpha ); value = exp ( temp ); } return value; } /******************************************************************************/ double r8_max ( double x, double y ) /******************************************************************************/ /* Purpose: R8_MAX returns the maximum of two R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 07 May 2006 Author: John Burkardt Parameters: Input, double X, Y, the quantities to compare. Output, double R8_MAX, the maximum of X and Y. */ { double value; if ( y < x ) { value = x; } else { value = y; } return value; } /******************************************************************************/ double r8_min ( double x, double y ) /******************************************************************************/ /* Purpose: R8_MIN returns the minimum of two R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 07 May 2006 Author: John Burkardt Parameters: Input, double X, Y, the quantities to compare. Output, double R8_MIN, the minimum of X and Y. */ { double value; if ( y < x ) { value = y; } else { value = x; } return value; } /******************************************************************************/ double r8_mnor_pdf ( int ndim, double mu[], double cov[], double logdet, double rvals[] ) /******************************************************************************/ /* Purpose: R8_MNOR_PDF evaluates the PDF of a multivariate normal distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int NDIM, the spatial dimension. Input, double MU[NDIM], the mean values. Input, double COV[NDIM*NDIM], the covariance matrix. (This matrix should be symmetric, positive definite. It is overwritten during the action of this double.) Input, double LOGDET, the logarithm of the determinant of the covariance matrix. Input, double RVALS[NDIM], the point where the PDF is evaluated. Output, double R8_MNOR_PDF, the value of the PDF at RVALS. */ { double *covinv; int i; int j; double pi = 3.141592653589793; double rsum; double value; /* Invert the covariance matrix. */ covinv = ( double * ) malloc ( ndim * ndim * sizeof ( double ) ); for ( i = 0; i < ndim; i++ ) { for ( j = 0; j <= i; j++ ) { covinv[j+i*ndim] = cov[j+i*ndim]; } for ( j = i + 1; j < ndim; j++ ) { covinv[j+i*ndim] = cov[i+j*ndim]; } } /* Calculate the PDF of the multivariate normal distribution. */ rsum = 0.0; for ( j = 0; j < ndim; j++ ) { for ( i = 0; i < ndim; i++ ) { rsum = rsum + ( rvals[i] - mu[i] ) * ( rvals[j] - mu[j] ) * covinv[i+j*ndim]; } } rsum = - 0.5 * ( double ) ( ndim ) * log ( 2.0 * pi ) - 0.5 * logdet - 0.5 * rsum; value = exp ( rsum ); free ( covinv ); return value; } /******************************************************************************/ double r8_normal_pdf ( double av, double sd, double rval ) /******************************************************************************/ /* Purpose: R8_NORMAL_PDF evaluates the PDF of a normal distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double AV, the mean value. Input, double SD, the standard deviation. 0.0 < SD. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_NORMAL_PDF, the value of the PDF at RVAL. */ { double pi = 3.141592653589793; double rtemp; double value; if ( sd <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_NORMAL_PDF - Fatal error!\n" ); fprintf ( stderr, " Standard deviation must be positive.\n" ); exit ( 1 ); } rtemp = ( rval - av ) * ( rval - av ) * 0.5 / ( sd * sd ); value = exp ( - rtemp ) / sd / sqrt ( 2.0 * pi ); return value; } /******************************************************************************/ double r8_normal_sample ( double av, double sd ) /******************************************************************************/ /* Purpose: R8_NORMAL_SAMPLE generates a normal random deviate. Discussion: This procedure generates a single random deviate from a normal distribution with mean AV, and standard deviation SD. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Reference: Joachim Ahrens, Ulrich Dieter, Extensions of Forsythe's Method for Random Sampling from the Normal Distribution, Mathematics of Computation, Volume 27, Number 124, October 1973, page 927-937. Parameters: Input, double AV, the mean. Input, double SD, the standard deviation. Output, double R8_NORMAL_SAMPLE, a random deviate from the distribution. */ { double value; value = sd * r8_normal_01_sample ( ) + av; return value; } /******************************************************************************/ double r8_normal_01_sample ( ) /******************************************************************************/ /* Purpose: R8_NORMAL_01_SAMPLE samples the standard normal probability distribution. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. The Box-Muller method is used, which is efficient, but generates two values at a time. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 April 2013 Author: John Burkardt Parameters: Output, double R8_NORMAL_01_SAMPLE, a normally distributed random value. */ { const double pi = 3.14159265358979323; double r1; double r2; static int used = 0; double x; static double y = 0.0; /* If we've used an even number of values so far, generate two more, return one, and save one. */ if ( ( used % 2 )== 0 ) { r1 = r8_uniform_01_sample ( ); r2 = r8_uniform_01_sample ( ); x = sqrt ( -2.0 * log ( r1 ) ) * cos ( 2.0 * pi * r2 ); y = sqrt ( -2.0 * log ( r1 ) ) * sin ( 2.0 * pi * r2 ); } else { x = y; } used = used + 1; return x; } /******************************************************************************/ int r8_round_i4 ( double x ) /******************************************************************************/ /* Purpose: R8_ROUND_I4 rounds an R8, returning an I4. Example: X Value 1.3 1 1.4 1 1.5 1 or 2 1.6 2 0.0 0 -0.7 -1 -1.1 -1 -1.6 -2 Licensing: This code is distributed under the GNU LGPL license. Modified: 25 March 2013 Author: John Burkardt Parameters: Input, double X, the value. Output, int R8_ROUND_I4, the rounded value. */ { int value; if ( x < 0.0 ) { value = - floor ( - x + 0.5 ); } else { value = floor ( x + 0.5 ); } return value; } /******************************************************************************/ double r8_scinvchi_pdf ( double df, double s, double rval ) /******************************************************************************/ /* Purpose: R8_SCINVCHI_PDF: PDF for a scaled inverse chi-squared distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double DF, the degrees of freedom. 0.0 < DF. Input, double S, the scale factor. 0.0 < S. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_SCINVCHI_PDF, the value of the PDF at RVAL. inverse-chi-square distribution. */ { double temp1; double temp2; double value; if ( df <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_SCINVCHI_PDF - Fatal error!\n" ); fprintf ( stderr, " Degrees of freedom must be positive.\n" ); exit ( 1 ); } if ( s <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_SCINVCHI_PDF - Fatal error!\n" ); fprintf ( stderr, " Scale parameter must be positive.\n" ); exit ( 1 ); } if ( rval <= 0.0 ) { value = 0.0; } else { temp2 = df * 0.5; temp1 = temp2 * log ( temp2 ) + temp2 * log ( s ) - ( temp2 * s / rval ) - ( temp2 + 1.0 ) * log ( rval ) - r8_gamma_log ( temp2 ); value = exp ( temp1 ); } return value; } /******************************************************************************/ double r8_uniform_pdf ( double lower, double upper, double rval ) /******************************************************************************/ /* Purpose: R8_UNIFORM_PDF evaluates the PDF of a uniform distribution. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double LOWER, UPPER, the lower and upper range limits. LOWER < UPPER. Input, double RVAL, the point where the PDF is evaluated. Output, double R8_UNIFORM_PDF, the value of the PDF at RVAL. */ { double value; if ( upper <= lower ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_UNIFORM_PDF - Fatal error!\n" ); fprintf ( stderr, " For uniform PDF, the lower limit must be\n" ); fprintf ( stderr, " less than the upper limit\n" ); exit ( 1 ); } if ( rval < lower ) { value = 0.0; } else if ( rval <= upper ) { value = 1.0 / ( upper - lower ); } else { value = 0.0; } return value; } /******************************************************************************/ double r8_uniform_sample ( double low, double high ) /******************************************************************************/ /* Purpose: R8_UNIFORM_SAMPLE generates a uniform random deviate. Discussion: This procedure generates a real deviate uniformly distributed between LOW and HIGH. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Parameters: Input, double LOW, HIGH, the lower and upper bounds. Output, double R8_UNIFORM_SAMPLE, a random deviate from the distribution. */ { double value; value = low + ( high - low ) * r8_uniform_01_sample ( ); return value; } /******************************************************************************/ double r8_uniform_01_sample ( ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01_SAMPLE generates a uniform random deviate from [0,1]. Discussion: This function should be the only way that the package accesses random numbers. Setting OPTION to 0 accesses the R8_UNI_01() function in RNGLIB, for which there are versions in various languages, which should result in the same values being returned. Setting OPTION to 1 in the C version calls the system RNG "rand()". Licensing: This code is distributed under the GNU LGPL license. Modified: 29 April 2013 Author: John Burkardt. Parameters: Output, double R8_UNIFORM_01_SAMPLE, a random deviate. */ { static int option = 0; double value; if ( option == 1 ) { value = ( double ) rand ( ) / ( double ) RAND_MAX; } else { value = r8_uni_01 ( ); } return value; } /******************************************************************************/ double *r8block_zero_new ( int l, int m, int n ) /******************************************************************************/ /* Purpose: R8BLOCK_ZERO_NEW returns a new zeroed R8BLOCK. Discussion: An R8BLOCK is a triple dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 April 2013 Author: John Burkardt Parameters: Input, int L, M, N, the number of rows, columns, and levels. Output, double R8BLOCK_ZERO_NEW[L*M*N], the new zeroed matrix. */ { double *a; int i; int j; int k; a = ( double * ) malloc ( l * m * n * sizeof ( double ) ); for ( k = 0; k < n; k++ ) { for ( j = 0; j < m; j++ ) { for ( i = 0; i < l; i++ ) { a[i+j*l+k*l*m] = 0.0; } } } return a; } /******************************************************************************/ int r8mat_pofa ( double a[], int lda, int n ) /******************************************************************************/ /* Purpose: R8MAT_POFA factors a real symmetric positive definite matrix. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 April 2013 Author: Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart, C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, LINPACK User's Guide, SIAM, (Society for Industrial and Applied Mathematics), 3600 University City Science Center, Philadelphia, PA, 19104-2688. ISBN 0-89871-172-X Parameters: Input/output, double A[LDA*N]. On input, the symmetric matrix to be factored. Only the diagonal and upper triangle are used. On output, an upper triangular matrix R so that A = R'*R where R' is the transpose. The strict lower triangle is unaltered. If INFO /= 0, the factorization is not complete. Input, int LDA, the leading dimension of the array A. Input, int N, the order of the matrix. Output, int R8MAT_POFA, error flag. 0, for normal return. K, signals an error condition. The leading minor of order K is not positive definite. */ { double dot; int i; int info; int j; int k; double s; double t; for ( j = 1; j <= n; j++ ) { s = 0.0; for ( k = 1; k <= j-1; k++ ) { dot = 0.0; for ( i = 0; i < k - 1; i++ ) { dot = dot + a[i+(k-1)*lda] * a[i+(j-1)*lda]; } t = a[k-1+(j-1)*lda] - dot; t = t / a[k-1+(k-1)*lda]; a[k-1+(j-1)*lda] = t; s = s + t * t; } s = a[j-1+(j-1)*lda] - s; if ( s <= 0.0 ) { info = j; return info; } a[j-1+(j-1)*lda] = sqrt ( s ); } info = 0; return info; } /******************************************************************************/ double *r8mat_zero_new ( int m, int n ) /******************************************************************************/ /* Purpose: R8MAT_ZERO_NEW returns a new zeroed R8MAT. Licensing: This code is distributed under the GNU LGPL license. Modified: 26 September 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Output, double R8MAT_ZERO_NEW[M*N], the new zeroed matrix. */ { double *a; int i; int j; a = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a[i+j*m] = 0.0; } } return a; } /******************************************************************************/ double *r8vec_copy_new ( int n, double a1[] ) /******************************************************************************/ /* Purpose: R8VEC_COPY_NEW copies an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], the vector to be copied. Output, double R8VEC_COPY_NEW[N], the copy of A1. */ { double *a2; int i; a2 = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return a2; } /******************************************************************************/ void r8vec_heap_d ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_HEAP_D reorders an R8VEC into a descending heap. Discussion: An R8VEC is a vector of R8's. A heap is an array A with the property that, for every index J, A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices 2*J+1 and 2*J+2 are legal). Diagram: A(0) / \ A(1) A(2) / \ / \ A(3) A(4) A(5) A(6) / \ / \ A(7) A(8) A(9) A(10) Licensing: This code is distributed under the GNU LGPL license. Modified: 31 May 2009 Author: John Burkardt Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the size of the input array. Input/output, double A[N]. On input, an unsorted array. On output, the array has been reordered into a heap. */ { int i; int ifree; double key; int m; /* Only nodes (N/2)-1 down to 0 can be "parent" nodes. */ for ( i = (n/2)-1; 0 <= i; i-- ) { /* Copy the value out of the parent node. Position IFREE is now "open". */ key = a[i]; ifree = i; for ( ; ; ) { /* Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position IFREE. (One or both may not exist because they equal or exceed N.) */ m = 2 * ifree + 1; /* Does the first position exist? */ if ( n <= m ) { break; } else { /* Does the second position exist? */ if ( m + 1 < n ) { /* If both positions exist, take the larger of the two values, and update M if necessary. */ if ( a[m] < a[m+1] ) { m = m + 1; } } /* If the large descendant is larger than KEY, move it up, and update IFREE, the location of the free position, and consider the descendants of THIS position. */ if ( key < a[m] ) { a[ifree] = a[m]; ifree = m; } else { break; } } } /* When you have stopped shifting items up, return the item you pulled out back to the heap. */ a[ifree] = key; } return; } /******************************************************************************/ double *r8vec_mnor_sample ( double parm[] ) /******************************************************************************/ /* Purpose: R8VEC_MNOR_SAMPLE generates a multivariate normal deviate. Discussion: The method is: 1) Generate P independent standard normal deviates - Ei ~ N(0,1) 2) Using Cholesky decomposition find A so that A'*A = COVM 3) A' * E + MEANV ~ N(MEANV,COVM) Note that PARM contains information needed to generate the deviates, and is set up by SETGMN. PARM(1) contains the size of the deviates, P PARM(2:P+1) contains the mean vector. PARM(P+2:P*(P+3)/2+1) contains the upper half of the Cholesky decomposition of the covariance matrix. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2013 Author: Original FORTRAN77 version by Barry Brown, James Lovato. C version by John Burkardt. Parameters: Input, double PARM[P*(P+3)/2+1], parameters set by SETGMN.! Output, double R8VEC_MNOR_SAMPLE[P], a random deviate from the distribution. */ { double ae; int i; int icount; int j; int p; double *work; double *x; p = ( int ) ( parm[0] ); /* Generate P independent normal deviates. */ work = ( double * ) malloc ( p * sizeof ( double ) ); for ( i = 0; i < p; i++ ) { work[i] = r8_normal_01_sample ( ); } /* Compute X = MEANV + A' * WORK */ x = ( double * ) malloc ( p * sizeof ( double ) ); for ( i = 0; i < p; i++ ) { icount = 0; ae = 0.0; for ( j = 0; j <= i; j++ ) { icount = icount + j; ae = ae + parm[i+j*p-icount+p+1] * work[j]; } x[i] = ae + parm[i+1]; } free ( work ); return x; } /******************************************************************************/ void r8vec_sort_heap_a ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_SORT_HEAP_A ascending sorts an R8VEC using heap sort. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 31 May 2009 Author: John Burkardt Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the number of entries in the array. Input/output, double A[N]. On input, the array to be sorted; On output, the array has been sorted. */ { int n1; double temp; if ( n <= 1 ) { return; } /* 1: Put A into descending heap form. */ r8vec_heap_d ( n, a ); /* 2: Sort A. The largest object in the heap is in A[0]. Move it to position A[N-1]. */ temp = a[0]; a[0] = a[n-1]; a[n-1] = temp; /* Consider the diminished heap of size N1. */ for ( n1 = n - 1; 2 <= n1; n1-- ) { /* Restore the heap structure of the initial N1 entries of A. */ r8vec_heap_d ( n1, a ); /* Take the largest object from A[0] and move it to A[N1-1]. */ temp = a[0]; a[0] = a[n1-1]; a[n1-1] = temp; } return; } /******************************************************************************/ double r8vec_sum ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_SUM returns the sum of an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A[N], the vector. Output, double R8VEC_SUM, the sum of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a[i]; } return value; } /******************************************************************************/ void r8vec_transpose_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_TRANSPOSE_PRINT prints an R8VEC "transposed". Discussion: An R8VEC is a vector of R8's. Example: A = (/ 1.0, 2.1, 3.2, 4.3, 5.4, 6.5, 7.6, 8.7, 9.8, 10.9, 11.0 /) TITLE = 'My vector: ' My vector: 1.0 2.1 3.2 4.3 5.4 6.5 7.6 8.7 9.8 10.9 11.0 Licensing: This code is distributed under the GNU LGPL license. Modified: 12 November 2010 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; int ihi; int ilo; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); if ( n <= 0 ) { printf ( " (Empty)\n" ); return; } for ( ilo = 0; ilo < n; ilo = ilo + 5 ) { ihi = i4_min ( ilo + 5, n ); for ( i = ilo; i < ihi; i++ ) { printf ( " %12g", a[i] ); } printf ( "\n" ); } return; } /******************************************************************************/ double *r8vec_zero_new ( int n ) /******************************************************************************/ /* Purpose: R8VEC_ZERO_NEW creates and zeroes an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 25 March 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, double R8VEC_ZERO_NEW[N], a vector of zeroes. */ { double *a; int i; a = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a[i] = 0.0; } return a; } /******************************************************************************/ void restart_read ( int chain_num, double fit[], int gen_num, int par_num, char *restart_read_filename, double z[] ) /******************************************************************************/ /* Purpose: RESTART_READ reads parameter sample data from a restart file. Discussion: Only a single generation (presumably the last one) was written to the file. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 May 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int CHAIN_NUM, the total number of chains. Output, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of each sample. Input, int GEN_NUM, the total number of generations. Input, int PAR_NUM, the total number of parameters. Input, char *RESTART_READ_FILENAME, the name of the restart file. Output, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. */ { int chain_index; int gen_index; int index; int k; char *line; size_t n; int par_index; FILE *restart_unit; n = 255; line = ( char * ) malloc ( n + 1 ); restart_unit = fopen ( restart_read_filename, "rt" ); if ( ! restart_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RESTART_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open the file \"%s\".\n", restart_read_filename ); exit ( 1 ); } /* Assume only one generation. */ gen_index = 0; /* Read and ignore line 1. */ fgets ( line, 255, restart_unit ); for ( chain_index = 0; chain_index < chain_num; chain_index++ ) { index = chain_index+chain_num*gen_index; fscanf ( restart_unit, "%d%lf", &k, fit+index ); for ( par_index = 0; par_index < par_num; par_index++ ) { index = par_index+par_num*chain_index+par_num*chain_num*gen_index; fscanf ( restart_unit, "%lf", z+index ); } } fclose ( restart_unit ); free ( line ); return; } /******************************************************************************/ void restart_write ( int chain_num, double fit[], int gen_num, int par_num, char *restart_write_filename, double z[] ) /******************************************************************************/ /* Purpose: RESTART_WRITE writes a restart file. Discussion: Only data for the final generation is written. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int CHAIN_NUM, the total number of chains. Input, double FIT[CHAIN_NUM*GEN_NUM], the likelihood of each sample. Input, int GEN_NUM, the total number of generations. Input, int PAR_NUM, the total number of parameters. Input, char *RESTART_WRITE_FILENAME, the name of the restart file. Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. */ { int c; int p; FILE *restart_unit; restart_unit = fopen ( restart_write_filename, "wt" ); if ( !restart_unit ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RESTART_WRITE - Fatal error!\n" ); fprintf ( stderr, " Could not open the file \"%s\".\n", restart_write_filename ); exit ( 1 ); } fprintf ( restart_unit, "DREAM.C:Parameter_values_for_restart.\n" ); for ( c = 0; c < chain_num; c++ ) { fprintf ( restart_unit, "%d %14.7g", c, fit[c+(gen_num-1)*chain_num] ); for ( p = 0; p < par_num; p++ ) { fprintf ( restart_unit, " %14.7g", z[p+c*par_num+(gen_num-1)*par_num*chain_num] ); } fprintf ( restart_unit, "\n" ); } fclose ( restart_unit ); printf ( "\n" ); printf ( "RESTART_WRITE:\n" ); printf ( " Created restart file \"%s\".\n", restart_write_filename ); return; } /******************************************************************************/ double *sample_candidate ( int chain_index, int chain_num, double cr[], int cr_index, int cr_num, int gen_index, int gen_num, double jumprate_table[], int jumpstep, double limits[], int pair_num, int par_num, double z[] ) /******************************************************************************/ /* Purpose: SAMPLE_CANDIDATE generates candidate parameter samples. Licensing: This code is distributed under the GNU LGPL license. Modified: 30 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Reference: Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, Dave Higdon, Accelerating Markov Chain Monte Carlo Simulation by Differential Evolution with Self-Adaptive Randomized Subspace Sampling, International Journal of Nonlinear Sciences and Numerical Simulation, Volume 10, Number 3, March 2009, pages 271-288. Parameters: Input, int CHAIN_INDEX, the chain index. 0 <= CHAIN_INDEX < CHAIN_NUM. Input, int CHAIN_NUM, the total number of chains. Input, double CR[CR_NUM], the CR values. Input, int CR_INDEX, the index of the chosen CR value. 0 <= CR_INDEX < CR_NUM. Input, int CR_NUM, the total number of CR values. Input, int GEN_INDEX, the current generation. 0 <= GEN_INDEX < GEN_NUM. Input, int GEN_NUM, the total number of generations. Input, double JUMPRATE_TABLE[PAR_NUM], the jumprate table. Input, int JUMPSTEP, forces a "long jump" every JUMPSTEP generations. Input, double LIMITS[2*PAR_NUM], limits for the parameters. Input, int PAIR_NUM, the number of pairs of crossover chains. Input, int PAR_NUM, the total number of parameters. Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. Output, double SAMPLE_CANDIDATE[PAR_NUM], a candidate parameter sample. Local parameters: Input, int JUMP_DIM[JUMP_NUM], the dimensions in which a jump is to be made. Local, int JUMP_NUM, the number of dimensions in which a jump will be made. 0 <= JUMP_NUM <= PAR_NUM. Local, double JUMPRATE, the jump rate. */ { double av; double b; double *diff; double *eps; int i; int *jump_dim; int jump_num; double jumprate; double *noise_e; int pair[2]; int *r; double r2; double sd; double *zp; /* Used to calculate E following a uniform distribution on (-B,+B). Because B is currently zero, the noise term is suppressed. */ b = 0.0; /* Pick N pairs of other chains for crossover. */ r = ( int * ) malloc ( 2 * pair_num * sizeof ( int ) ); for ( i = 0; i < pair_num; i++ ) { while ( 1 ) { r2 = r8_uniform_01_sample ( ); pair[0] = ( int ) ( r2 * ( double ) chain_num ); r2 = r8_uniform_01_sample ( ); pair[1] = ( int ) ( r2 * ( double ) chain_num ); if ( pair[0] != pair[1] && pair[0] != chain_index && pair[1] != chain_index ) { break; } } r[0+i*2] = pair[0]; r[1+i*2] = pair[1]; } /* Determine the jump rate. */ jump_dim = ( int * ) malloc ( par_num * sizeof ( int ) ); jumprate_choose ( cr, cr_index, cr_num, gen_index, jump_dim, &jump_num, &jumprate, jumprate_table, jumpstep, par_num ); /* Calculate E in equation 4 of Vrugt. */ noise_e = ( double * ) malloc ( par_num * sizeof ( noise_e ) ); for ( i = 0; i < par_num; i++ ) { noise_e[i] = b * ( 2.0 * r8_uniform_01_sample ( ) - 1.0 ); } /* Get epsilon value from multinormal distribution */ eps = ( double * ) malloc ( par_num * sizeof ( double ) ); av = 0.0; sd = 1.0E-10; for ( i = 0; i < par_num; i++ ) { eps[i] = r8_normal_sample ( av, sd ); } /* Generate the candidate sample ZP based on equation 4 of Vrugt. */ diff = diff_compute ( chain_num, gen_index, gen_num, jump_dim, jump_num, pair_num, par_num, r, z ); zp = ( double * ) malloc ( par_num * sizeof ( double ) ); for ( i = 0; i < par_num; i++ ) { zp[i] = z[i+chain_index*par_num+(gen_index-1)*par_num*chain_num]; } for ( i = 0; i < par_num; i++ ) { zp[i] = zp[i] + ( 1.0 + noise_e[i] ) * jumprate * diff[i] + eps[i]; } /* Enforce limits on the sample ZP. */ sample_limits ( limits, par_num, zp ); free ( diff ); free ( eps ); free ( jump_dim ); free ( noise_e ); free ( r ); return zp; } /******************************************************************************/ void sample_limits ( double limits[], int par_num, double zp[] ) /******************************************************************************/ /* Purpose: SAMPLE_LIMITS enforces limits on a sample variable. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, double LIMITS[2*PAR_NUM], the parameter limits. Input, int PAR_NUM, the total number of parameters. Input/output, double ZP[PAR_NUM], a variable, whose entries, if necessary, will be "folded" so that they lie within the limits. */ { int i; double w; for ( i = 0; i < par_num; i++ ) { w = limits[1+i*2] - limits[0+i*2]; if ( w == 0.0 ) { zp[i] = limits[0+i*2]; } else if ( w < 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SAMPLE_LIMITS - Fatal error!\n" ); fprintf ( stderr, " Upper limit less than lower limit.\n" ); exit ( 1 ); } else { while ( zp[i] < limits[0+i*2] ) { zp[i] = zp[i] + w; } while ( limits[1+i*2] < zp[i] ) { zp[i] = zp[i] - w; } } } return; } /******************************************************************************/ double *std_compute ( int chain_num, int gen_index, int gen_num, int par_num, double z[] ) /******************************************************************************/ /* Purpose: STD_COMPUTE computes the current standard deviations, for each parameter. Discussion: The computation encompasses all chains and generations up to the current ones. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 April 2013 Author: Original FORTRAN90 version by Guannan Zhang. C version by John Burkardt. Parameters: Input, int CHAIN_NUM, the total number of chains. Input, int GEN_INDEX, the current generation. 0 <= GEN_INDEX < GEN_NUM. Input, int GEN_NUM, the total number of generations. Input, int PAR_NUM, the total number of parameters. Input, double Z[PAR_NUM*CHAIN_NUM*GEN_NUM], the Markov chain sample data. Output, double STD_COMPUTE[PAR_NUM], the standard deviations. */ { int i; int j; int k; double mean; double *std; std = ( double * ) malloc ( par_num * sizeof ( double ) ); for ( i = 0; i < par_num; i++ ) { mean = 0.0; for ( k = 0; k <= gen_index; k++ ) { for ( j = 0; j < chain_num; j++ ) { mean = mean + z[i+j*par_num+k*par_num*chain_num]; } } mean = mean / ( double ) ( chain_num ) / ( double ) ( gen_index ); std[i] = 0.0; for ( k = 0; k <= gen_index; k++ ) { for ( j = 0; j < chain_num; j++ ) { std[i] = std[i] + pow ( z[i+j*par_num+k*par_num*chain_num] - mean, 2 ); } } std[i] = sqrt ( std[i] / ( double ) ( chain_num * gen_index - 1 ) ); } return std; }