/*------------------------------------------------------ Maximum likelihood estimation of migration rate and effectice population size using a Metropolis-Hastings Monte Carlo algorithm ------------------------------------------------------- R e p o r t e r R O U T I N E S reports things progress=True or verbose Peter Beerli beerli@scs.fsu.edu Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA Copyright 2003-2007 Peter Beerli, Tallahassee FL This software is distributed free of charge for non-commercial use and is copyrighted. Of course, we do not guarantee that the software works and are not responsible for any damage you may cause or have. $Id: reporter.c 754 2007-06-18 16:16:01Z beerli $ -------------------------------------------------------*/ /* \file reporter.c Routines that report progress and also calculates Gelman-Rubin convergence statistic */ #include "migration.h" #include "mcmc.h" #include "fst.h" #include "random.h" #include "tools.h" #include "broyden.h" #include "combroyden.h" #include "options.h" #include "sighandler.h" #include "migrate_mpi.h" #ifdef DMALLOC_FUNC_CHECK #include #endif void both_chain_means (MYREAL *mc, MYREAL *lc, MYREAL *tc, long len, long lastn, long n); void calc_gelmanb (MYREAL *gelmanb, MYREAL *mc, MYREAL *tc, MYREAL *lc, long len, long lastn, long n); void calc_gelmanw (MYREAL *gelmanw, world_fmt * world, MYREAL *mc, MYREAL *tc, long len, long lastn, long n); void calc_gelmanr (MYREAL *gelmanr, MYREAL *gelmanw, MYREAL *gelmanb, long len, long lastn, long n); void calc_average_biggest_gelmanr (MYREAL *gelmanr, long len, MYREAL *meanR, MYREAL *bigR); void print_gelmanr (MYREAL average, MYREAL biggest); MYREAL calc_s (long tthis, MYREAL *tc, world_fmt * world); MYREAL calc_s_bayes (long tthis, MYREAL *tc, world_fmt * world); void chain_means (MYREAL *thischainmeansm, world_fmt * world); void calc_gelmanw2 (MYREAL *gelmanw, MYREAL *s1, MYREAL *s2, long len); void all_chain_means (MYREAL *mc, MYREAL *chainmeans, long *nmeans, long len, long maxreplicate); void calc_allgelmanb (MYREAL *gelmanb, MYREAL *mc, MYREAL *chainmeans, long *nmeans, long len, long maxreplicate); void calc_allgelmanw2 (MYREAL *gelmanw, MYREAL *chain_s, long *nmeans, long len, long maxreplicate); void calc_allgelmanr (MYREAL *gelmanr, MYREAL *gelmanw, MYREAL *gelmanb, long *nmeans, long len, long maxreplicate); //public functions /// /// convergence indicator for ML runs void convergence_check (world_fmt * world, boolean progress) { static MYREAL *lastchainmeans, *chainmeans, *thischainmeans; static MYREAL *gelmanw, *gelmanb, *gelmanr; static boolean done = FALSE; static long len = 0; static long lastn = 0; static boolean first = TRUE; long maxreplicate = (world->options->replicate && world->options->replicatenum > 0) ? world->options->replicatenum : 1; long n = 0; if (world->chains == 1 && maxreplicate <= 1) return; if (world->start) first = TRUE; if (progress || world->options->gelman) { if (!done) { done = TRUE; // len defines the length of arrays that // have to hold all km, kt, p, and mindex means (ml) or parametes (bayes) if(world->options->bayes_infer) len = world->numpop2 + 1; else len = world->numpop2 + world->numpop * 3; lastchainmeans = (MYREAL *) mycalloc (len, sizeof (MYREAL)); thischainmeans = (MYREAL *) mycalloc (len, sizeof (MYREAL)); chainmeans = (MYREAL *) mycalloc (len, sizeof (MYREAL)); gelmanw = (MYREAL *) mycalloc (len, sizeof (MYREAL)); gelmanb = (MYREAL *) mycalloc (len, sizeof (MYREAL)); gelmanr = (MYREAL *) mycalloc (len, sizeof (MYREAL)); } n = -1; // a dummy value because the gelman_ functions expect a dummy last value memset (thischainmeans, 0, sizeof (MYREAL) * len); if (first) { first = FALSE; chain_means (lastchainmeans, world); return; } else { chain_means (thischainmeans, world); both_chain_means (chainmeans, lastchainmeans, thischainmeans, len, lastn, n); calc_gelmanb (gelmanb, chainmeans, thischainmeans, lastchainmeans, len, lastn, n); calc_gelmanw (gelmanw, world, thischainmeans, lastchainmeans, len, lastn, n); calc_gelmanr (gelmanr, gelmanw, gelmanb, len, lastn, n); calc_average_biggest_gelmanr (gelmanr, len, &world->convergence->gelmanmeanRall, &world->convergence->gelmanmaxRall); memcpy (lastchainmeans, thischainmeans, sizeof (MYREAL) * len); lastn = n; } } } void report_values(MYREAL *vec, long len, char text[]) { long i; fprintf(stdout,"%s ",text); for(i=0;iconvergence->chain_means; long *chain_counts = world->convergence->chain_counts; MYREAL *chain_s = world->convergence->chain_s; long *nmeans; if (world->chains == 1 && maxreplicate <= 1) return; // len defines the length of arrays that // have to hold all km, kt, p, and mindex means len = world->numpop2 + 1; nmeans = (long *) mycalloc (maxreplicate, sizeof (long)); gelmanw = (MYREAL *) mycalloc (len, sizeof (MYREAL)); gelmanb = (MYREAL *) mycalloc (len, sizeof (MYREAL)); gelmanr = (MYREAL *) mycalloc (len, sizeof (MYREAL)); chain_averages = (MYREAL *) mycalloc (len, sizeof (MYREAL)); // can be done independently before this function // evaluate the within_variances for each monitored parameter: wp[i] = 1/(n-1) Sum[p[i]- mean[p[i]]] // evaluate the total_within_variance w = 1/m Sum[wp[i]] // needs to be done here // evaluate between_variance for each monitored parameter: bp[i] = Sum[p[j]-mean[p[j]]] // where m is the number of indpendent chains, m = 2 ... many, shall we evaluate all pairs? // evaluate between_variance for all parameters b = n/(m-1) .... [check program? // calculate expected variance V = (1-1/n) W + 1/n B // Gelman-Rubin statistic sqrt(R) = Sqrt[V/W] for(i = 0; i < maxreplicate; i++) { chain_i = &chain_means[i*len]; s_i = &chain_s[i*len]; n_i = chain_counts[i]; nmeans[i] = n_i; //report_values(chain_i, len, "chain_i"); // setting whether this replicate is already recorded or not // would allow to use asynchronous MPI stuff if(chain_i[0] <= 0.) continue; for(j = 0; j < i; j++) { chain_j = &chain_means[j*len]; s_j = &chain_s[j*len]; n_j = chain_counts[j]; //report_values(chain_j, len, "chain_j"); if(chain_j[0] <= 0.) continue; both_chain_means (chain_averages, chain_i, chain_j, len, n_i, n_j); //report_values(chain_averages, len, "averages"); calc_gelmanb (gelmanb, chain_averages, chain_i, chain_j, len, n_i, n_j); //report_values(gelmanb, len, "gelmanb"); calc_gelmanw2 (gelmanw, s_i, s_j, len); //report_values(gelmanw, len, "gelmanw"); calc_gelmanr (gelmanr, gelmanw, gelmanb, len, n_i, n_j); //report_values(gelmanr, len, "gelmanr"); calc_average_biggest_gelmanr (gelmanr, len, &world->convergence->gelmanmeanmaxR[j * maxreplicate + i], &world->convergence->gelmanmeanmaxR[i * maxreplicate + j]); } } all_chain_means(chain_averages,chain_means, nmeans, len, maxreplicate); calc_allgelmanb(gelmanb,chain_averages, chain_means, nmeans, len, maxreplicate); calc_allgelmanw2 (gelmanw, chain_s, nmeans, len, maxreplicate); calc_allgelmanr (gelmanr, gelmanw, gelmanb, nmeans, len, maxreplicate); calc_average_biggest_gelmanr (gelmanr, len, &world->convergence->gelmanmeanRall, &world->convergence->gelmanmaxRall); myfree(gelmanw); myfree(gelmanb); myfree(gelmanr); myfree(chain_averages); myfree(nmeans); } void both_chain_means (MYREAL *mc, MYREAL *lc, MYREAL *tc, long len, long lastn, long n) { long i; for (i = 0; i < len; i++) { mc[i] = (lc[i] * lastn + tc[i] * n) / (n + lastn); } } /// /// calculates the overall means of all replicate chains void all_chain_means (MYREAL *mc, MYREAL *chainmeans, long *nmeans, long len, long maxreplicate) { long i; long j; long nsum = 0; MYREAL sum = 0.; for (i = 0; i < len; i++) { sum = 0.; nsum = 0; for(j=0; j < maxreplicate; j++) { if(nmeans[j]>1) { sum += chainmeans[j * len + i] * nmeans[j]; nsum += nmeans[j]; } } mc[i] = sum / nsum; } } void calc_gelmanb (MYREAL *gelmanb, MYREAL *mc, MYREAL *tc, MYREAL *lc, long len, long lastn, long n) { long i; for (i = 0; i < len; i++) { gelmanb[i] = /* 1/(2-1)* */ (pow ((lc[i] - mc[i]), 2.) + (pow ((tc[i] - mc[i]), 2.))); } } void calc_allgelmanb (MYREAL *gelmanb, MYREAL *mc, MYREAL *chainmeans, long *nmeans, long len, long maxreplicate) { long i; long j; //long nsum = 0;; MYREAL sum = 0.; MYREAL val; for (i = 0; i < len; i++) { sum = 0.; //nsum = 0; for(j=0; j < maxreplicate; j++) { if(nmeans[j]>1) { val = chainmeans[j * len + i] - mc[i]; sum += val * val; //nsum += nmeans[j]; } } gelmanb[i] = sum / (maxreplicate -1); //gelmanb[i] = nn * (pow ((lc[i] - mc[i]), 2.) + (pow ((tc[i] - mc[i]), 2.))); } } /// /// collect the autocorrelation values and the ESS values void collect_eff_values(world_fmt *world) { const long nn = world->numpop2 + world->bayes->mu; long i; long maxrep = (world->options->replicate == TRUE) ? ((world->options->replicatenum > 0) ? world->options->replicatenum : world->options->lchains ) : 1; for(i=0;iauto_archive[i] += world->autocorrelation[i]/maxrep; world->ess_archive[i] += world->effective_sample[i]; } } /// /// prints order of parameters for header files in bayesfile and bayesallfile void print_param_eff_size(FILE *file, world_fmt *world, long numparam, int offset, MYREAL *autocorr, MYREAL *effsample) { long pa; long pa0; long numpop = world->numpop; long numpop2 = world->numpop2; long frompop, topop; long frompop2, topop2; long pos=0; char *custm2 = world->options->custm2; bayes_fmt * bayes = world->bayes; char *buffer; if(!world->options->bayes_infer) return; buffer = (char*) mycalloc(numparam * LINESIZE, sizeof(char)); pos = sprintf(buffer,"%*.*sParameter Autocorrelation Effective Sample size\n", offset, offset, " "); pos += sprintf(buffer+pos,"%*.*s--------- --------------- ---------------------\n", offset, offset, " "); for(pa0=0;pa0map[pa0][1] != INVALID) { if(pa0 == bayes->map[pa0][1]) pa = pa0; else pa = bayes->map[pa0][1]; if(pa0 < numpop) { if((pa0 == pa) && (custm2[pa0] == '*')) pos += sprintf(buffer+pos,"%*.*s%s_%li % -5.5f %10.2f\n", offset, offset, " ", "Theta",pa0+1,autocorr[pa],effsample[pa]); else pos += sprintf(buffer+pos,"%*.*s%s_%li=%li [%c] % -5.5f %10.2f\n", offset, offset, " ", "Theta",pa0+1, pa+1, custm2[pa0],autocorr[pa],effsample[pa]); } else { // do we estimate mutation rate changes? if(pa0 == numpop2) { pos += sprintf(buffer+pos,"%*.*s%s % -5.5f %10.2f\n", offset, offset, " ", "Rate",autocorr[pa],effsample[pa]); } else { m2mm(pa0,numpop,&frompop,&topop); if((pa0==pa) && (custm2[pa0]=='*')) { pos += sprintf(buffer+pos,"%*.*s%s_(%li,%li) % -5.5f %10.2f\n", offset, offset, " ", "M", frompop+1, topop+1,autocorr[pa],effsample[pa]); } else { m2mm(pa,numpop,&frompop2,&topop2); pos += sprintf(buffer+pos,"%*.*s%s_(%li,%li) [%c] % -5.5f %10.2f\n", offset, offset, " ", "M", frompop+1, topop+1, custm2[pa0], autocorr[pa], effsample[pa]); } } } } } FPRINTF(file,"%s\n", buffer); myfree(buffer); } /// /// returns variance and other measurement of a single chain used in Bayesian burnin_chain() /// the standard deviation and the mean are calculated using a result of /// B. P. Welford 1962 (from D. E. Knuth Art of computer programming Vol 2 page 232 /// the recursion form of the autocorrelation r was developed with Koffi Sampson /// April 2007, see mathematica code below /// newff[xx_] := Module[{}, /// nt = Length[xx]; /// r = 0; /// x1 = xx[[1]]; /// mo = x1; /// mn = x1; /// n = 2 ; /// s = 0; /// While[n <= nt, /// xo = xx[[n - 1]]; /// xn = xx[[n]]; /// mn = mn + (xn - mn)/n; /// s = s + (xn - mo) (xn - mn); /// r = r + xo xn + mo (n mo - x1 - xo) - mn ((n + 1) mn - x1 - xn); /// mo = mn; /// n++; /// ]; /// r/s /// ] /// The effective sample is calculated as the length the sample n * (1-r)/(1+r) /// MYREAL single_chain_var(world_fmt *world, long T, MYREAL *variance, MYREAL *autoc, MYREAL *effsample) { static MYREAL *mean; static MYREAL *S; static MYREAL *r; static MYREAL *xold; static MYREAL *xstart; static boolean done=FALSE; static long n = 0; static long nn = 0; long i; long j; MYREAL x = 0.0; MYREAL xo = 0.0; MYREAL x1 = 0.0; MYREAL m = 0.0; MYREAL mo = 0.0; MYREAL v = 0.0; MYREAL *delta; if(!done) { nn = world->numpop2 + world->bayes->mu; mean = mycalloc(nn,sizeof(MYREAL)); S = mycalloc(nn,sizeof(MYREAL)); r = mycalloc(nn,sizeof(MYREAL)); xold = mycalloc(nn,sizeof(MYREAL)); xstart = mycalloc(nn,sizeof(MYREAL)); done=TRUE; } // reset static variable for each chain if(world==NULL) { n=0; memset(mean,0,nn * sizeof(MYREAL)); memset(S,0,nn * sizeof(MYREAL)); memset(r,0,nn * sizeof(MYREAL)); memset(xold,0,nn * sizeof(MYREAL)); memset(xstart,0,nn * sizeof(MYREAL)); return 0; } // hack: change syntax so that an init really works if(T==0) return 0; delta = mycalloc(nn,sizeof(MYREAL)); if(n==0) { n += 1; for(i=0; i < nn; i++) { if(world->bayes->map[i][1] == INVALID) continue; else { j = world->bayes->map[i][1]; } x = world->param0[j]; delta[j] = x - mean[j]; mean[j] += delta[j]/n; S[j] += delta[j] * (x - mean[j]); v += S[j]; } memcpy(xstart,world->param0,sizeof(MYREAL) * world->numpop2); memcpy(xold,world->param0,sizeof(MYREAL) * world->numpop2); if(world->bayes->mu) { xstart[world->numpop2] = world->options->mu_rates[world->locus]; xold[world->numpop2] = world->options->mu_rates[world->locus]; } myfree(delta); } else { // n is at least 1 n += 1; for(i=0; i < nn; i++) { if(world->bayes->map[i][1] == INVALID) continue; else { j = world->bayes->map[i][1]; } x = world->param0[j]; xo = xold[j]; x1 = xstart[j]; delta[j] = x - mean[j]; mo = mean[j]; mean[j] += delta[j]/n; m = mean[j]; S[j] += delta[j] * (x - m); r[j] += xo * x + mo * (n * mo - x1 - xo) - m * ((n+1) * m - x1 - x); autoc[j] = S[j] > 0.0 ? r[j]/S[j]: 1.0; effsample[j] = n * (1. - autoc[j])/(1. + autoc[j]); v += S[j]; } // printf("n=%li effsample=%g atoc=%g r=%g s=%g mean=%g theta1=%g\n", n,effsample[0],autoc[0],r[0],S[0],mean[0],world->param0[0]); memcpy(xold,world->param0,sizeof(MYREAL) * world->numpop2); if(world->bayes->mu) { xold[world->numpop2] = world->options->mu_rates[world->locus]; } myfree(delta); } *variance = v / n; return *variance; } /// calculate effective sample size /// which is the N' = N (1-r)/(1+r) /// where r is the autocorrelation coefficient for lag 1. /// with r(1) = sum[(x[i]-xmean)(x[i+1]-xmean)/standarddeviation(x),i,0,n-1] void calc_gelmanw (MYREAL *gelmanw, world_fmt * world, MYREAL *mc, MYREAL *tc, long len, long lastn, long n) { long i; MYREAL s1, s2; for (i = 0; i < len; i++) { s1 = calc_s (i, tc, world); s2 = calc_s (i, mc, world); gelmanw[i] = 0.5 * (s1 + s2); } } void calc_gelmanw2 (MYREAL *gelmanw, MYREAL *s1, MYREAL *s2, long len) { long i; for (i = 0; i < len; i++) { gelmanw[i] = 0.5 * (s1[i] + s2[i]); } } /// /// calculates the overall s^2 of all replicate chains void calc_allgelmanw2 (MYREAL *gelmanw, MYREAL *chain_s, long *nmeans, long len, long maxreplicate) { long i; long j; MYREAL sum = 0.; for (i = 0; i < len; i++) { sum = 0.; for(j=0; j < maxreplicate; j++) { if(nmeans[j]>1) { sum += chain_s[j * len + i]; } } gelmanw[i] = sum / maxreplicate; } } void calc_gelmanr (MYREAL *gelmanr, MYREAL *gelmanw, MYREAL *gelmanb, long len, long lastn, long n) { long i; MYREAL nn = (n + lastn) / 2.; MYREAL sqplus = 0.; // MYREAL v; for (i = 0; i < len; i++) { sqplus = ((nn - 1.) / nn * gelmanw[i] + gelmanb[i]); //v = 1.5 * sqplus; gelmanr[i] = sqrt (sqplus / gelmanw[i]);// - (nn-1)/(2. * nn)); } } void calc_allgelmanr (MYREAL *gelmanr, MYREAL *gelmanw, MYREAL *gelmanb, long *nmeans, long len, long maxreplicate) { long i; long j; long nn=0; MYREAL sqplus; // MYREAL v; for(j=0; j < maxreplicate; j++) { nn += nmeans[j]; } nn /= maxreplicate; for (i = 0; i < len; i++) { sqplus = (nn - 1.) / nn * gelmanw[i] + gelmanb[i]; // v = (maxreplicate+1)/maxreplicate * sqplus; gelmanr[i] = sqrt (sqplus / gelmanw[i]);// - (nn-1)/(maxreplicate * nn)); } } void calc_average_biggest_gelmanr (MYREAL *gelmanr, long len, MYREAL *meanR, MYREAL *bigR) { long i; MYREAL average = 0; MYREAL biggest = 0.; for (i = 0; i < len; i++) { if (biggest < gelmanr[i]) biggest = gelmanr[i]; average += gelmanr[i]; } if (len > 0) *meanR = average / len; else *meanR = average; *bigR = biggest; } /// /// calculate the variance of a chain tc is the parameter average MYREAL calc_s (long tthis, MYREAL *tc, world_fmt * world) { timearchive_fmt * atl = &world->atl[world->rep][world->locus]; long T = atl->T; long numpop = world->numpop; // long numpop2 = numpop * numpop; long startp; long startl; long startkm; long i; long j; MYREAL xx; MYREAL s = 0.0; startkm = numpop; startp = startkm + numpop; startl = startp + numpop; if (tthis < startkm) { i = tthis; for (j = 0; j < T; j++) { xx = atl->tl[j].kt[i] - tc[i]; s += xx * xx; } if(s>0.) s /= T - 1.; return s; } else { if (tthis < startp) { i = tthis - startkm; for (j = 0; j < T; j++) { xx = atl->tl[j].km[i] - tc[i]; s += xx * xx; } s /= T - 1.; return s; } else { if (tthis < startl) { i = tthis - startp; for (j = 0; j < T; j++) { xx = atl->tl[j].p[i] - tc[i]; s += xx * xx; } s /= T - 1.; return s; } else { i = tthis - startl; for (j = 0; j < T; j++) { xx = atl->tl[j].mindex[i] - tc[i]; s += xx * xx; } s /= T - 1.; return s; } } } } /// /// calculate the variance of a chain tc is the parameter average /// by supplying the length of the chain MYREAL calc_s2 (long tthis, MYREAL *tc, world_fmt * world, long T) { timearchive_fmt * atl = &world->atl[world->rep][world->locus]; long numpop = world->numpop; // long numpop2 = numpop * numpop; long startp; long startl; long startkm; long i; long j; MYREAL xx; MYREAL s = 0.0; startkm = numpop; startp = startkm + numpop; startl = startp + numpop; if (tthis < startkm) { i = tthis; for (j = 0; j < T; j++) { xx = atl->tl[j].kt[i] - tc[i]; s += xx * xx; } if(s>0.) s /= T - 1.; return s; } else { if (tthis < startp) { i = tthis - startkm; for (j = 0; j < T; j++) { xx = atl->tl[j].km[i] - tc[i]; s += xx * xx; } s /= T - 1.; return s; } else { if (tthis < startl) { i = tthis - startp; for (j = 0; j < T; j++) { xx = atl->tl[j].p[i] - tc[i]; s += xx * xx; } s /= T - 1.; return s; } else { i = tthis - startl; for (j = 0; j < T; j++) { xx = atl->tl[j].mindex[i] - tc[i]; s += xx * xx; } s /= T - 1.; return s; } } } } /// /// calculate the variance of a chain tc is the parameter average MYREAL calc_s_bayes (long tthis, MYREAL *tc, world_fmt * world) { long T = world->convergence->chain_counts[world->rep]; long nn = 2+world->numpop2 + (world->bayes->mu ? 1 : 0) ; // long pnum = world->bayes->numparams; MYREAL * params = world->bayes->params; long i; long j; MYREAL xx; MYREAL s = 0.0; i = tthis; for (j = 0; j < T; j++) { xx = params[j * nn + i] - tc[i]; s += xx * xx; } if(s>0.) s /= T - 1.; return s; } /// /// calculates means for the different events on a tree for the Gelman-Rubin statistic /// this version of the statistics operates on the compressed treedata and not the parameters /// void chain_means_ml (MYREAL *thischainmeans, world_fmt * world) { timearchive_fmt * atl = &world->atl[world->rep][world->locus]; long T = atl->T; long numpop = world->numpop; long numpop2 = numpop * numpop; long startp; long startl; long startkm; long i; long j; startkm = numpop; startp = startkm + numpop; startl = startp + numpop; for (j = 0; j < T; j++) { for (i = 0; i < numpop; i++) { thischainmeans[i] += atl->tl[j].kt[i]; thischainmeans[i + startkm] += atl->tl[j].km[i]; thischainmeans[i + startp] += atl->tl[j].p[i]; } for (i = 0; i < numpop2; i++) { thischainmeans[i + startl] += atl->tl[j].mindex[i]; } } for (i = 0; i < numpop; i++) { thischainmeans[i] /= T; thischainmeans[i + startkm] /= T; thischainmeans[i + startp] /= T; } for (i = startl; i < numpop2 + startl; i++) { thischainmeans[i] /= T; } } /// /// calculates means for the different parameters for the Gelman-Rubin statistic /// for bayesian inference void chain_means_bayes (MYREAL *thischainmeans, world_fmt * world) { long T = world->convergence->chain_counts[world->rep]; long i; long j; MYREAL *params = world->bayes->params; long nn = 2+world->numpop2 + (world->bayes->mu ? 1 : 0) ; for (j = 0; j < T; j++) { for (i = 2; i < nn; i++) { thischainmeans[i-2] += params[j * nn + i]; } } for (i = 2; i < nn; i++) { thischainmeans[i-2] /= T; } } /// /// chain_means void chain_means (MYREAL *thischainmeans, world_fmt * world) { if(world->options->bayes_infer) { world->convergence->chain_counts[world->rep] = world->bayes->paramnum; chain_means_bayes (thischainmeans, world); } else { world->convergence->chain_counts[world->rep] = world->atl[world->rep][world->locus].T; chain_means_ml (thischainmeans, world); } } /// /// calculates the average over void chain_s(MYREAL *cs, MYREAL *cm, world_fmt *world, long replicate, long T) { long len; long start; long stop; long i; if(world->options->bayes_infer) { len = world->numpop2+1; start = replicate * len; stop = start + len; for(i = start; i < stop; i++) { cs[i] = calc_s_bayes (i-start, &cm[start], world); } } else { len = world->numpop2 + 3*world->numpop; start = replicate * len; stop = start + len; for(i = start; i < stop; i++) { cs[i] = calc_s (i-start, &cm[start], world); } } } /// /// report convergence statistic to screen or file void convergence_progress(FILE *file, world_fmt *world) { long i; long j; long maxreplicate = (world->options->replicate && world->options->replicatenum > 0) ? world->options->replicatenum : 1; MYREAL *gelmanmeanR = world->convergence->gelmanmeanmaxR; //upper half of matrix MYREAL *gelmanmaxR = world->convergence->gelmanmeanmaxR; //lower half of matrix fprintf(file,"\n\nGelman-Rubin convergence statistic\n"); fprintf(file,"----------------------------------\n\n"); fprintf(file,"Values close to 1.0, especially values < 1.2 are a sign of convergence of the\n"); fprintf(file,"chains. On very short chains this statistic does not work well\n\n"); if(world->options->gelmanpairs) { fprintf(file,"Pairwise evaluation of replicates\n[above diagonal: mean; below diagonal: maximum value]\n"); fprintf(file,"%4li ", 1L); for(i=1; i < maxreplicate; i++) { fprintf(file,"%6li ", i+1); } fprintf(file,"\n"); for(i=0; i < maxreplicate; i++) { for(j=0; j < maxreplicate; j++) { if(i==j) fprintf(file,"------ "); else { if(i < j) { if(gelmanmeanR[i * maxreplicate + j] <=0.) { fprintf(file," -N- "); } else { fprintf(file," %1.3f ",gelmanmeanR[i * maxreplicate + j]); } } else { if(gelmanmaxR[i * maxreplicate + j] <=0.) { fprintf(file," -N- "); } else { fprintf(file,"*%1.3f ",gelmanmaxR[i * maxreplicate + j]); } } } } fprintf(file,"\n"); } fprintf(file,"\n\n"); } fprintf(file,"Overall convergence:\n"); fprintf(file,"Mean sqrt(R) = %f\n", world->convergence->gelmanmeanRall); fprintf(file,"Maximum sqrt(R) = %f\n\n\n", world->convergence->gelmanmaxRall); }