Bayesian Fusion MCMC
 Performs a hierarchical Bayesian linear regression analysis for heteroscedastic Gaussian measurement errors in both independent (x) and dependent (y) coordinates. 
 Allows for additive Gaussian intrinsic scatter in the y coordinate.
Assumes an informative Gaussian prior on true x coord (xt).
Analytic integration over true x and y coord.
for Mathematica 8
See Section 2.3 of the ''Supplement to Bayesian Logical Data Analysis for the Physical Sciences"

P. C. Gregory
Department of Physics and Astronomy
University of British Columbia

May 14, 2016

All users of this software are requested to include an acknowledgement 
and the following reference in their publications.
      Ref: P. C. Gregory, "Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with (Mathematica^®) Support,"  
      Cambridge University Press 2010
      Includes the 'Supplement to Bayesian Logical Data Analysis for the Physical Sciences,' 
      in the Resources section at
http://www.cambridge.org/pl/academic/subjects/statistics-probability/statistics-physical-sciences-and-engineering/bayesian-logical-data-analysis-physical-sciences-comparative-approach-mathematica-support

Disclaimer

Disclaimer of Warranty
We make no warranties, express or implied, that this Mathematica software is free of error, or is consistent with any particular standard of merchantability, or that it will meet your requirements for any particular application.This software should not be relied on for solving a problem whose incorrect solution could result in injury to a person or loss of property. If you do use the software in such a manner, it is at your own risk. The author and publisher disclaim all liability for direct or consequential damages resulting from your use of this software.

Introduction

Fusion Markov chain Monte Carlo (FMCMC) is a new general purpose Metropolis algorithm to facilitate Bayesian nonlinear model fitting. It incorporates parallel tempering, simulated annealing and genetic crossover operations. Each of these features facilitate the detection of a global minimum in chi-squared in a highly multi-modal environment. By combining all three, the algorithm greatly increases the probability of realizing this goal. FMCMC was originally developed to provide a multi-planet periodogram for searching for multiple planets in the precision radial velocity. That problem is characterized by sparse data and frequently yields multi-modal solutions that are explored via parallel tempering. In this notebook, FMCMC is used for linear regression analysis with measurement errors in both coordinates. This analysis analytically marginalizes over the unknown true coordinates.

The FMCMC is controlled by a unique adaptive control system that automates the tuning of the proposal distributions for efficient exploration of the model parameter space even when the parameters are highly correlated. This controlled statistical fusion approach has the potential to integrate other relevant  statistical tools as desired. The FMCMC is implemented in a Mathematica notebook. On an 8 core PC, the automatic parallelization results in a speed increase of 7 times.

The roots of this notebook go back to a semi-automatic parallel tempering (Section 12.8 of my book, 'Bayesian Logical Data Analysis for the Physical Sciences', Cambridge University Press, 2005) version of Markov chain Monte Carlo (MCMC) intended to analyze a set of stellar radial velocity data for evidence of extrasolar planets. For a more detailed overview of FMCMC in action see my article ``Extra-solar planets via Bayesian Fusion MCMC,'' in `Astrostatistical Challenges for the New Astronomy', J. M. Hilbe (ed), Springer, 2013, pp. 121-148. The current latest implementation in this Mathematica notebook employs a significantly improved control system (the stage 1 CS has been replaced by a more efficient stage 0).
The latest algorithm details can be found in Chapter 1 of, 'Supplement to Bayesian Logical Data Analysis for the Physical Sciences,' which is available from the Cambridge University website for my text book, "Bayesian Logical Data Analysis for the Physical Sciences," under the Resources section.
http://www.cambridge.org/pl/academic/subjects/statistics-probability/statistics-physical-sciences-and-engineering/bayesian-logical-data-analysis-physical-sciences-comparative-approach-mathematica-support

FMCMC runs multiple chains in parallel (typically 1-8) each with a different tempering parameter β. The number of chains (nβ values) can easily be modified as can the tempering levels. The default set of β values have worked well for a variety of data sets. Every npara interations (by default napar = 10),  proposals to swap the parameter values of a random pair of adjacent chains is considered. This prevents the FMCMC from becoming trapped in one probability peak region.

FMCMC joins BUGS (Bayesian inference Using Gibbs Sampling) and JAGS (Just Another Gibbs Sampler) as a powerful tool to facilitate Bayesian analysis.

1. SETUP

Remove["Global`*"]

Off[General::"spell"]

Off[General::"spell1"]

Needs["MultivariateStatistics`"]

Needs["ComputationalGeometry`"]

Needs["ErrorBarPlots`"]

2. LAUNCH SLAVES

LaunchKernels[]

3. DEFINE AND PLOT THE DATA WITH ERROR BARS

(*yvec1=y1vec*)

Raw measured y values.

yvec1={-3.5176371923433583`,-0.5103290583213148`,1.7172773377218313`,-0.4885056611547153`,2.881223100054913`,4.722706101035343`,-6.890446270665659`,1.662983759309704`,-4.237306928823851`,-0.36893877202014014`,-0.15082459692255012`,-0.7953391981200502`,-4.697044887902264`,5.476259880958446`,-1.5542407416654775`}

yscale=StandardDeviation[yvec1]

Rescale the measured y values.

yvec=(yvec1)/yscale

(*xvec1=x1vec*)

Raw measured x values.

xvec1={0.11731743852654608`,-2.5002885313531666`,-0.7019003999642486`,-4.222189623544026`,5.283218831913044`,6.5125815533467195`,-12.336002318962509`,8.392385690238529`,-4.523638026268021`,-0.6644059400772453`,-5.447685137122657`,-0.12896971469656648`,-4.304043503684108`,13.848263630630631`,-1.8534703757895508`}

xscale=StandardDeviation[xvec1]

Rescale the measured x values.

xvec=(xvec1)/xscale

sigmax=1.;sigmay=1.0;

Rescale the sigmax and sigmay measurement uncertainies.

σx=sigmax/xscale;σy=sigmay/yscale;nd=Length[xvec]

xerrvec=Table[σx,{nd}]

yerrvec=Table[σy,{nd}]

Rescale the mean and standard deviation of the Gaussian used to draw the raw true x values.

trueMu=0.0/xscale

trueTau=7/xscale

Here are the rescaled measurements and errors for plotting.

derr=Table[{{xvec[[i]],yvec[[i]]},ErrorBar[{-σx,σx},{-σy,σy}]},{i,1,nd}]

plxy1=ErrorListPlot[derr,Frame->True,PlotRange->{All,{-2.9,2.9}},Axes->False,PlotLabel->Style["",{"Helvetica-Bold",12}],PlotStyle->{Black,Thickness[0.003]},FrameLabel->{Style["x",{"Helvetica-Bold",16}],Style["y",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Here is the original raw straight line intercept and slope.

A=0.0;B=0.5;

Here is the raw true x values before rescaling.

xtvecRaw={0.27615947685146963`,-2.608790204578791`,0.34217263164449846`,-2.9651977758825065`,4.930101740012054`,6.081036426222702`,-12.763402274572924`,6.622446831000894`,-5.317762164032756`,0.4774641274465768`,-5.165408399772968`,-0.7010243149465465`,-3.091166685453353`,14.150327400268095`,-2.3204822994711396`};

Here are the points representing the rescaled true straight line.

xmin1=1.1*Min[xtvecRaw];xmax1=1.2*Max[xtvecRaw];
ymin1=A+B*xmin1
ymax1=A+B*xmax1
pnts1={{(xmin1)/xscale,(ymin1)/yscale},{(xmax1)/xscale,(ymax1)/yscale}}

trueSlope=(pnts1[[2,2]]-pnts1[[1,2]])/(pnts1[[2,1]]-pnts1[[1,1]])

Here is the rescaled slope which is labeled trueSlope for the purpose of comparing with inferences made about the slope from our analysis of the rescaled data set.

trueSlope=0.5*xscale/yscale

The rescaled intercept.

trueIntercept=0.0/yscale

trueSigma=1.0/yscale

plxyline=ListPlot[pnts1,Frame->True,PlotRange->All,Axes->False,Joined->True,PlotLabel->Style["",{"Helvetica-Bold",12}],PlotStyle->{Gray,Thickness[0.006]},FrameLabel->{Style["x",{"Helvetica-Bold",14}],Style["y",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

plxy1A=Show[{plxy1,plxyline},AspectRatio->0.5]

plxy1a=Show[{plxy1,plxyline},Graphics[{Text["(a)",{0,-2.4}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,DisplayFunction->Identity]

plxy1b=Show[{plxy1,plxyline},Graphics[{Text["(b)",{0,-2.4}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,DisplayFunction->Identity]

The rescaled true x values.

xtvec=(xtvecRaw)/xscale

The original raw true y values that included intrinsic scatter.

ytvecRaw={-0.1984716077695378`,-0.7339448858812833`,-0.1782696707577021`,-1.327950566427275`,4.383000563983801`,4.661730450696535`,-5.839836676699235`,2.1955529171902572`,-2.9267717741165544`,0.5350159424802654`,-0.9787935753116568`,-0.09070580260060801`,-3.8674209207625783`,7.421015629262308`,-0.782063516211839`};

The rescaled true y values.

ytvec=(ytvecRaw)/yscale

Plot of the rescaled true x & y values.

plxtyt=ListPlot[Transpose[{xtvec,ytvec}],Frame->True,PlotRange->{All,{-2.9,2.9}},Axes->False,PlotLabel->Style["",{"Helvetica-Bold",12}],PlotStyle->{Black,PointSize[0.015],Thickness[0.003]},FrameLabel->{Style["(Subscript[x, t])",{"Helvetica-Bold",16}],Style["(Subscript[y, t])",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

plxtyt1=Show[{plxtyt,plxyline},AspectRatio->0.5]

plxtyt1a=Show[{plxtyt,plxyline},Graphics[{Text["(a)",{0,-2.4}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,DisplayFunction->Identity]

plxtytxy=Show[GraphicsGrid[{{plxtyt1a,plxy1b}},Spacings->{Scaled[-0.25],Scaled[0.05]}]]

(*Export["c:/Users/Philip/Documents/greg13/MCMCregression/fig/CorYesStLineFitIntrinsicScatYesM7informGausXtPriorAnalyticX5fast22Aug13data.eps", plxtytxy,ImageSize->72*9]*)

4. INITIALIZE PROGRAM VALUES (depends on model)

npara=10; (* number of iterations in one parallel step *)

nthin=1*npara; (* store MCMC iterations at a thinning interval of nthin *)

n=50000*nthin; (* number of MCMC iterations *)

loglikeoffset=100.0; (* offset added to loglike calculation to keep the calculations within the machine precision range which is of importance in the Marginal Likelihood part of the calculations via the Nested Restricted Mont Carlo Method. *)

nβ=1;(* the number of choices of tempering parameter *)

βlis={1.0};

(*βlis={0.09,0.13,0.20,0.30,0.42,0.55,0.74,1.0};*)
(* a sample 8 tempering level set of β values *)

n2β=nβ; (* 2nd stage CS operates only on σall for top n2β chains *)

ns=1*npara; (* ns is number of iterations between swap proposals *)

(* Note:the product nβ\[Cross]ns/(swap acceptance rate) should be ≈ 100 to 200*)

np=5; (* np is the number of model parameters computational explored *)

ncs=5; (* ncs = number of model parameter proposal distributions to be varied *)

nl=np+6;(* item in output list containing prior1*(((like1))^β) (don't change this) *)

n2l=np+5;(* item in output list containing Log[like] (don't change this) *)

facc1=0.5; (* fraction of iterations correlated proposal sampler (proposePcor) used *)

swg=1; (* swg=1 if gene swap included, swg = 0 otherwise *)

smp=1;(* smp=0 or 1, smp=1 better multiple peak detection, reduced acceptance rate*)

5. STARTING PARAMETER VALUES (depends on model)

avst=Table[0,{np}]; (* initialize table of start parameter values *)
name=Table[0,{np}]; (* initialize table of parameter names *)
unit=Table[0,{np}]; (* initialize table of parameter units *)

On the next line user specifies the parameter labels as strings (within " ") the way they will appear in plots.

name={"α","β","σ","μ","τ"}

On the next line user specifies the parameter units as strings (within " ") the way they will appear in plots.

unit={" "," "," "," "," "}

On the next line user specifies the parameter start values.

avst={0.05,1.7,0.8,0.1,0.3}

On the next line user specifies the parameter types.

typeA={0,0,0,0,0};
typeB=Table[-2,{nd-1}];typeC=Table[-3,{nd-1}];
(*type=Join[typeA,typeB,typeC]*)
type=typeA; 
(* If type = -3, parameter is a true y coord η *)
(* If type = -2, parameter is a true x coord ξ *)
(* If type = -1, parameter is a period *)
(* type = 1, angular parameter, wrap around at 0 and 1 *)
(* type = 2, angular parameter, wrap around at 0 and 4π *)
(* type = 3, angular parameter, wrap around at 0 and 4π *)
(* type = 4, angular parameter, wrap around at -2π and 2π *)
(* type = 5, parameter is velocity amplitude (Subscript[K, 1]) *)
(* type = 0, all other parameters *)

6. PRIOR SEARCH RANGES (depends on model)

avmin=Table[0,{np}];(* initialize table of parameter lower bounds *)
avMax=Table[0,{np}];(* initialize table of parameter upper bounds *)

Next line provides a reminder of the parameter names (user don't change).

name

On the next line user specifies the lower bound on search range for each parameter.

avmin={-5.0,-99.0,0.0,-5.0,0.0}

On the next line user specifies the upper bound on search range for each parameter.

avMax={5.0,99.0,5.0,5.0,5.0}

7. PRIOR RANGE FOR BAYES FACTOR CALCULATION (depends on model)

bfmin=Table[0,{np}];(* Bayes factor parameter lower bounds *)
bfMax=Table[0,{np}];(* Bayes factor parameter upper bounds *)

Next line provides a reminder of the parameter names (user don't change).

name

On the next line user specifies the lower bound on each parameter for Bayes Factor calculation.

bfmin={-5.0,-99.0,0.0,-5.0,0.0}

On the next line user specifies the upper bound on each parameter for Bayes Factor calculation.

bfMax={5.0,99,5.0,5.0,5.0}

8. CONTROL SYSTEM PARAMETERS  (user don't change)

sw0=1;(* =1 when 0th stage control system is on and =0 when turned off.*)

sw5=1;(* =0 never execute 2rd stage CS, =1 execute 2nd stage CS *)

sw5a=0;(* =0 do not allow iterative repeats of 2rd stage CS, =1 allow iterative repeats *)

If[sw0>0,n0=5000,n0=1000]; (* duration of stage 0 CS in MCMC iterations. Note sw0=1 means stage 0 CS active *)

deltaN=5000; (* intervals in i used to decide on whether to rerun 2nd stage CS *)

sw8 =1; (* =0 never test, =1 execute test, by default set = 0 *)
(*** Test checks if start parameters with latest extra noise parameter are better ***)
(*** Only perform test once and only if sw8 > 0 and ncsys < n ***)
(*** Only used when initial parameter set input thought to be best ***)
(*** but extra noise parameter set to a large value to explore more widely ***)

n2β=nβ; (* 2nd stage CS operates only on σall for top n2β chains *)

ns=1*npara; (* ns is number of iterations between swap proposals *)

(* Note:the product nβ\[Cross]ns/(swap acceptance rate) should be ≈ 100 to 200*)

λ=0.25; (* target acceptance ratio *)

n3=1000;(* number of iterations per parameter in stage 2 *)

scmax=10.0; (* max factor to reduce/increase σ *)

alfm=0.6426+0.1507*Exp[-np/15](* Scaling exponent for calulation of accept rate for individual parameters 2nd stage CS *)

kalf1=0.87;(* Empirical factor for scaling exponent *)

kalf=If[nβ>1,kalf1*Table[1.0-(nβ-i)*0.2/(nβ-1),{i,1,nβ}],{kalf1}] (* Empirical factor for scaling exponent *)

accept=Table[(((λ))^(1/(np^(kalf[[kn]]*alfm)))),{kn,1,nβ},{kj,1,ncs}](* Calulation of accept rate for individual parameters in 2nd stage CS *)

(*accept=Table[(((λ))^(1/(ncs^((kalf^((βlis[[kn]])^(-0.3)))*alfm)))),{kn,1,nβ},{kj,1,ncs}]*)

accept[[1,1]]

9. STARTING SIGMAS  (don't normally need to change this)

During the post burn-in period the FMCMC program makes new parameter set proposals in two ways. For a fraction of the iterations the program employs an independent 'I' Gaussian proposal scheme (a separate Gaussian for each parameter). For the remaining fraction, specified by the progam value facc1, proposals are drawn using the 'C' proposal scheme which embeds the parameter correlations. Because the 'C' proposals reflect the parameter correlations, large jumps are possible allowing for more efficient movement in parameter space than can be achieved by the 'I' proposals alone. In item (1) above facc1 is set to a typical value of 0.5, so half the time we employ 'I' proposals and half the time 'C' proposals. At each iteration the choice of whether to use a 'C' proposal is determined by whether the value of a random number in the interval between 0 and 1 is lass than facc1. The combination of 'I' and 'C' proposals ensure that the whole parameter space can be reached and that the FMCMC chain is aperiodic. How is this achieved? 

During the burn-in phase a multi-stage adaptive control system automates the selection of efficient 'I' and 'C' proposal distributions. The control system is controlled by three switches, sw0, sw5, and facc1. In normal operation sw0=1, sw5=1 and facc1 = 0.5. After burn-in the 'I' and 'C' proposal distributions are not changed unless the adaptive control system detects that a new significantly higher probability peak in parameter space has emerged. If this happens the proposal distributions are updated and the burn-in phase extended. 

At the start, an initial vector of 'I' proposal sigma called, σst, is specified. In σst each proposal sigma value is set equal to 15% of the prior search range for that particular parameter. Since the prior search range is specified in step (3) above, this is automated by the following command:

σst=Table[If[type[[i]]==-1,0.15*(1./avmin[[1]]-1./avMax[[1]]),0.15*(avMax[[i]]-avmin[[i]])],{i,1,np}];

NB: a parameter with a type designation = -1 is a frequency so we need to convert the period values specified by avmin and avMax to frequencies.

σst becomes the starting sigma set for all the chains (specified by σall below) regardless of their tempering β value.

If[sw0==1,σall=Table[σst[[i]],{k,1,nβ},{i,1,np}]];

Another possible option is to fix the set of 'I' proposal sigma to desirable values (e.g., from a previous run) and only allow the 'C' proposal set to evolve during the burn-in. This is indicated by setting  sw0= 0,  sw5=0. If that is your intention, then remove the (* *) from the beginning and end of the next command and substitute your desired set for the nβ tempering chains by replacing the list in red brackets. Note, when a list is copied and pasted from a Mathematica output line, real numbers are given precisely as stored.

(*If[sw0==0&&sw5==0,σall={{0.0005390381230073959`,0.3459481022652124`,4 π,3.096968414314065`,2 π,0.0006085192595099706`,0.40487934750024235`,4 π,2.9649928523886966`,2 π,23.671932120052507`,8.727836408628711`},{0.0005782964831892351`,0.6013225678800871`,11.878416308529`,3.512175824917636`,2 π,0.0020058096409386285`,0.574609649209688`,4 π,5.321329808056671`,2 π,11.169708423072647`,5.6233735374939116`},{0.0005604873345521284`,0.498855611651304`,4 π,2.6527106195377543`,2 π,0.0009189134624128917`,0.46834126334979953`,4 π,2.230223614457575`,2 π,7.937755417489954`,4.004753274870652`},{0.000051653231543274674`,0.16770065159513692`,0.15429953863190327`,2.0421276745723795`,1.9454611963821122`,0.00018763138834972494`,0.20552109345197359`,0.16566749646159962`,2.741978621292261`,1.1188795279363213`,1.6663047611321666`,1.4698073555329132`},{0.0000376229347698429`,0.10864299048805837`,0.12296056576155409`,1.540937487379388`,1.2221235692158612`,0.00005529829151756749`,0.11968822157067664`,0.12960782188503808`,1.9442833298635613`,0.7623008636254842`,1.1306941111754456`,1.1832937278866873`},{0.00002490299062850884`,0.051920071008914895`,0.07827811543544246`,1.1758550244964348`,0.42096386277980763`,0.00003277276693020877`,0.05824792481719135`,0.08216436074754138`,1.3677908192047616`,0.3416807321074294`,0.7947630168596174`,0.6203984759377142`},{0.000021458164217817645`,0.0401403496342861`,0.07187975875817842`,0.8709083426923744`,0.4809539635153016`,0.00002673522030681429`,0.0460049440941334`,0.062380532931659224`,1.1071009288092482`,0.30637201781143375`,0.6370651104551156`,0.5554600564634931`},{0.000017140571815839838`,0.03387652697756123`,0.058023345242256386`,0.7494871289329318`,0.3780866273389161`,0.000021076322882366344`,0.034740243306850424`,0.05829845173728058`,0.923327754673528`,0.27167636162274206`,0.559863190870571`,0.453508273203633`}}]; *) 
(* initialize set of 'I' σ for all nβ parallel chains *)

10. GENE BOUNDARIES

ngene=1;

genst=Table[0,{ngene}]; (* initialize table of gene start boundaries *)

genend=Table[0,{ngene}]; (* initialize table of gene end boundaries *)

genst[[1]]= 1; (* start of gene 1 *)

genend[[1]]= np; (* start of gene 1 *)

11. NORMALIZATION FACTORS FOR MODEL PARAMETER PRIORS (depends on model)

mJlim=Mean[xerrvec];(* constant for modified scale invariant prior for τ *)

mJslim=2*Mean[yerrvec];(* constant for scale invariant prior for intrinsic scatter *)

c=-Log[(bfMax[[1]]-bfmin[[1]])*(bfMax[[2]]-bfmin[[2]])*Log[((mJslim+bfMax[[3]])/mJslim)]*(bfMax[[4]]-bfmin[[4]])*Log[((mJlim+bfMax[[5]])/mJlim)]]

12. DEFINE STEPS TO COMPUTE LIKELIHOOD FOR CURRENT MODEL (depends on model)

Specify a procedure for calculation of basic work = prior×likelihood for model. 
Uses the format ({{Module[{a, b, … }, proc] , a procedure with local variables a, b, …  }})

Evaluate constant term in the likelihood.

logfact=-nd*Log[2*π]

yx={yvec,xvec}

workP:=
Module[{a1,a2,a3,a4,a5,x,y},a1=avs[[1]];a2=avs[[2]];a3=avs[[3]];a4=avs[[4]];a5=avs[[5]];
ζ={a1+a2*a4,a4};yx-ζ;
vmat=Table[{{(a2^2)*(a5^2)+(a3^2)+((yerrvec[[i]])^2),a2*(a5^2)},{a2*(a5^2),(a5^2)+((xerrvec[[i]])^2)}},{i,1,nd}];
logprior=c -Log[a3+mJslim]-Log[a5+mJslim]; (* the prior for parameters *)
loglike=loglikeoffset+logfact-0.5*(∑_(i=1)^nd)(Log[Det[vmat[[i]]]])-0.5*(∑_(i=1)^nd)((Transpose[yx-ζ][[i]].Inverse[vmat[[i]]].Transpose[{Transpose[yx-ζ][[i]]}])[[1]]);
   (* likelihood fn for particular parameters *)  
ordP=FromDigits[Ordering[Extract[avs,Position[type,-1]]]];
pvec=Join[avs,{ordP,swpmv}] 
]

13. MODULE THAT PROPOSES THE NEW PARAMETER SET  (don't normally need to change this)

Propose new parameter values uses Switch and type[[ii]]to handle different parameter types. 
Update Period via frequency proposal.
Check if proposed values lie outside prior borders, if so reflect parameter or wrap angular parameters.

Parallel version that proposes a change to all of the avs parameter values by embedding this in a Table command.
Differs from the earlier propose (which was employed in a Do loop for non parallel version) in the ending because proposeP used within a Table or ParallelTable command.

proposeP:=
Module[{avmint,avMaxt},
Switch[type[[ii]],
0,
av0=RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-3,
av0=RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-2,
av0=RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-1,
If[smp>0&&i>ncsys&&RandomReal[]<0.25,avm1=1.0/RandomReal[NormalDistribution[1.0/values[[ii]],2*σalln[[ii]]]];avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]],avm1=1.0/RandomReal[NormalDistribution[1.0/values[[ii]],σalln[[ii]]]];avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]]],
-11,
If[smp>0&&i>ncsys&&RandomReal[]<0.25,avm1=1.0/RandomReal[NormalDistribution[1.0/values[[ii]],2*σalln[[ii]]]];avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]],avm1=1.0/RandomReal[NormalDistribution[1.0/values[[ii]],σalln[[ii]]]];avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]]],
1,
avs[[ii]]=Mod[RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]],1],
2,
avs[[ii]]=Mod[RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]],4π],
3,
avs[[ii]]=Mod[RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]],2*π],
4,
avs[[ii]]=Mod[RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]],4π,-(2π)],
5,
av5=RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]];avs[[ii]]=av5;avmint=avmin[[ii]];avMaxt=avMax[[ii]]*(((avmin[[ii-3]]/avs[[ii-3]]))^(1/3))*(((1-((avs[[ii-2]])^2)))^(-0.5));
Which[av5<avmint,avs[[ii]]=Mod[2*avmint-av5,(avMaxt-avmint),avmint],av5>avMaxt,avs[[ii]]=Mod[2*avMaxt-av5,(avMaxt-avmint),avmint]]
];avs[[ii]]
]

Parallel version that proposes a change to all avs parameter values based on the difference in two randomly selected parameter vectors previously stored in the correlated sample buffer.

proposePcor:=
Module[{avmint,avMaxt},
sampdiff=corscalek*(samp1-samp2);
Switch[type[[ii]],
0,
av0=values[[ii]]+sampdiff[[ii]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-3,
av0=values[[ii]]+sampdiff[[ii]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-2,
av0=values[[ii]]+sampdiff[[ii]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-1,
avm1=1.0/(1.0/values[[ii]]+sampdiff[[ii]]);avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]],
-11,
avm1=1.0/(1.0/values[[ii]]+sampdiff[[ii]]);avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]],
1,
samp1m=Mod[samp1[[ii]],1];samp2m=Mod[samp2[[ii]],1];sampd1=(samp1m-samp2m);
samp1mpi=Mod[(samp1m+0.5),1];samp2mpi=Mod[(samp2m+0.5),1];sampd1pi=(samp1mpi-samp2mpi);
sampdiff1=If[Abs[sampd1]≤Abs[sampd1pi],sampd1,sampd1pi];
avs[[ii]]=Mod[(values[[ii]]+corscalek*sampdiff1),1],
2,
samp1m=Mod[samp1[[ii]],2*π];samp2m=Mod[samp2[[ii]],2*π];sampd1=(samp1m-samp2m);
samp1mpi=Mod[(samp1m+π),2*π];samp2mpi=Mod[(samp2m+π),2*π];sampd1pi=(samp1mpi-samp2mpi);
sampdiff1=If[Abs[sampd1]≤Abs[sampd1pi],sampd1,sampd1pi];
avs[[ii]]=Mod[(values[[ii]]+corscalek*sampdiff1),4π],
3,
samp1m=Mod[samp1[[ii]],2*π];samp2m=Mod[samp2[[ii]],2*π];sampd1=(samp1m-samp2m);
samp1mpi=Mod[(samp1m+π),2*π];samp2mpi=Mod[(samp2m+π),2*π];sampd1pi=(samp1mpi-samp2mpi);
sampdiff1=If[Abs[sampd1]≤Abs[sampd1pi],sampd1,sampd1pi];
avs[[ii]]=Mod[(values[[ii]]+corscalek*sampdiff1),2*π],
4,
samp1m=Mod[samp1[[ii]],2*π];samp2m=Mod[samp2[[ii]],2*π];sampd1=(samp1m-samp2m);
samp1mpi=Mod[(samp1m+π),2*π];samp2mpi=Mod[(samp2m+π),2*π];sampd1pi=(samp1mpi-samp2mpi);
sampdiff1=If[Abs[sampd1]≤Abs[sampd1pi],sampd1,sampd1pi];
avs[[ii]]=Mod[(values[[ii]]+corscalek*sampdiff1),4π,-(2π)],
5,
av5=values[[ii]]+sampdiff[[ii]];avs[[ii]]=av5;avmint=avmin[[ii]];avMaxt=avMax[[ii]]*(((avmin[[ii-3]]/avs[[ii-3]]))^(1/3))*(((1-((avs[[ii-2]])^2)))^(-0.5));
Which[av5<avmint,avs[[ii]]=Mod[2*avmint-av5,(avMaxt-avmint),avmint],av5>avMaxt,avs[[ii]]=Mod[2*avMaxt-av5,(avMaxt-avmint),avmint]]
];avs[[ii]]
]

Parallel version that proposes a change to one of the avs parameter values.

proposeP1:=Module[{avmint,avMaxt},
Switch[type[[ii]],
0,
av0=RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-3,
av0=RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-2,
av0=RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]];avs[[ii]]=av0;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[av0<avmint,avs[[ii]]=Mod[2*avmint-av0,(avMaxt-avmint),avmint],av0>avMaxt,avs[[ii]]=Mod[2*avMaxt-av0,(avMaxt-avmint),avmint]],
-1,
If[smp>0&&i>ncsys&&RandomReal[]<0.25,avm1=1.0/RandomReal[NormalDistribution[1.0/values[[ii]],2*σalln[[ii]]]];avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]],avm1=1.0/RandomReal[NormalDistribution[1.0/values[[ii]],σalln[[ii]]]];avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]]],
-11,
If[smp>0&&i>ncsys&&RandomReal[]<0.25,avm1=1.0/RandomReal[NormalDistribution[1.0/values[[ii]],2*σalln[[ii]]]];avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]],avm1=1.0/RandomReal[NormalDistribution[1.0/values[[ii]],σalln[[ii]]]];avs[[ii]]=avm1;avmint=avmin[[ii]];avMaxt=avMax[[ii]];Which[avm1<avmint,avs[[ii]]=Mod[2*avmint-avm1,(avMaxt-avmint),avmint],avm1>avMaxt,avs[[ii]]=Mod[2*avMaxt-avm1,(avMaxt-avmint),avmint]]],
1,
avs[[ii]]=Mod[RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]],1],
2,
avs[[ii]]=Mod[RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]],4π],
3,
avs[[ii]]=Mod[RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]],2*π],
4,
avs[[ii]]=Mod[RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]],4π,-(2π)],
5,
av5=RandomReal[NormalDistribution[values[[ii]],σalln[[ii]]]];avs[[ii]]=av5;avmint=avmin[[ii]];avMaxt=avMax[[ii]]*(((avmin[[ii-3]]/avs[[ii-3]]))^(1/3))*(((1-((avs[[ii-2]])^2)))^(-0.5));
Which[av5<avmint,avs[[ii]]=Mod[2*avmint-av5,(avMaxt-avmint),avmint],av5>avMaxt,avs[[ii]]=Mod[2*avMaxt-av5,(avMaxt-avmint),avmint]]
];avs
]

14. MODULE THAT ACCEPTS THE NEW PARAMETER SET & correlation buffer values (don't change this)

acceptP:=Module[{},
pdfβi=logprior+β*loglike; (* factor needed to decide on tempering swap *)
pdfβiA=logprior+loglike;
testr1=pdfβi-logpdftempk[[nl]];
If[testr1>-uflow,If[testr1<uflow,r1=Exp[testr1],r1=1.0],r1=-1.0];  (* compute acceptance rule r1*)
u1=RandomReal[]; (* if uniform random ≤ r1 then accept proposal *)
If[u1≤r1,If[i>n0,nc2k=nc2k+1;nak=nak+1,If[sw0>0,nc2k=nc2k+1]];
pdfβiB1=Join[{β},pvec,{i},{loglike},{pdfβi}];
If[swpmv<1&&RandomReal[]<facc2,accdif=Table[If[type[[ia]]==-1||type[[ia]]==-11,1.0/avs[[ia]],avs[[ia]]],{ia,1,np}];accmov=Append[accmov,Join[accdif,{ordP}]]];
If[swpmv>0,nacork=nacork+1];
If[pdfβi>Last[maxlogpdfβik],maxlogpdfβik=pdfβiB1];If[pdfβiA>Last[maxlogpdfAlla],maxlogpdfAlla=Join[{β},pvec,{i},{loglike},{pdfβiA}]];
pdfβiB={{pdfβiA},{loglike},pdfβiB1,{nc2k},{nak},maxlogpdfβik,maxlogpdfAlla,accmov,{nacork},{npropPk}}
,pdfβiC=logpdftempk;pdfβiC[[n2l-1]]=i;pdfβiC[[n2l-2]]=swpmv-2;
pdfβiB={{pdfβiA},{loglike},pdfβiC,{nc2k},{nak},maxlogpdfβik,maxlogpdfAlla,accmov,{nacork},{npropPk}}
];pdfβiB
]

For use in second stage Control System use acceptP1.

acceptP1:=Module[{pdfβi,pdfβiA},
pdfβi=logprior+β*loglike; (* factor needed to decide on tempering swap *)
pdfβiA=logprior+loglike;
testr1=pdfβi-logpdftempk[[nl]];
If[testr1>-uflow,If[testr1<uflow,r1=Exp[testr1],r1=1.0],r1=-1.0]; (* compute acceptance rule r1*)
u1=RandomReal[]; (* if uniform random ≤ r1 then accept proposal *)
If[u1≤r1,nc2k=nc2k+1;
pdfβiB1=Join[{β},pvec,{i},{loglike},{pdfβi}];
If[pdfβi>Last[maxlogpdfβik],maxlogpdfβik=pdfβiB1];If[pdfβiA>Last[maxlogpdfAlla],maxlogpdfAlla=Join[{β},pvec,{i},{loglike},{pdfβiA}]];
pdfβiB={{pdfβiA},{loglike},pdfβiB1,{nc2k},maxlogpdfβik,maxlogpdfAlla}
,pdfβiC=logpdftempk;pdfβiC[[n2l-1]]=i;
pdfβiB={{pdfβiA},{loglike},pdfβiC,{nc2k},maxlogpdfβik,maxlogpdfAlla}
];pdfβiB
]

15. MODULE THAT ADJUSTS PROPOSAL σ IN THE 2nd Stage CS (user don't change this)

controlSysB:=
Module[{ip2,σallT,σallT2,σallTrem,σallT2T,σallT2TT},σallT=Transpose[σall];σallT2=Take[σallT,ncs];σallTrem=Take[σallT,-(np-ncs)];
σallT2T=Transpose[σallT2];
Do[σallT2T[[k]]=σallT2T[[k]]*(Sqrt[((N[acceptx[[k]]/n3]+0.01)/(accept[[k]]))(((1-accept[[k]]))/((1-N[acceptx[[k]]/n3]+0.01)))]),
{k,nβlow,nβ}];
σallT2TT=Transpose[σallT2T];σall=Transpose[Join[σallT2TT,σallTrem]];
(**************************************************************)
(* Control upper limit of σall to be ≤ parameter search range *)
Do[(* iteration over parameters *)
Do[Switch[type[[ip2]],
-3,If[σall[[k,ip2]]>avMax[[ip2]],σall[[k,ip2]]=avMax[[ip2]]],
-2,If[σall[[k,ip2]]>avMax[[ip2]],σall[[k,ip2]]=avMax[[ip2]]],
0,If[σall[[k,ip2]]>avMax[[ip2]],σall[[k,ip2]]=avMax[[ip2]]],
-1,If[σall[[k,ip2]]>(1./avmin[[ip2]]-1./avMax[[ip2]]),σall[[k,ip2]]=(1./avmin[[ip2]]-1./avMax[[ip2]])],
1,If[σall[[k,ip2]]>avMax[[ip2]],σall[[k,ip2]]=avMax[[ip2]]],
2,If[σall[[k,ip2]]>avMax[[ip2]],σall[[k,ip2]]=avMax[[ip2]]],
3,If[σall[[k,ip2]]>avMax[[ip2]],σall[[k,ip2]]=avMax[[ip2]]],
4,If[σall[[k,ip2]]>Abs[avMax[[ip2]]],σall[[k,ip2]]=Abs[avMax[[ip2]]]],
5,If[σall[[k,ip2]]>avMax[[ip2]],σall[[k,ip2]]=avMax[[ip2]]]
],{k,nβlow,nβ}
],(* complete β Do loop *)
{ip2,1,ncs}
];σall(* complete parameter Do loop *)
]

NotebookSave[]

16. STEP TO EXECUTE GREGORY'S FMCMC KERNEL (user don't change)

(* sw0 (1=on, 0=off) switch for stage 0 CS *)
(* sw5 (1=on, 0=off), if sw5=0 never executes stage 2 CS *)
(* sw6 (0 or 1) used in ParallelTable, enables calculation of i within each npara cycle *)

timerun=AbsoluteTiming[
accept1=Table[0,{i,1,nβ},{j,1,ncs}]; (* initialize storage for param. acceptance rates *)
accept2=Table[0,{i,1,nβ},{j,1,ncs}]; 
accept3=Table[0,{i,1,nβ},{j,1,ncs}];
corscale=Table[0.25,{nβ}]; (* initial value of corscale, a constant used in proposePcor *)
cs2cnt=1;(* initialize second stage CS counter *)
cs2cnt2=0;(* initialize counter to track how often 2nd stage CS is repeated a second time in a row *) 
f0=Table[If[i≤ncs,1.0,0.5*scmax],{k,1,nβ},{i,1,np}]; (*stg 0 CS σall[[i]] unchanged for i> ncs*)
facc2=0.5; (* fraction of accepted moves from proposeP added to accmov array *)
nnt=n/nthin;  (* no. of thinned iterations, n needs to be integer divisible by nthin *)
it=0;(* initialize counter for thinned iteration number *)
itmin=nnt;  (* determine "it" iteration # when movarray[[nβ]] has nmov entries *)
itsub=nnt;  (* determine "it" iteration # when movarray[[nβ]] has nmov entries post ncsys *)
kvals=Table[i,{i,1,nβ}];
logpdf=.;  (* remove any existing definition of file *)
logpdf=Table[0,{i,1,2},{j,1,nβ}]; (* initialize storage *) 
logpdfβi=.;  (* remove any existing definition of output file *)
logpdfβi=Table[0,{i,1,nnt},{j,1,nβ}]; (* initialize storage *)
maxlogpdfβi=Table[0,{j,1,nβ}]; (* initialize storage *)
maxlogpdfAllold={};(* record previous values of maxlogpdfAll *)
nmov=300; (* size of buffer for recently accepted parameter moves *)movarray=Table[Table[Table[0,{np+1}],{nmov}],{nβ}];(* an array of successful moves *)
na=Table[0,{nβ}];  (* initialize acceptance array *)
nacatot=Table[0,{nβ}];(* initialize counter for parameter sets added to movarray[[nβ]] *)
naccsub=Table[0,{nβ}];(* counter for parameter sets added to movarray[[nβ]] post ncsys *)
nacor=Table[0,{nβ}]; (* initialize acceptance array for correlated move sampler *)
nacorsub=Table[0,{nβ}]; (* initialize nacorsub array, reset =nacor each time ncsys reset *)
nagain=0; (* initial value of integer regarding repeats of stage 2 CS, =0 means no *)
nas=0; (* initiate an integer used in computing actual swap rate *)
If[sw0==0&&sw5==0,ncsys=n0,ncsys=n]; (* otherwise ncsys=iteration # when control system turned off *)
ncsys2=n;
ncsysold={}; (* records previous values of ncsys=iteration # *)
nc1=1;(* maj/minor cycle counter initially set to 1 *)
nc2=Table[0,{nβ}];(* # of prop. acceptances in each subcycle *)
ncsysacorsubscaleAll={};(* record previous values of corscale *)
npropP=Table[0,{nβ}]; (* initialize array for total number of correlated proposals made *)
npropPsub=Table[0,{nβ}]; (* initialize array for correlated proposals current start *)
pdfβiAmax2old={};(* record previous best logprior+loglike in β=1 chain in 2nd stage CS *)
(* specify some constants & switches used in the control system *) 
sw=Table[0,{nβ}];(* 0 in major cycles, +1 in minor cycles if er[[k]]<0 and -1 if er[[k]]>0 *)
sw4=0;(* =0 at start to skip over stage 2 CS, changed to =1 when stage 0 or 1 CS complete*)
swap={};  (* initialize diagnostic array for swaps *)
uflow=Log[$MaxNumber/10.0];(* set number to avoid underflow messages *) 
σallold={};(* record previous best σall *)
σ0monitor={};  (* initialize monitor of parameter proposal σ's for stage 0  *)
σmonitor={{1,σall}};  (* initialize monitor of parameter proposal σ's  *)
DistributeDefinitions["Global`"];
(* Evaluate the probability for the i=0 sample for each tempering level*)
Do[(* iteration over β *) β=βlis[[k]];avs=avst;
 (* initiate parameter values *)
(* ******************** Basic working *********************** *)
workP;
(* ********************************************************* *)
pdf1= loglike;
pdfβ1= logprior+β*loglike; (* factor needed to decide on tempering swap *)
logpdf[[1,k]]=Join[{β},pvec,{0,pdf1,pdfβ1}],{k,1,nβ} 
];  (* close first Do loop iterating over β *)
maxlogpdfβi=logpdf[[1]];
maxlogpdfAll=logpdf[[1,nβ]];(* initialize this array in each processor (chain) *)
Do[ (* start calculation for ith iteration beginning at i = 1 *)
iV=i+1-npara;sw6=Table[0,{nβ}];logpdftemp=logpdf[[1]];
With[{iV=iV,logpdftemp=logpdftemp,σall=σall,maxlogpdfβi=maxlogpdfβi,maxlogpdfAll=maxlogpdfAll,nc2=nc2,na=na,nacor=nacor,npropP=npropP,sw6=sw6,movarray=movarray,corscale=corscale,nacatot=nacatot},
workvals=
ParallelTable[β=βlis[[k]];sw6k=sw6[[k]];maxlogpdfAlla=maxlogpdfAll;movarrayk=movarray[[k]];nacatotk=nacatot[[k]];nacatotk1=0;npropPk=npropP[[k]];corscalek=corscale[[k]];accmov={};
Do[
If[RandomReal[]<facc1&&nacatotk≥nmov, (* If yes go the proposePcor route *)
If[sw6k<1,i=iV;logpdftempk=logpdftemp[[k]];values=Take[logpdftempk,{2,np+1}];avs=values;σalln=σall[[k]];maxlogpdfβik=maxlogpdfβi[[k]];nc2k=nc2[[k]];nak=na[[k]];nacork=nacor[[k]];sw6k=sw6k+1;
ordP=Take[logpdftempk,{nl-4}][[1]];movordP=Transpose[movarrayk][[np+1]];nordP=Flatten[Position[movordP,ordP]];movarrayordP=Take[movarrayk[[nordP]]];(*movarrayordP=Select[movarrayk,#[[np+1]]==ordP&];*)
nmovordP=Length[movarrayordP];
If[nmovordP≥5,samp1=movarrayordP[[RandomInteger[{1,nmovordP}]]];samp2=movarrayordP[[RandomInteger[{1,nmovordP}]]];nacatotk1=nacatotk,nacatotk1=0
],
i=i+1;logpdftempk=pdfβiB[[3]];values=Take[logpdftempk,{2,np+1}];avs=values;ordP=Take[logpdftempk,{nl-4}][[1]];movordP=Transpose[movarrayk][[np+1]];nordP=Flatten[Position[movordP,ordP]];movarrayordP=Take[movarrayk[[nordP]]];(*movarrayordP=Select[movarrayk,#[[np+1]]==ordP&];*)
nmovordP=Length[movarrayordP];
If[nmovordP≥5,samp1=movarrayordP[[RandomInteger[{1,nmovordP}]]];samp2=movarrayordP[[RandomInteger[{1,nmovordP}]]];nacatotk1=nacatotk,nacatotk1=0]
];
If[nacatotk1≥nmov,Do[proposePcor,{ii,1,np}];swpmv=1;npropPk=npropPk+1,Do[proposeP,{ii,1,np}];swpmv=0];workP;acceptP,
If[sw6k<1,i=iV;logpdftempk=logpdftemp[[k]];values=Take[logpdftempk,{2,np+1}];avs=values;σalln=σall[[k]];maxlogpdfβik=maxlogpdfβi[[k]];nc2k=nc2[[k]];nak=na[[k]];nacork=nacor[[k]];sw6k=sw6k+1,
i=i+1;logpdftempk=pdfβiB[[3]];values=Take[logpdftempk,{2,np+1}];avs=values
];
Do[proposeP,{ii,1,np}];swpmv=0;workP;acceptP],
{isub,1,npara}
];pdfβiB,
{k,1,nβ} 
];
];
tworkvals=Transpose[workvals];logpdf[[2]]=tworkvals[[3]];nc2=Flatten[Transpose[tworkvals[[4]]]];na=Flatten[Transpose[tworkvals[[5]]]];nacor=Flatten[Transpose[tworkvals[[9]]]];npropP=Flatten[Transpose[tworkvals[[10]]]];maxlogpdfβi=tworkvals[[6]];lastworkvals=tworkvals[[7]];
If[facc1>0,accarray=tworkvals[[8]];naca=Table[Length[accarray[[k]]],{k,1,nβ}];nacatot=nacatot+naca;If[nacatot[[nβ]]≥nmov&&it<itmin,itmin=it];If[i>ncsys,naccsub=naccsub+naca]];(* naccsub=Table[0,{nβ}] each time ncsys reset *)
dpropP=npropP-npropPsub;
If[naccsub[[nβ]]≥nmov&&Min[dpropP]≥50&&it<itsub,itsub=it;
acorsub=(nacor-nacorsub)/dpropP;If[Min[acorsub]>0.08,corscale=corscale*((((((acorsub+0.01))/0.25)*(0.75/((1-acorsub+0.01)))))^0.25)];ncsysacorsubscale={ncsys,ncsys2,itsub*nthin,acorsub,corscale};ncsysacorsubscaleAll=Append[ncsysacorsubscaleAll,ncsysacorsubscale];
If[acorsub[[nβ]]>0.28||acorsub[[nβ]]<0.2,naccsub=Table[0,{nβ}];npropPsub=npropP;nacorsub=nacor;ncsys2=i;itsub=nnt]
];
acceptvec=Sort[workvals,OrderedQ[{#2,#1}]&];
kbestvec=Sort[lastworkvals,#2[[nl]]<#1[[nl]]&][[1]];
pdfβiAmax=kbestvec[[nl]];
If[pdfβiAmax>maxlogpdfAll[[nl]],maxlogpdfAllold=Append[maxlogpdfAllold,maxlogpdfAll];maxlogpdfAll=kbestvec];If[facc1>0,Do[If[Length[accarray[[k]]]>0,movarray[[k]]=Take[Join[movarray[[k]],accarray[[k]]],-nmov]],{k,1,nβ}] 
];
If[i>n0,Goto[jump],
If[sw0>0,
Which[i==IntegerPart[n0/5],σ0monitor=Append[σ0monitor,{nc2,σall}];σall=f0*2.0*σall/scmax;nc2=Table[0,{nβ}],i==IntegerPart[2*n0/5],σ0monitor=Append[σ0monitor,{nc2,σall}];σall=f0*2.0*σall/scmax;nc2=Table[0,{nβ}],i==IntegerPart[3*n0/5],σ0monitor=Append[σ0monitor,{nc2,σall}];σall=f0*2.0*σall/scmax;nc2=Table[0,{nβ}],i==IntegerPart[4*n0/5],σ0monitor=Append[σ0monitor,{nc2,σall}];σall=f0*2.0*σall/scmax;nc2=Table[0,{nβ}],i==n0,σ0monitor=Append[σ0monitor,{nc2,σall}];nc2=Table[0,{nβ}];sw4=1;tσ0mon=Transpose[σ0monitor];ttσ0mon1=Transpose[tσ0mon[[1]]];ttσ0mon2=Transpose[tσ0mon[[2]]];near=Table[Nearest[ttσ0mon1[[k]],IntegerPart[n0/20]],{k,1,nβ}];posnear=Flatten[Table[Last[Position[ttσ0mon1[[k]],near[[k,1]]]],{k,1,nβ}]];σall=Table[Flatten[Take[ttσ0mon2[[k]],{posnear[[k]]}]],{k,1,nβ}];σmonitor=Append[σmonitor,{n0,σall}]
]
]
];
Label[jump];
(* Consider a swap between 2 tempering levels if i an integer multiple of ns *)
If[nβ<2,Goto[genecrossover]];
If[IntegerQ[i/ns],(* next pick the 2 levels k1 & k2 to swap *)
(* following command chooses k1 from 1 to nβ-1 and sets k2=k1-1 *)
k1=RandomInteger[{1,(nβ-1)}];
k2=k1+1;
(* evaluate swap acceptance rule*)
loglike1=logpdf[[2,k1,n2l]];
loglike2=logpdf[[2,k2,n2l]];
β1=logpdf[[2,k1,1]];β2=logpdf[[2,k2,1]];
logprior1=logpdf[[2,k1,nl]]-β1*loglike1;
logprior2=logpdf[[2,k2,nl]]-β2*loglike2;
p12=logprior1+β2*loglike1;
p21=logprior2+β1*loglike2;
r=Exp[(β2-β1)*(loglike1-loglike2)];  
u2=RandomReal[];
If[u2≤r (* if uniform random ≤ r then accept swap *),nas=nas+1;swap=Append[swap,{k1+k2,1,i}];
temp=logpdf[[2,k1]];
Do[logpdf[[2,k1,ji]]=logpdf[[2,k2,ji]],{ji,2,np+2}];
logpdf[[2,k1,n2l]]=loglike2;logpdf[[2,k1,nl]]=p21;
Do[logpdf[[2,k2,ji]]=temp[[ji]],{ji,2,np+2}];
logpdf[[2,k2,n2l]]=loglike1;logpdf[[2,k2,nl]]=p12,
swap=Append[swap,{k1+k2,0,i}]]  (* complete accept If swap statement *)
] ; (* complete first If swap statement *)
(* Gene swap/crossover between top chain(s) and values in maxlogpdfAll *)
Label[genecrossover];
If[swg>0&&IntegerQ[i/(2*nthin)],
If[nβ<2,
If[ngene<2,
topvalues=acceptvec[[1]];avtop=Take[topvalues[[3]],{2,np+1}];avcurmax=Take[maxlogpdfAll,{2,np+1}];
workvals=
Table[βnow=topvalues[[3,1]];
nbreak=RandomInteger[{genst[[1]]+1,genend[[1]]-1}];
If[RandomReal[]<0.5,
Do[avs[[ii]]=avcurmax[[ii]],{ii,1,nbreak}];
Do[avs[[ii]]=avtop[[ii]],{ii,nbreak+1,np}],
Do[avs[[ii]]=avtop[[ii]],{ii,1,nbreak}];
Do[avs[[ii]]=avcurmax[[ii]],{ii,nbreak+1,np}]
];
workP;pdfβiA=logprior+loglike; 
pdfβiB={{pdfβiA},Join[{1+βnow},pvec,{i},{loglike},{pdfβiA}]},
{kg,1,ngene} 
],(* end of ngene=1 & nβ=1 option *)
topvalues=acceptvec[[1]];avtop=Take[topvalues[[3]],{2,np+1}];avcurmax=Take[maxlogpdfAll,{2,np+1}];
With[{i=i,genst=genst,genend=genend,avcurmax=avcurmax,avtop=avtop,topvalues=topvalues},
workvals=
ParallelTable[
avs=avcurmax;βnow=topvalues[[3,1]];
Do[avs[[ii]]=avtop[[ii]],{ii,genst[[kg]],genend[[kg]]}];
workP;pdfβiA=logprior+loglike; 
pdfβiB={{pdfβiA},Join[{1+βnow},pvec,{i},{loglike},{pdfβiA}]},
{kg,1,ngene} 
]
]
],(* end of ngene>1 & nβ=1 option *)
If[ngene<2,
topvalues=acceptvec[[1]];avtop=Take[topvalues[[3]],{2,np+1}];
secvalues=acceptvec[[2]];avsec=Take[secvalues[[3]],{2,np+1}];
avcurmax=Take[maxlogpdfAll,{2,np+1}];
workvals=
Table[
Which[kg≤ngene,βnow=topvalues[[3,1]];
nbreak=RandomInteger[{genst[[1]]+1,genend[[1]]-1}];
If[RandomReal[]<0.5,
Do[avs[[ii]]=avcurmax[[ii]],{ii,1,nbreak}];
Do[avs[[ii]]=avtop[[ii]],{ii,nbreak+1,np}],
Do[avs[[ii]]=avtop[[ii]],{ii,1,nbreak}];
Do[avs[[ii]]=avcurmax[[ii]],{ii,nbreak+1,np}]
];
workP;pdfβiA=logprior+loglike;
pdfβiB={{pdfβiA},Join[{2+βnow},pvec,{i},{loglike},{pdfβiA}]},
kg≤2*ngene,βnow=secvalues[[3,1]];
nbreak=RandomInteger[{genst[[1]]+1,genend[[1]]-1}];
If[RandomReal[]<0.5,
Do[avs[[ii]]=avcurmax[[ii]],{ii,1,nbreak}];
Do[avs[[ii]]=avsec[[ii]],{ii,nbreak+1,np}],
Do[avs[[ii]]=avsec[[ii]],{ii,1,nbreak}];
Do[avs[[ii]]=avcurmax[[ii]],{ii,nbreak+1,np}]
];
workP;pdfβiA=logprior+loglike;
pdfβiB={{pdfβiA},Join[{2+βnow},pvec,{i},{loglike},{pdfβiA}]}
],
{kg,1,2*ngene}
], (* end of ngene=1 & nβ>1 option *)
topvalues=acceptvec[[1]];avtop=Take[topvalues[[3]],{2,np+1}];
secvalues=acceptvec[[2]];avsec=Take[secvalues[[3]],{2,np+1}];
avcurmax=Take[maxlogpdfAll,{2,np+1}];
With[{i=i,genst=genst,genend=genend,avcurmax=avcurmax,avtop=avtop,
topvalues=topvalues,avsec=avsec,secvalues=secvalues},
workvals=
ParallelTable[
If[kg≤ngene,avs=avcurmax;βnow=topvalues[[3,1]];
Do[avs[[ii]]=avtop[[ii]],{ii,genst[[kg]],genend[[kg]]}];
workP;pdfβiA=logprior+loglike;
pdfβiB={{pdfβiA},Join[{1+βnow},pvec,{i},{loglike},{pdfβiA}]},
avs=avcurmax;βnow=secvalues[[3,1]];kgg=kg-ngene;
Do[avs[[ii]]=avsec[[ii]],{ii,genst[[kgg]],genend[[kgg]]}];
workP;pdfβiA=logprior+loglike;
pdfβiB={{pdfβiA},Join[{2+βnow},pvec,{i},{loglike},{pdfβiA}]}
],
{kg,1,2*ngene}
]
]
]
];
acceptvec1=Sort[workvals,OrderedQ[{#2,#1}]&];
kbestvec=acceptvec1[[1]];pdfβiAmax=kbestvec[[1,1]];If[pdfβiAmax>Last[maxlogpdfAll],maxlogpdfAllold=Append[maxlogpdfAllold,maxlogpdfAll];
maxlogpdfAll=kbestvec[[2]]];
];(* end of gene swap/crossover If statement *)
If[sw5<1||sw4<1,(* If sw4==0 skip 2nd stage control system *)
If[IntegerQ[i/nthin],it=it+1;logpdfβi[[it]]=logpdf[[2]]];logpdf[[1]]=logpdf[[2]],
(************************************************************************)
(*** Control system 2nd stage, acceptance for individual parameters ***)
(************************************************************************)
If[IntegerQ[i/nthin],it=it+1;logpdfβi[[it]]=logpdf[[2]]];
(* ***************** Find best chain values ***************** *)
bestvec=Take[maxlogpdfAll,{2,np+1}];
swpmv=0;
Do[(* iteration over β *) β=βlis[[k]];avs=bestvec;
 (* initiate parameter values *)
(* ******************** Basic working ********************** *)
workP;
(* ********************************************************* *)
pdfj= loglike;
pdfβj= logprior+β*loglike; 
If[pdfβj>Last[maxlogpdfβi[[k]]],
maxlogpdfβi[[k]]=Join[{β},pvec,{1,pdfj,pdfβj}]],
{k,1,nβ} 
];  (* close first Do loop iterating over β *)
(* *********************************************************** *)
logpdf[[1]]=maxlogpdfβi;
logpdfinow2=maxlogpdfβi;
(************************************************************************)
(*** Now test if start parameters with current extra noise parameter are better ***)
(*** Only perform test once and only if sw8 > 0 and ncsys < n ***)
(*** Only used when initial parameter set input thought to be best ***)
(*** but extra noise parameter set to a large value to explore more widely ***)
If[sw8>0,
If[ncsys<n,
Do[(* iteration over β *) β=βlis[[k]];avs=avst;avs[[np]]=maxlogpdfβi[[k]][[np+1]];
 (* initiate parameter values *)
(* ******************** Basic working ********************** *)
workP;
(* ********************************************************* *)
pdfj= loglike;
pdfβj= logprior+β*loglike; 
If[pdfβj>Last[maxlogpdfβi[[k]]],
maxlogpdfβi[[k]]=Join[{β},pvec,{1,pdfj,pdfβj}]],
{k,1,nβ} 
];  (* close first Do loop iterating over β *)
(* *********************************************************** *)
logpdf[[1]]=maxlogpdfβi;
logpdfinow2=maxlogpdfβi;
sw8=0;
]
];
nβlow=nβ-n2β+1;
k2vals=Table[iip,{iip,nβlow,nβ}];
Do[(* iteration over parameters *)
n2c2=Table[0,{nβ}];(* # of para. prop. acceptances in each subcycle *)
iiV=ip2;sw7=Table[0,{nβ}];
With[{i=i,n3=n3,iiV=iiV,logpdfinow2=logpdfinow2,σall=σall,maxlogpdfβi=maxlogpdfβi,maxlogpdfAll=maxlogpdfAll,n2c2=n2c2,nβlow=nβlow,sw7=sw7,swpmv=swpmv},
workvals=
ParallelTable[β=βlis[[k]];sw7k=sw7[[k]];maxlogpdfAlla=maxlogpdfAll;
Do[If[sw7k<1,ii=iiV;logpdftempk=logpdfinow2[[k]];values=Take[logpdftempk,{2,np+1}];avs=values;σalln=σall[[k]];maxlogpdfβik=maxlogpdfβi[[k]];nc2k=n2c2[[k]];sw7k=sw7k+1,
logpdftempk=pdfβiB[[3]];values=Take[logpdftempk,{2,np+1}];avs=values];
proposeP1;workP;acceptP1,
{n3}
];pdfβiB,
{k,nβlow,nβ} 
]
];
tworkvals=Transpose[workvals];n2c2=Flatten[Transpose[tworkvals[[4]]]];maxlogpdfβi=tworkvals[[5]];lastworkvals=tworkvals[[6]];
kbestvec=Sort[lastworkvals,#2[[nl]]<#1[[nl]]&][[1]];
pdfβiAmax2=kbestvec[[nl]];
If[pdfβiAmax2>maxlogpdfAll[[nl]],maxlogpdfAllold=Append[maxlogpdfAllold,maxlogpdfAll];maxlogpdfAll=kbestvec];
Do[accept1[[k,ip2]]=n2c2[[k]],{k,1,nβ}]  
,{ip2,1,ncs}
];(* complete parameter Do loop *)
acceptx=accept1;
σold1=σall;
controlSysB;
σmonitor=Append[σmonitor,{i+100+(cs2cnt-1)*200,σall}];
(**************************************************************)
Do[(* iteration over parameters *)
n2c2=Table[0,{nβ}];(* # of para. prop. acceptances in each subcycle *)
iiV=ip2;sw7=Table[0,{nβ}];
With[{i=i,n3=n3,iiV=iiV,logpdfinow2=logpdfinow2,σall=σall,maxlogpdfβi=maxlogpdfβi,maxlogpdfAll=maxlogpdfAll,n2c2=n2c2,nβlow=nβlow,sw7=sw7,swpmv=swpmv},
workvals=
ParallelTable[β=βlis[[k]];sw7k=sw7[[k]];maxlogpdfAlla=maxlogpdfAll;
Do[If[sw7k<1,ii=iiV;logpdftempk=logpdfinow2[[k]];values=Take[logpdftempk,{2,np+1}];avs=values;σalln=σall[[k]];maxlogpdfβik=maxlogpdfβi[[k]];nc2k=n2c2[[k]];sw7k=sw7k+1,
logpdftempk=pdfβiB[[3]];values=Take[logpdftempk,{2,np+1}];avs=values];
proposeP1;workP;acceptP1,
{n3}
];pdfβiB,
{k,nβlow,nβ} 
]
];
tworkvals=Transpose[workvals];n2c2=Flatten[Transpose[tworkvals[[4]]]];maxlogpdfβi=tworkvals[[5]];lastworkvals=tworkvals[[6]];
kbestvec=Sort[lastworkvals,#2[[nl]]<#1[[nl]]&][[1]];
pdfβiAmax2=kbestvec[[nl]];
If[pdfβiAmax2>maxlogpdfAll[[nl]],maxlogpdfAllold=Append[maxlogpdfAllold,maxlogpdfAll];maxlogpdfAll=kbestvec];
Do[accept2[[k,ip2]]=n2c2[[k]],{k,1,nβ}] 
,{ip2,1,ncs}
];(* complete parameter Do loop *)
acceptx=accept2;
σold2=σall;
controlSysB;
σmonitor=Append[σmonitor,{i+200+(cs2cnt-1)*200,σall}];
(**************************************************************)
If[sw5a>0&&i>(n0+deltaN),
Do[(* iteration over parameters *)
n2c2=Table[0,{nβ}];(* # of para. prop. acceptances in each subcycle *)
iiV=ip2;sw7=Table[0,{nβ}];
With[{i=i,n3=n3,iiV=iiV,logpdfinow2=logpdfinow2,σall=σall,maxlogpdfβi=maxlogpdfβi,maxlogpdfAll=maxlogpdfAll,n2c2=n2c2,nβlow=nβlow,sw7=sw7,swpmv=swpmv},
workvals=ParallelTable[β=βlis[[k]];sw7k=sw7[[k]];maxlogpdfAlla=maxlogpdfAll;
Do[If[sw7k<1,ii=iiV;logpdftempk=logpdfinow2[[k]];values=Take[logpdftempk,{2,np+1}];avs=values;σalln=σall[[k]];maxlogpdfβik=maxlogpdfβi[[k]];nc2k=n2c2[[k]];sw7k=sw7k+1,
logpdftempk=pdfβiB[[3]];values=Take[logpdftempk,{2,np+1}];avs=values];
proposeP1;workP;acceptP1,
{n3}
];pdfβiB,
{k,nβlow,nβ} 
]
];
tworkvals=Transpose[workvals];n2c2=Flatten[Transpose[tworkvals[[4]]]];maxlogpdfβi=tworkvals[[5]];lastworkvals=tworkvals[[6]];
kbestvec=Sort[lastworkvals,#2[[nl]]<#1[[nl]]&][[1]];
pdfβiAmax2=kbestvec[[nl]];
If[pdfβiAmax2>maxlogpdfAll[[nl]],maxlogpdfAllold=Append[maxlogpdfAllold,maxlogpdfAll];maxlogpdfAll=kbestvec];
Do[accept3[[k,ip2]]=n2c2[[k]],{k,1,nβ}]
,{ip2,1,ncs}
];(* complete parameter Do loop *)
If[Min[N[accept3[[nβ]]/n3]]<0.3||Max[N[accept3[[nβ]]/n3]]>0.95,nagain=1,nagain=0]
];(* complete If statement for accept3 calculation *)
If[cs2cnt<2&&nagain>0,sw4=1;cs2cnt=cs2cnt+1;cs2cnt2=cs2cnt2+1,sw4=0;cs2cnt=1;nacorsub=nacor;npropPsub=npropP];
If[ncsys==n,ncsys=i;na=Table[0,{nβ}];naccsub=Table[0,{nβ}];nacorsub=nacor;npropPsub=npropP;itsub=nnt];
]; (* complete If statement for stage 2 CS *)
If[sw5>0&&i>ncsys&&IntegerQ[(i-ncsys)/deltaN]&&(Last[maxlogpdfAll]-pdfβiAmax2)>3.0,sw4=1;ncsysold=Append[ncsysold,ncsys];σallold=Append[σallold,σall];pdfβiAmax2old=Append[pdfβiAmax2old,pdfβiAmax2];na=Table[0,{nβ}];ncsys=i;ncsys2=ncsys;naccsub=Table[0,{nβ}];nacorsub=nacor;npropPsub=npropP;itsub=nnt,
];(* complete If statement *)
,{i,npara,n,npara} (* complete iteration Do loop *)
];
acceptratio=If[ncsys<n,N[na/(n-ncsys)],N[na/n]]  (* evaluate acceptance ratio *)
]

(* If facc1 >0, correlated proposal sampler employed and nburn = itsub, otherwise nburn = ncsys/nthin *)

The line above gives the execution time and average acceptance rate for each chain after the control system was turned off. The control system switched off at an iteration number specified by nburn. If facc1 > 0, correlated proposal sampler employed and nburn = itsub, otherwise nburn = ncsys/nthin

itsub

ncsys

ncsysold

na

nacor

nacatot

naccsub

dpropP

If[facc1>0&&naccsub[[nβ]]≥nmov&&Min[dpropP]≥50,acorsub=N[(nacor-nacorsub)/dpropP]]

corscale

If[facc1>0,N[ncsysacorsubscaleAll]]

itmin

bestvec

maxlogpdfAll

Length[maxlogpdfAllold]

maxlogpdfAllold

maxlogpdfβi

logpdfinow2

pdfβiAmax2old

MemoryInUse[]

n2c2

na

σall

σmonitor

timerun

The line above gives the execution time and average acceptance rate for each chain after the control system was turned off.

At each iteration a number of parameter values and diagnostics are recorded. The set which achieves the highest value 
of  prior x (likelihood) is called maxlogpdfAll regardless of the β value of the chain it occured in. Because of the included genetic crossover operation this will be a reasonable proxy for the MAP value. See Section 23 below for explanation of all terms except the first item. The first item is only relevant to the program designer.

maxlogpdfAll

(*maxlogpdfAll={2.`,-0.10657370589745688`,0.8770792896201273`,0.3772879923057287`,-0.019089508982076026`,0.9332873854466412`,0,0,202280,68.06793425148757`,55.73376006298533`};*)

17. Output MCMC results

To output the whole MCMC data set all tempering levels.

(*ng=Length[logpdfβi];stmp=OpenWrite["c:/Users/Philip/Gregory3planHD11964ParallelAllM7.dat",PageWidth->Infinity];
Do[Write[stmp,logpdfβi[[i]]],{i,1,ng}];Close[stmp]*)

tr=Transpose[logpdfβi];tr5=tr[[nβ]];ng=Length[tr5]

To output only the β = 1 tempering level results (need to remove the (* *) at the ends of lines).

tr=Transpose[logpdfβi];tr5=tr[[nβ]];ng=Length[tr5];stmp=OpenWrite["c:/Users/Philip/Documents/greg13/MCMCregression/regressLineAnalytic_Public_14May16M8ACER.dat",PageWidth->Infinity];
Do[Write[stmp,tr5[[i]]],{i,1,ng}];Close[stmp]

NotebookSave[]

18. Here is how to input results (need to remove the (* *) at the ends of lines).

dataMC=ReadList["c:/Users/Philip/Documents/greg13/MCMCregression/regressLineAnalytic_Public_14May16M8ACER.dat"];nnt=Length[dataMC];n=nthin*nnt;ng=nnt

dataMC[[5000]]

Next if only β = 1 tempering level input, populate all tempering levels with this so analysis algorithms will work.

logpdfβi=Transpose[Table[dataMC,{k,1,nβ}]];

19. Examine the automatic anealing of the intrinsic scatter parameter σ

tmaxlogpdfAll=Transpose[maxlogpdfAllold];lisP1=Join[{{-0.01,avst[[3]]}},Transpose[{tmaxlogpdfAll[[n2l-1]]*(10^(-5)),tmaxlogpdfAll[[4]]}],{{maxlogpdfAll[[n2l-1]]*(10^(-5)),maxlogpdfAll[[4]]}}];
plPsa=ListLogPlot[{lisP1},Frame->True,PlotRange->{{-0.01,0.5},All},PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iterations  (\[Cross] (10^5))",{"Helvetica-Bold",14}],Style["Intrinsic scatter",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

20. MCMC performance diagnostics & parameter proposal σ's

The next command gives the acceptance rate for each chain (tempering level).

acceptratio

The MCMC iteration number (unthinned) at which control system was turned off for the last time.

ncsys

The output from the next command gives the swap acceptance rate between adjacent tempering levels. First item is the swap rate between the two with lowest β tempering levels, etc.

If[nβ>1&&n>ncsys,swaplen={};swaprate={};nsw0=1;swaps=Sort[swap];spl=Split[Transpose[swaps][[1]]];Do[nswi=Length[spl[[i]]];swlst=Take[swaps,{nsw0,nswi+nsw0-1}];swpr=N[(Plus@@(Transpose[swlst][[2]]))/nswi];swaprate=Append[swaprate,swpr];swaplen=Append[swaplen,nswi];nsw0=nsw0+nswi,{i,1,nβ-1}];swaprate,]

The following is a matrix of Gaussian parameter proposal σ's chossen by the control system. There should be nβ lines, one for each tempering level. The top row corresponds to β0 and the bottom to β = 1. The columns give the proposal σ's for the parameters.

minσ

MatrixForm[σall]

A number of specialized diagnostic, only of use to designer.

The line below lists the most probable parameter set and some diagnotics recorded at each MCMC iteration.

maxlogpdfAll

Earlier terminatons of the control system.

ncsysold

A control system diagnostic.

swminσ

The acceptance or non cceptance of a swap between two adjacent chains is appended to the variable "swap" which has the format {k1+k2,1,i} for acceptance or  {k1+k2,0,i} for non acceptance. The next line of code sorts the swaps according to k1+k2 value, then splits the k1+k2 component of the sorted swaps, determines how many (nswi) swaps were between the upper to chains.

If[nβ>1,swaps=Sort[swap];spl=Split[Transpose[swaps][[1]]];nswi=Length[spl[[nβ-1]]];swlst=Take[swaps,-nswi]];

Look at the first 20 swaps between top two chains.

If[nβ>1,Take[swlst,20]]

Look at the last 20 swaps between top two chains.

If[nβ>1,Take[swlst,-20]]

MatrixForm[σold1]

MatrixForm[σold2]

MatrixForm[N[accept1/n1]]

Times@@(accept1[[nβ]]/N[n1])

MatrixForm[N[accept2/n1]]

Times@@(accept2[[nβ]]/N[n1])

MatrixForm[N[accept3/n1]]

cs2cnt2

Plotting results of MCMC

21. Setup some graph axes labels

Klabel=
Style[StringForm["``  (m (s^(-1)))",
Style["K",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

Vlabel=
Style[StringForm["``  (m (s^(-1)))",
Style["V",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

elabel=
Style[StringForm["Eccentricity  ``",
Style["e",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

ωlabel=
Style[StringForm["Longitude of periastron  ``",
Style["ω"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

ω1label=
Style[StringForm["``",
Style["(Subscript[ω, 1])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

ω2label=
Style[StringForm["``",
Style["(Subscript[ω, 2])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

ω3label=
Style[StringForm["``",
Style["(Subscript[ω, 3])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

ω4label=
Style[StringForm["``",
Style["(Subscript[ω, 4])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

ω5label=
Style[StringForm["``",
Style["(Subscript[ω, 5])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

blabel=
Style[StringForm["Noise scale parameter  ``",
Style["b",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

χlabel=
Style[StringForm["``",
Style["χ"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

χ1label=
Style[StringForm["``",
Style["(Subscript[χ, 1])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

χ2label=
Style[StringForm["``",
Style["(Subscript[χ, 2])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

χ3label=
Style[StringForm["``",
Style["(Subscript[χ, 3])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

χ4label=
Style[StringForm["``",
Style["(Subscript[χ, 4])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

χ5label=
Style[StringForm["``",
Style["(Subscript[χ, 5])"]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

asinilabel=
Style[StringForm["``",
Style["a  (AU)",Italic],Style["i",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

asinilabel1=
Style[StringForm["``",
Style["\!\(\*SubscriptBox[\(a\), \(\(1\)\(\\ \\ \\ \)\)]\)(AU)",Italic],Style["i",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

asinilabel2=
Style[StringForm["``",
Style["\!\(\*SubscriptBox[\(a\), \(\(2\)\(\\ \\ \\ \)\)]\)(AU)",Italic],Style["i",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

asinilabel3=
Style[StringForm["``",
Style["\!\(\*SubscriptBox[\(a\), \(\(3\)\(\\ \\ \\ \)\)]\)(AU)",Italic],Style["i",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

asinilabel4=
Style[StringForm["``",
Style["\!\(\*SubscriptBox[\(a\), \(\(4\)\(\\ \\ \)\)]\)(AU)",Italic],Style["i",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

asinilabel5=
Style[StringForm["``",
Style["\!\(\*SubscriptBox[\(a\), \(\(5\)\(\\ \\ \\ \)\)]\)(AU)",Italic],Style["i",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

Msinilabel=
Style[StringForm["`` sin ``",
Style["M",Italic],Style["i    ((Subscript[M, J]))",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

Msinilabel1=
Style[StringForm["`` sin ``",
Style["(Subscript[M, 1])",Italic],Style["i    ((Subscript[M, J]))",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

Msinilabel2=
Style[StringForm["`` sin ``",
Style["(Subscript[M, 2])",Italic],Style["i    ((Subscript[M, J]))",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

Msinilabel3=
Style[StringForm["`` sin ``",
Style["(Subscript[M, 3])",Italic],Style["i    ((Subscript[M, J]))",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

Msinilabel4=
Style[StringForm["`` sin ``",
Style["(Subscript[M, 4])",Italic],Style["i    ((Subscript[M, J]))",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

Msinilabel5=
Style[StringForm["`` sin ``",
Style["(Subscript[M, 5])",Italic],Style["i    ((Subscript[M, J]))",Italic]],
FontFamily->Times,FontWeight->Bold,FontSize->13]

22. Extracting model parameter samples for the β = 1 MCMC chain

MatrixForm[logpdfβi[[1000]]]

All the results are stored by the variable logpdfβi.
First let's look at an example of the output for iteration 1000 given by logpdfβi[[1000]] for a three planet search.
{{1.`,-0.05137493019965652`,0.7726223866082063`,0.33473268083098506`,0.08497398310245971`,0.8648919062576799`,0,-2,10000,67.16844349549996`,54.926044164384045`}}
Here is MatrixForm[logpdfβi[[1000]]];
(({{1.`, -0.05137493019965652`, 0.7726223866082063`, 0.33473268083098506`, 0.08497398310245971`, 0.8648919062576799`, 0, -2, 10000, 67.16844349549996`, 54.926044164384045`}}))
The first row contains the result for β = 1.0, because in this run we only used one β value.
We want to extract the β = 1.0 values. Do this with tr = Transpose[logpdfβi][[nβ]];
The first item in each row is the value of β, 
The next 5 items are the current values of the parameters {"α", "β", "σ", "μ", "τ"}.
The next item is the order of the periods (not relevant in this problem so just reads 0).
The next number indicates  (# == 0 accepted I proposals, # == 1 accepted C proposals, # == -2 rejected I proposals, # == -1 rejected C proposals).
The next item is the iteration number.
The penultimate item is the Log of the likelihood, that is used in computing the 
     global likelihood of the model for the Bayes factor calculation.
The last item is the value of   prior*(((like))^β) for the particular set of parameter values.
To obtain columns of each quantity use trn = Transpose[tr];
Finally we obtain the complete list of each parameter iterations with  a1good=trn[[2]]; a2good=trn[[3]]; etc.

Extract parameter iterations for the k = nβ member of the tempering levels.

k=nβ;tr=Transpose[logpdfβi][[k]];trn=Transpose[tr];ng=Length[trn[[1]]]

23. Examine parameter iterations for the β = 1 (k = nβ) MCMC chain

Start by choosing some display parameters.

niter=30;(* with niter = 100, plot only every 100th point *)

nburn1=10000;(* skip first nburn1 iterations for 2nd (upper right) plot *)

(*pl=Table[0,{i,np+2+nplan}];*)(* initiate table of graphs to be plotted *)

pl=Table[0,{i,np+2}];(* initiate table of graphs to be plotted *)

fn=Table[i,{i,1,ng}];(* initiate a variable representing the iteration number *)

Here is where we create the plots for the np parameters

xiterticks={25000,50000,75000,100000};

(*xiterticks={10000,20000};*)

Do[ label=StringForm[name[[i]]];If[i==99,pl[[i+2]]=ListLogPlot[Take[Transpose[{fn,trn[[i+1]]}],{1,ng,niter}],PlotRange->All,Joined->False,Frame->True,Axes->False,PlotStyle->{Black,PointSize[0.008]},FrameTicks->{xiterticks,Automatic,Automatic,Automatic},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iteration",{"Helvetica-Bold",12}],Style[label,{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}],pl[[i+2]]=ListPlot[Take[Transpose[{fn,trn[[i+1]]}],{1,ng,niter}],PlotRange->All,Joined->False,Frame->True,Axes->False,PlotStyle->{Black,PointSize[0.008]},FrameTicks->{xiterticks,Automatic,Automatic,Automatic},AspectRatio->0.5,
PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iteration",{"Helvetica-Bold",12}],Style[label,{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
],{i,1,np}
];

Next create a plot of (Subscript[Log, 10]) (Prior × Likelihood) vs iteration.

conv=Log[10,Exp[1.0]];pl[[1]]=ListPlot[Take[Transpose[{fn,conv*trn[[nl]]}],{1,ng,niter}],PlotRange->All,Joined->False,Frame->True,Axes->False,PlotStyle->{Black,PointSize[0.008]},FrameTicks->{xiterticks,Automatic,Automatic,Automatic},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iteration  chain  (β = "<>ToString[βlis[[k]]]<>")",{"Helvetica-Bold",12}],Style["(Subscript[Log, 10]) (Prior × Like)",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}];

Now create a plot of (Subscript[Log, 10]) (Prior × Likelihood) dropping first nburn1 iterations.

pl[[2]]=ListPlot[Take[Take[Transpose[{fn,conv*trn[[nl]]}],-(ng-nburn1)],{1,ng-nburn1,niter}],PlotRange->All,Joined->False,Frame->True,PlotStyle->{Black,PointSize[0.008]},FrameTicks->{xiterticks,Automatic,Automatic,Automatic},Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iteration  chain  (β = "<>ToString[βlis[[k]]]<>")",{"Helvetica-Bold",12}],Style["(Subscript[Log, 10]) (Prior × Like)",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}];

Next line shows the final plots. Note: for a nearly circular orbit the two parameters χ and ω are highly correlated
and are separately not well determined. Even then, the combination 2 π χ + ω is well constrained.

plIter=Show[GraphicsGrid[If[EvenQ[np+2],Table[{pl[[ii]],pl[[ii+1]]},{ii,1,np+2,2}],Join[Table[{pl[[ii]],pl[[ii+1]]},{ii,1,IntegerPart[(np+2)/2]*2,2}],{{pl[[IntegerPart[(np+2)/2]*2+1]]}}]],Spacings->{Scaled[-0.15],Scaled[-0.05]}],DisplayFunction->$DisplayFunction]

plIter2=Show[GraphicsGrid[{{pl[[1]],pl[[3]]},{pl[[4]],pl[[5]]},{pl[[6]],pl[[7]]}},Spacings->{Scaled[-0.15],Scaled[0.05]}]]

(*Export["c:/Users/Philip/Documents/greg13/MCMCregression/fig/CorYesStLineFitIntrinsicScatYesM7informGausXtPriorAnalyticX5fast26Aug13traces.eps", plIter2,ImageSize->72*9]*)

NotebookSave[]

Examine parameter iterations for the  k = nβ  MCMC chain for all iterations.

Plot (Subscript[Log, 10]) (Prior×Likelihood) vs intrinsic scatter parameter σ
for different chains specified by the k value

k=nβ;cutoff=-10;niter=10;trn=Transpose[Take[Transpose[logpdfβi][[k]],{1,nnt}]];v1good=trn[[4]];iter=trn[[n2l-1]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,v1good,iter}];psel=Sort[ppf,OrderedQ[{#2, #1}]&];ppsel=Select[psel,#[[1]]>cutoff&];nsel=Length[ppsel];pppsel=Take[ppsel,{1,nsel,niter}];tr=Transpose[pppsel];(* above line orders the items in descending period *)
lisv1=Transpose[{tr[[2]],tr[[1]]}];
xiterticks={0.0,0.25,0.5,1.0,1.5};
plv1Col=ListPlot[lisv1,Frame->True,PlotRange->{{0,1.5},All},PlotStyle->{{Black},{Red},{Blue},{Green},{Magenta},{Orange}},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameTicks->{xiterticks,Automatic,Automatic,Automatic},FrameLabel->{Style["σ",{"Helvetica-Bold",16}],Style["(Subscript[Log, 10]) (Prior\[Cross]Likelihood)",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Plot (Subscript[Log, 10]) (Prior×Likelihood) vs  slope parameter
for chain specified by the k value

k=nβ;cutoff=-10;niter=10;trn=Transpose[Take[Transpose[logpdfβi][[k]],{1,nnt}]];v1good=trn[[3]];iter=trn[[n2l-1]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,v1good,iter}];psel=Sort[ppf,OrderedQ[{#2, #1}]&];ppsel=Select[psel,#[[1]]>cutoff&];nsel=Length[ppsel];pppsel=Take[ppsel,{1,nsel,niter}];tr=Transpose[pppsel];(* above line orders the items in descending period *)
lisv1=Transpose[{tr[[2]],tr[[1]]}];
xiterticks={0.0,0.25,0.5,0.75,1.0,1.25,1.5,2};
plv1Col=ListPlot[lisv1,Frame->True,PlotRange->{{0,2},All},PlotStyle->{{Black},{Red},{Blue},{Green},{Magenta},{Orange}},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameTicks->{xiterticks,Automatic,Automatic,Automatic},FrameLabel->{Style["β",{"Helvetica-Bold",16}],Style["(Subscript[Log, 10]) (Prior\[Cross]Likelihood)",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Plot intrinsic scatter parameter σ  vs  slope parameter
for chain specified by the k value

k=nβ;cutoff=-10;niter=30;trn=Transpose[Take[Transpose[logpdfβi][[k]],{1,nnt}]];v1good=trn[[3]];v2good=trn[[4]];iter=trn[[n2l-1]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,v1good,v2good,iter}];psel=Sort[ppf,OrderedQ[{#2, #1}]&];ppsel=Select[psel,#[[1]]>cutoff&];nsel=Length[ppsel];pppsel=Take[ppsel,{1,nsel,niter}];tr=Transpose[pppsel];(* above line orders the items in descending period *)
lisP1e=Transpose[{tr[[2]],tr[[3]]}];
(*xiterticks={1,7.2,28,106,186,374,1000,3000,10000,30000,100000,300000};*)
plEccPCol=ListPlot[{lisP1e},Frame->True,PlotRange->{{0,2},{0,1.5}},PlotStyle->{{Black},{Red},{Blue},{Green},{Magenta},{Orange}},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],(*FrameTicks->{xiterticks,Automatic,Automatic,Automatic},*)FrameLabel->{Style["β",{"Helvetica-Bold",14}],Style["σ",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

NotebookSave[]

24. Set the burn-in value nburn.

The burn - in period is the greater of : (a) no. of iterations at which all the parameter distributions have converged, and (b) the point at which the control system is turned off.

ncsys is the iteration number at which the control system achieved its goal of automatically 
finding a good set of parameter proposal distribution σ's, and was turned off. 
Note: since only 1 out of every nthin iterations is kept the effective iteration at which the 
         control system cuts off is Ceiling[N[ncsys]/nthin].

ncsys

minimumBurnin=Ceiling[ncsys/nthin]

(* If facc1 >0, correlated proposal sampler employed and nburn = itsub, otherwise nburn = ncsys/nthin *)

itsub

nburn =If[facc1>0.0,itsub,minimumBurnin](* Default choice of burn-in = minimumBurnin *)

Examine the parameter iterations (traces) above to see it is safe to set nburn,
the iteration number at which the proposal σ control system is turned off.

ncsys=5000;

nburn=2566;

More detailed examination of 'I' and 'C' proposal acceptance rates.

Calculation of acceptance rates for  β = 1 chain.
 (# == 0 accepted I proposals, # == 1 accepted C proposals, # == -2 rejected I proposals, # == -1 rejected C proposals)

ncor1=Length[Select[Take[Transpose[Transpose[logpdfβi][[nβ]]][[nl-3]],{nburn,nnt}],#==-1&]]

ncor2=Length[Select[Take[Transpose[Transpose[logpdfβi][[nβ]]][[nl-3]],{nburn,nnt}],#==1&]]

Acceptance fraction for C proposals = acceptC

acceptC=N[ncor2/(ncor1+ncor2)]

ncor3=Length[Select[Take[Transpose[Transpose[logpdfβi][[nβ]]][[nl-3]],{nburn,nnt}],#==-2&]]

ncor4=Length[Select[Take[Transpose[Transpose[logpdfβi][[nβ]]][[nl-3]],{nburn,nnt}],#==0&]]

Acceptance fraction for I proposals = acceptI

acceptI=N[ncor4/(ncor4+ncor3)]

Overall acceptance fraction for I & C proposals = acceptIC

acceptIC=N[(ncor2+ncor4)/(ncor1+ncor2+ncor3+ncor4)]

Plot showing the distribution of the four categories of acceptance and rejection

ListPlot[Take[Transpose[Transpose[logpdfβi][[nβ]]][[nl-3]],{nburn,nnt,50}],Frame->True]

Examine parameter iterations for the  k = nβ  MCMC chain post burn-in.

(*k=nβ;tr=dataMC;trn=Transpose[tr];ng=Length[trn[[1]]]*)

(*k=nβ;tr=Take[Transpose[logpdfβi][[k]],76052];trn=Transpose[tr];ng=Length[trn[[1]]]*)

Length[tr]

k=nβ;tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt}];trn=Transpose[tr];ng=Length[trn[[1]]]

Start by choosing some display parameters.

niter=100;(* with niter = 100, plot only every 100th point *)

nburn1=10000;(* skip first nburn1 iterations for 2nd (upper right) plot *)

(*pl=Table[0,{i,np+2+nplan}];*)(* initiate table of graphs to be plotted *)

pl=Table[0,{i,np+2+2*nplan}];(* initiate table of graphs to be plotted *)

fn=Table[i,{i,1,ng}];(* initiate a variable representing the iteration number *)

Here is where we create the plots for the np parameters

(*xiterticks={25000,50000,75000};*)

xiterticks={10000,20000};

Do[ label=StringForm[name[[i]]];
If[i==99,pl[[i+2]]=ListLogPlot[Take[Transpose[{fn,trn[[i+1]]}],{1,ng,niter}],PlotRange->All,Joined->False,Frame->True,Axes->False,PlotStyle->{Black,PointSize[0.008]},FrameTicks->{xiterticks,Automatic,Automatic,Automatic},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iteration",{"Helvetica-Bold",12}],Style[label,{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}],pl[[i+2]]=ListPlot[Take[Transpose[{fn,trn[[i+1]]}],{1,ng,niter}],PlotRange->All,Joined->False,Frame->True,Axes->False,PlotStyle->{Black,PointSize[0.008]},FrameTicks->{xiterticks,Automatic,Automatic,Automatic},AspectRatio->0.5,
PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iteration",{"Helvetica-Bold",12}],Style[label,{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
],{i,1,np}
];

Next create a plot of (Subscript[Log, 10]) (Prior × Likelihood) vs iteration.

conv=Log[10,Exp[1.0]];pl[[1]]=ListPlot[Take[Transpose[{fn,conv*trn[[nl]]}],{1,ng,niter}],PlotRange->All,Joined->False,Frame->True,Axes->False,PlotStyle->{Black,PointSize[0.008]},FrameTicks->{xiterticks,Automatic,Automatic,Automatic},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iteration  chain  (β = "<>ToString[βlis[[k]]]<>")",{"Helvetica-Bold",12}],Style["(Subscript[Log, 10]) (Prior × Like)",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}];

Now create a plot of (Subscript[Log, 10]) (Prior × Likelihood) dropping first nburn1 iterations.

pl[[2]]=ListPlot[Take[Take[Transpose[{fn,conv*trn[[nl]]}],-(ng-nburn1)],{1,ng-nburn1,niter}],PlotRange->All,Joined->False,Frame->True,PlotStyle->{Black,PointSize[0.008]},FrameTicks->{xiterticks,Automatic,Automatic,Automatic},Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iteration  chain  (β = "<>ToString[βlis[[k]]]<>")",{"Helvetica-Bold",12}],Style["(Subscript[Log, 10]) (Prior × Like)",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}];

Next line shows the final plots. Note: for a nearly circular orbit the two parameters χ and ω are highly correlated
and are separately not well determined. Even then, the combination 2 π χ + ω is well constrained.

plIter=Show[GraphicsGrid[If[EvenQ[np+2],Table[{pl[[ii]],pl[[ii+1]]},{ii,1,np+2,2}],Join[Table[{pl[[ii]],pl[[ii+1]]},{ii,1,IntegerPart[(np+2)/2]*2,2}],{{pl[[IntegerPart[(np+2)/2]*2+1]]}}]],Spacings->{Scaled[-0.15],Scaled[-0.05]}],DisplayFunction->$DisplayFunction]

Plot (Subscript[Log, 10]) (Prior×Likelihood) vs intrinsic scatter parameter σ
for chain specified by the k value

Specify the interval in iteration number to plot.

cutoff=-10;

interval=30;xiterticks={};

k=nβ;cutoff=-10;niter=40;trn=Transpose[Take[Transpose[logpdfβi][[k]],{nburn,nnt}]];v1good=trn[[4]];iter=trn[[n2l-1]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,v1good,iter}];psel=Sort[ppf,OrderedQ[{#2, #1}]&];ppsel=Select[psel,#[[1]]>cutoff&];nsel=Length[ppsel];pppsel=Take[ppsel,{1,nsel,niter}];tr=Transpose[pppsel];(* above line orders the items in descending period *)
lisv1=Transpose[{tr[[2]],tr[[1]]}];
xiterticks={0.0,0.25,0.5,1.0,1.5};
plv1=ListPlot[lisv1,Frame->True,PlotRange->{{0,1.5},All},PlotStyle->{{Black},{Red},{Blue},{Green},{Magenta},{Orange}},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameTicks->{xiterticks,Automatic,Automatic,Automatic},FrameLabel->{Style["σ",{"Helvetica-Bold",16}],Style["(Subscript[Log, 10]) (Prior\[Cross]Likelihood)",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

(*Export["c:/Users/Philip/Documents/greg13/MCMCregression/fig/CorYesStLineFitIntrinsicScatYesM7informGausXtPriorAnalyticX5fast26Aug13LikeVsSigma.eps", plv1,ImageSize->72*7]*)

plv1L=Show[plv1,Graphics[{Text["(a)",{1.3,-6.5}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5]

Plot (Subscript[Log, 10]) (Prior×Likelihood) vs  slope parameter
for chain specified by the k value

k=nβ;cutoff=-10;niter=10;trn=Transpose[Take[Transpose[logpdfβi][[k]],{nburn,nnt}]];v1good=trn[[3]];iter=trn[[n2l-1]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,v1good,iter}];psel=Sort[ppf,OrderedQ[{#2, #1}]&];ppsel=Select[psel,#[[1]]>cutoff&];nsel=Length[ppsel];pppsel=Take[ppsel,{1,nsel,niter}];tr=Transpose[pppsel];(* above line orders the items in descending period *)
lisv1=Transpose[{tr[[2]],tr[[1]]}];
xiterticks={0.0,0.25,0.5,0.75,1.0,1.25,1.5,2};
plv1Col=ListPlot[lisv1,Frame->True,PlotRange->{{0,2},All},PlotStyle->{{Black},{Red},{Blue},{Green},{Magenta},{Orange}},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameTicks->{xiterticks,Automatic,Automatic,Automatic},FrameLabel->{Style["β",{"Helvetica-Bold",16}],Style["(Subscript[Log, 10]) (Prior\[Cross]Likelihood)",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Plot intrinsic scatter parameter σ  vs  slope parameter
for chain specified by the k value

k=nβ;cutoff=-10;niter=30;trn=Transpose[Take[Transpose[logpdfβi][[k]],{nburn,nnt}]];v1good=trn[[3]];v2good=trn[[4]];iter=trn[[n2l-1]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,v1good,v2good,iter}];psel=Sort[ppf,OrderedQ[{#2, #1}]&];ppsel=Select[psel,#[[1]]>cutoff&];nsel=Length[ppsel];pppsel=Take[ppsel,{1,nsel,niter}];tr=Transpose[pppsel];(* above line orders the items in descending period *)
lisP1e=Transpose[{tr[[2]],tr[[3]]}];
(*xiterticks={1,7.2,28,106,186,374,1000,3000,10000,30000,100000,300000};*)
plEccPCol=ListPlot[{lisP1e},Frame->True,PlotRange->{{0,2},{0,1.5}},PlotStyle->{{Black},{Red},{Blue},{Green},{Magenta},{Orange}},AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],(*FrameTicks->{xiterticks,Automatic,Automatic,Automatic},*)FrameLabel->{Style["β",{"Helvetica-Bold",14}],Style["σ",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Plot slope vs intercept 
for chain specified by the k value

Specify the interval in iteration number to plot.

cutoff=-10;

interval=30;xiterticks={};

Timing[plPerProb=Table[0,{nβ}];k=nβ;conv=Log[10,Exp[1.0]];
tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt}];(* extract values for the k=nβ chain *)
str=RotateLeft[Sort[RotateRight[tr,{0,1}],OrderedQ[{#2, #1}]&],{0,1}];(* rotate so prob is first item, sort in descending order, rotate back *)
besta1=str[[1,2]];
trn=Transpose[str];a1good=trn[[2]];a2good=trn[[3]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,a1good,a2good}];pp=Select[ppf,#[[1]]>cutoff&];nppf=Length[pp];pint=Take[pp,{1,nppf,interval}];tpint=Transpose[pint];ppsel=Transpose[{tpint[[2]],tpint[[3]]}];nsel=Length[ppsel];pselA=Sort[ppsel,OrderedQ[{#2, #1}]&];(* above line orders the items in descending period *)
pntcol=ListPlot[pselA,PlotRange->All];
plPerProb[[k]]=Show[pntcol,
Frame->True,AspectRatio->0.5,Axes->False,
FrameLabel->{Style["α",{"Helvetica-Bold",12}],Style["β",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
]

Plot the slope versus intercept correlation.

(*xiterticks={1,3,5.37,30,100,300,1000,3000,10000,30000};*)plEccPCol=ListPlot[{pselA},Frame->True,Axes->False,PlotStyle->{{Black},{Green},{Red},{Magenta},{Orange},{Black}},PlotRange->All,AspectRatio->0.5,FrameLabel->{Style["α",{"Helvetica-Bold",16}],Style[" β",{"Helvetica-Bold",16}]},PlotStyle->{{Blue},{Orange},{Red},{Green},{Black},{Magenta}},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Plot intrinsic scatter μ vs slope 
for different chains specified by the k value

Specify the interval in iteration number to plot.

cutoff=-10;

interval=50;xiterticks={};

One plot for each of the nβ tempering levels.

Length[pselA]

Timing[plPerProb=Table[0,{nβ}];k=nβ;conv=Log[10,Exp[1.0]];
tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt}];(* extract values for the k=nβ chain *)
str=RotateLeft[Sort[RotateRight[tr,{0,1}],OrderedQ[{#2, #1}]&],{0,1}];(* rotate so prob is first item, sort in descending order, rotate back *)
besta1=str[[1,2]];
trn=Transpose[str];a1good=trn[[3]];a2good=trn[[5]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,a1good,a2good}];pp=Select[ppf,#[[1]]>cutoff&];nppf=Length[pp];pint=Take[pp,{1,nppf,interval}];tpint=Transpose[pint];ppsel=Transpose[{tpint[[2]],tpint[[3]]}];nsel=Length[ppsel];pselA=Sort[ppsel,OrderedQ[{#2, #1}]&];(* above line orders the items in descending period *)
pntcol=ListPlot[pselA,PlotRange->All];
plPerProb[[k]]=Show[pntcol,
Frame->True,AspectRatio->0.5,Axes->False,
FrameLabel->{Style["β",{"Helvetica-Bold",12}],Style["μ",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
]

(*FrameTicks->{xiterticks,Automatic,Automatic,Automatic}*)

Plot the μ versus slope correlation.

(*xiterticks={1,3,5.37,30,100,300,1000,3000,10000,30000};*)plEccPCol=ListPlot[{pselA},Frame->True,Axes->False,PlotStyle->{{Black},{Green},{Red},{Magenta},{Orange},{Black}},PlotRange->All,AspectRatio->0.5,FrameLabel->{Style["β",{"Helvetica-Bold",16}],Style["μ",{"Helvetica-Bold",16}]},PlotStyle->{{Blue},{Orange},{Red},{Green},{Black},{Magenta}},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Plot intrinsic scatter σ vs slope 
for different chains specified by the k value

Specify the interval in iteration number to plot.

cutoff=-2;

interval=30;xiterticks={};

One plot for each of the nβ tempering levels.

Length[pselA]

(*FrameTicks->{xiterticks,Automatic,Automatic,Automatic}*)

Timing[plPerProb=Table[0,{nβ}];k=nβ;conv=Log[10,Exp[1.0]];
tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt}];(* extract values for the k=nβ chain *)
str=RotateLeft[Sort[RotateRight[tr,{0,1}],OrderedQ[{#2, #1}]&],{0,1}];(* rotate so prob is first item, sort in descending order, rotate back *)
besta1=str[[1,2]];
trn=Transpose[str];a1good=trn[[3]];a2good=trn[[4]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,a1good,a2good}];pp=Select[ppf,#[[1]]>cutoff&];nppf=Length[pp];pint=Take[pp,{1,nppf,interval}];tpint=Transpose[pint];ppsel=Transpose[{tpint[[2]],tpint[[3]]}];nsel=Length[ppsel];pselA=Sort[ppsel,OrderedQ[{#2, #1}]&];(* above line orders the items in descending period *)
pntcol=ListPlot[pselA,PlotRange->All];
plPerProb[[k]]=Show[pntcol,
Frame->True,AspectRatio->0.5,Axes->False,
FrameLabel->{Style["β",{"Helvetica-Bold",12}],Style["σ",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
]

Plot  σ vs slope correlation.

(*xiterticks={1,3,5.37,30,100,300,1000,3000,10000,30000};*)plEccPCol=ListPlot[{pselA},Frame->True,Axes->False,PlotStyle->{{Black},{Green},{Red},{Magenta},{Orange},{Black}},PlotRange->All,AspectRatio->0.5,FrameLabel->{Style["β",{"Helvetica-Bold",12}],Style["σ",{"Helvetica-Bold",12}]},PlotStyle->{{Blue},{Orange},{Red},{Green},{Black},{Magenta}},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Plot intrinsic scatter μ vs τ 
for different chains specified by the k value

Specify the interval in iteration number to plot.

cutoff=-10;

interval=30;xiterticks={};

One plot for each of the nβ tempering levels.

Length[pselA]

(*FrameTicks->{xiterticks,Automatic,Automatic,Automatic}*)

Timing[plPerProb=Table[0,{nβ}];k=nβ;conv=Log[10,Exp[1.0]];
tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt}];(* extract values for the k=nβ chain *)
str=RotateLeft[Sort[RotateRight[tr,{0,1}],OrderedQ[{#2, #1}]&],{0,1}];(* rotate so prob is first item, sort in descending order, rotate back *)
besta1=str[[1,2]];
trn=Transpose[str];a1good=trn[[5]];a2good=trn[[6]];prob=trn[[nl]];ng=Length[trn[[1]]];probMax=Max[prob];logprob = conv*(prob - probMax);ppf=Transpose[{logprob,a1good,a2good}];pp=Select[ppf,#[[1]]>cutoff&];nppf=Length[pp];pint=Take[pp,{1,nppf,interval}];tpint=Transpose[pint];ppsel=Transpose[{tpint[[2]],tpint[[3]]}];nsel=Length[ppsel];pselA=Sort[ppsel,OrderedQ[{#2, #1}]&];(* above line orders the items in descending period *)
pntcol=ListPlot[pselA,PlotRange->All];
plPerProb[[k]]=Show[pntcol,
Frame->True,AspectRatio->0.5,Axes->False,
FrameLabel->{Style["μ",{"Helvetica-Bold",12}],Style["τ",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
]

Plot μ vs τ  correlation.

(*xiterticks={1,3,5.37,30,100,300,1000,3000,10000,30000};*)plEccPCol=ListPlot[{pselA},Frame->True,Axes->False,PlotStyle->{{Black},{Green},{Red},{Magenta},{Orange},{Black}},PlotRange->All,AspectRatio->0.5,FrameLabel->{Style["μ",{"Helvetica-Bold",16}],Style["τ",{"Helvetica-Bold",16}]},PlotStyle->{{Blue},{Orange},{Red},{Green},{Black},{Magenta}},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

Plot (Subscript[Log, 10])[Prior × Likelihood] versus MCMC iteration.

k=nβ;niter=100;tr=Transpose[Take[Transpose[logpdfβi][[k]],{1,nnt,niter}]];
lisLogProb=Join[{{-0.03,Min[tr[[nl]]]}},Transpose[{tr[[n2l-1]]*(10^(-6)),tr[[nl]]}]];plProb=ListPlot[lisLogProb,Frame->True,PlotRange->All,Axes->False,AspectRatio->0.5,PlotStyle->Black,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Iterations  (\[Cross] (10^6))",{"Helvetica-Bold",14}],Style["(Subscript[Log, 10])[Prior \[Cross] Like]",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

NotebookSave[]

Plot marginal  distribution for model parameters

Extract parameter iterations beyond burn-in for plotting. Specify the tempering level, normally k = nβ.

k=nβ;tr=Take[Transpose[logpdfβi][[k]],-(nnt-nburn)];trn=Transpose[tr];trnn=Drop[trn,{1}];ng=Length[trnn[[1]]]

Plot histograms of all the parameter iterations.

plh=Table[1,{i,np}];Do[label=StringForm[name[[i]]<>unit[[i]]];
plh[[i]]=Histogram[trnn[[i]],PlotRange->{{Min[trnn[[i]]],Max[trnn[[i]]]},All},Frame->True,ChartStyle->EdgeForm[],Axes->False,AspectRatio->0.5,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style[label,{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}],{i,1,np}
];

plhis=Show[GraphicsGrid[If[EvenQ[np],Table[{plh[[ii]],plh[[ii+1]]},{ii,1,np,2}],Join[Table[{plh[[ii]],plh[[ii+1]]},{ii,1,IntegerPart[(np)/2]*2,2}],{{plh[[IntegerPart[(np)/2]*2+1]]}}]],Spacings->{Scaled[-0.2],Scaled[-0.02]}]]

Plot marginal  distribution for model parameters for most common peak.

Extract parameter iterations beyond burn-in with period values between pLow and pHi.
By default pLow = avmin[[1]] and pHi = avMax[[1]], the search range.

trn=.;k=nβ;
If[ncsys<n,a1Low=-2.0;a1Hi=2.0;tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt,1}];trPsel=Select[tr,#[[2]]>a1Low&&#[[2]]<a1Hi&];lsel=Length[trPsel];If[lsel>0,trn=Transpose[trPsel];trnn=Drop[trn,{1}]];ng=Length[trnn[[1]]]]

Details behind generating the marginal distributions (using nai bins).

Specify the number of points nai to bin the continuous parameters and set up some tables we will need for some parameter distribution statistics.

nai=150;plp=Table[1,{i,np}];meanai=Table[0,{j,1,np}];medianai=Table[0,{j,1,np}];modeai=Table[0,{j,1,np}];minai=Table[0,{j,1,np}];maxai=Table[0,{j,1,np}];areai=Table[0,{j,1,np}];

meanai

Next block computes the marginal PDF for the continuous parameters.
The text in red is devoted to handling the parameters χ and ω in a special way since
they are Modulo 1 and 2π, respectively. A problem can arise in calculating credible regions
if the PDF has a bump at both ends of the range and a gap in the middle. The two ends are 
actually close to one another in a wrap around sense because they are angular coordinates.
The aim of the red code is to put them together to achieve the minimum size credible region.

Recal we specified search ranges initially of
avmin[[4]]=0; avMax[[4]]=1.0;avmin[[6]]=0; avMax[[4]]=2π; 
The first step is to determine the repeat length of the angular parameters
If[avMax[[i]]>1.,useMax=2.0*π,useMax=1.0];
  Then add a second (repeat) cycle of each with
aimarg2=Join[aimarg1,(aimarg1+avMax[[i]]/2.)];
  Then we make a binned version (2 nai1 bins in total) of the double cycle.
  The idea is to Select one cycle (nai1 bins) of each varying the starting phase. 
  We compute the variance of each Select and note the start phase strtpbin2[[i]]
  of the first Select which has the minimum variance.
  Finally we Select one cycle from aimarg2 starting at a value strtpbin2[[i]].

strtpbin2=Table[0,{np+2}];swstrt=0;nai1=100;
If[ncsys<n&&lsel>0,
Do[label=StringForm[name[[i]]<>unit[[i]]];aidiff=0;aimarg=trn[[i+1]];swstrt=1;
If[aidiff>0,,meanai[[i]]=Mean[aimarg];medianai[[i]]=Median[aimarg];low=Round[10000000*Min[aimarg]]/10000000;hi=Round[10000000*Max[aimarg]]/10000000;dai=(hi-low)/nai;aibin=BinCounts[aimarg,{low,hi,dai}];aimidpoints=Table[N[k],{k,low+dai/2,hi-dai/2,dai}];
aibintot=Plus@@aibin;yai=(aibin/(aibintot*dai));yaisum=Plus@@yai;pbin=Sort[Transpose[{yai/yaisum,aimidpoints}],OrderedQ[{#2,#1}]&];modeai[[i]]=pbin[[1,2]];minai[[i]]=1000000000;maxai[[i]]=-1000000000;area=0;j=1;
While[area≤0.683,area=area+ pbin[[j,1]];
If[pbin[[j,2]]<minai[[i]],minai[[i]]=pbin[[j,2]],Null];
If[pbin[[j,2]]>maxai[[i]],maxai[[i]]=pbin[[j,2]],Null];
j=j+1];areai[[i]]=area;
aimargbin=Transpose[{aimidpoints,yai}];
xscale=1/N[10^Floor[Log[10,(Max[aimidpoints]-Min[aimidpoints])]]];If[xscale>0,xaxisticks=Round[2*xscale*Table[low+(hi-low)*(1+2*(j-1))/4.0,{j,1,2}]]/(2*xscale),xaxisticks=Round[xscale*Table[low+(hi-low)*(1+2*(j-1))/4.0,{j,1,2}]]/(xscale)];
plp[[i]]=ListPlot[aimargbin,PlotRange->{{low,hi},{0,1.1*Max[yai]}},Frame->True,FrameTicks->{xaxisticks,Automatic,Automatic,Automatic}, Axes->False,PlotStyle->{Black,Thickness[0.007],PointSize[0.05]},Joined->True,AspectRatio->0.5,FrameLabel->{Style[label,{"Helvetica-Bold",16}],
		Style["PDF",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
],{i,1,np}
];
];

The line below gives the actual area within the credible region for each parameter achieved for the particular choice of the number of bins, nai.

If[ncsys<n&&lsel>0&&aidiff<1,N[areai],]

The 4th and 6th items of strtpbin2 are the lower boundaries of the phase windows for χ and ω 
that lead to the minimum credible regions for these two items..

strtpbin2

meanai

medianai

modeai

minai

maxai

Next step does the plotting.

If[ncsys<n&&lsel>0&&aidiff<1,plmarg1all=Show[GraphicsGrid[If[EvenQ[np],Table[{plp[[ii]],plp[[ii+1]]},{ii,1,np,2}],Join[Table[{plp[[ii]],plp[[ii+1]]},{ii,1,IntegerPart[(np)/2]*2,2}],{{plp[[IntegerPart[(np)/2]*2+1]]}}]],Spacings->{Scaled[-0.1],Scaled[0.05]}]]]

Display marginal  distributions together with true parameter values.

trueIntercept

plalphaLine=Show[plp[[1]],Graphics[{Thickness[0.005],Line[{{trueIntercept,0},{trueIntercept,6}}]}]]

trueSlope

plbetaLine=Show[plp[[2]],Graphics[{Thickness[0.005],Line[{{trueSlope,0},{trueSlope,5}}]}]]

plsigmaLine=Show[plp[[3]],Graphics[{Thickness[0.005],Line[{{trueSigma,0},{trueSigma,5}}]}]]

plmuLine=Show[plp[[4]],Graphics[{Thickness[0.005],Line[{{trueMu,0},{trueMu,5}}]}]]

pltauLine=Show[plp[[5]],Graphics[{Thickness[0.005],Line[{{trueTau,0},{trueTau,5}}]}]]

plmarg3=Show[GraphicsGrid[{{plalphaLine,plbetaLine},{plsigmaLine}},Spacings->{Scaled[-0.1],Scaled[0.05]}]]

plmarg5=Show[GraphicsGrid[{{plalphaLine,plbetaLine},{plsigmaLine,plmuLine},{pltauLine}},Spacings->{Scaled[-0.15],Scaled[0.05]}]]

(*Export["c:/Users/Philip/Documents/greg13/MCMCregression/fig/regressLineAnalytic_Public_14May16M8_Marginals.eps", plmarg5,ImageSize->72*9]*)

Plot marginal  distribution for model parameters for most common peak. 
(moving average smoothing shorter)

Extract parameter iterations beyond burn-in with period values between pLow and pHi.
By default pLow = avmin[[1]] and pHi = avMax[[1]], the search range.

trn=.;k=nβ;
If[ncsys<n,a1Low=-2.0;a1Hi=2.0;tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt,1}];trPsel=Select[tr,#[[2]]>a1Low&&#[[2]]<a1Hi&];lsel=Length[trPsel];If[lsel>0,trn=Transpose[trPsel];trnn=Drop[trn,{1}]];ng=Length[trnn[[1]]]]

Details behind generating the marginal distributions (using nai bins).

Specify the number of points nai to bin the continuous parameters and set up some tables we will need for some parameter distribution statistics.

nai=150;plpS=Table[1,{i,np}];meanai=Table[0,{j,1,np}];medianai=Table[0,{j,1,np}];modeai=Table[0,{j,1,np}];minai=Table[0,{j,1,np}];maxai=Table[0,{j,1,np}];areai=Table[0,{j,1,np}];

meanai

Next block computes the marginal PDF for the continuous parameters (the (Subscript[P, 1]) is discrete).
The text in red is devoted to handling the parameters χ and ω in a special way since
they are Modulo 1 and 2π, respectively. A problem can arise in calculating credible regions
if the PDF has a bump at both ends of the range and a gap in the middle. The two ends are 
actually close to one another in a wrap around sense because they are angular coordinates.
The aim of the red code is to put them together to achieve the minimum size credible region.

Recal we specified search ranges initially of
avmin[[4]]=0; avMax[[4]]=1.0;avmin[[6]]=0; avMax[[4]]=2π; 
The first step is to determine the repeat length of the angular parameters
If[avMax[[i]]>1.,useMax=2.0*π,useMax=1.0];
  Then add a second (repeat) cycle of each with
aimarg2=Join[aimarg1,(aimarg1+avMax[[i]]/2.)];
  Then we make a binned version (2 nai1 bins in total) of the double cycle.
  The idea is to Select one cycle (nai1 bins) of each varying the starting phase. 
  We compute the variance of each Select and note the start phase strtpbin2[[i]]
  of the first Select which has the minimum variance.
  Finally we Select one cycle from aimarg2 starting at a value strtpbin2[[i]].

strtpbin2=Table[0,{np+2}];swstrt=0;nai1=100;
If[ncsys<n&&lsel>0,
Do[label=StringForm[name[[i]]<>unit[[i]]];aidiff=0;aimarg=trn[[i+1]];swstrt=1;
If[type[[i]]==2||type[[i]]==4,aimarg1=Mod[aimarg,2*π];useMax=2π;
       aimarg2=Join[aimarg1,(aimarg1+useMax)];low1=0.0;hi1=N[2*useMax];dai1=(hi1-low1)/(2*nai1);
aibin1=BinCounts[aimarg2,{low1,hi1,dai1}];aimidpoints1=Table[N[k],{k,low1+dai1/2,hi1-dai1/2,dai1}];ddiff=Length[aimidpoints1]-Length[aibin1];If[ddiff>0,aidiff=ddiff,];
If[aidiff>0,,pbin1=Transpose[{aibin1,aimidpoints1}];varbin0=10000000.0;Do[pbin2=Select[pbin1,#[[2]]>((j-1)*dai1)&&#[[2]]<(((j-1)+nai1)*dai1)&];tpbin2=Transpose[pbin2];sumpbin2=Plus@@tpbin2[[1]];meanpbin2=tpbin2[[1]].tpbin2[[2]]/(sumpbin2);varpbin2=tpbin2[[1]].(((tpbin2[[2]]-meanpbin2))^2)/(sumpbin2);If[varpbin2<varbin0,strtpbin2[[i]]=tpbin2[[2,1]]-dai1/2;varbin0=varpbin2,],{j,1,nai1}];
aimarg=Select[aimarg2,#>strtpbin2[[i]]&&#<(strtpbin2[[i]]+nai1*dai1)&]]];
If[aidiff>0,,meanai[[i]]=Mean[aimarg];medianai[[i]]=Median[aimarg];low=Round[10000000*Min[aimarg]]/10000000;hi=Round[10000000*Max[aimarg]]/10000000;dai=(hi-low)/nai;aibin=BinCounts[aimarg,{low,hi,dai}];aimidpoints=Table[N[k],{k,low+dai/2,hi-dai/2,dai}];
aibintot=Plus@@aibin;yai=(aibin/(aibintot*dai));yaisum=Plus@@yai;pbin=Sort[Transpose[{yai/yaisum,aimidpoints}],OrderedQ[{#2,#1}]&];modeai[[i]]=pbin[[1,2]];minai[[i]]=1000000000;maxai[[i]]=-1000000000;area=0;j=1;
While[area≤0.683,area=area+ pbin[[j,1]];
If[pbin[[j,2]]<minai[[i]],minai[[i]]=pbin[[j,2]],Null];
If[pbin[[j,2]]>maxai[[i]],maxai[[i]]=pbin[[j,2]],Null];
j=j+1];areai[[i]]=area;
aimargbin=Transpose[{aimidpoints,yai}];xscale=1/N[10^Floor[Log[10,(Max[aimidpoints]-Min[aimidpoints])]]];If[xscale>0,xaxisticks=Round[2*xscale*Table[low+(hi-low)*(1+2*(j-1))/4.0,{j,1,2}]]/(2*xscale),xaxisticks=Round[xscale*Table[low+(hi-low)*(1+2*(j-1))/4.0,{j,1,2}]]/(xscale)];plpS[[i]]=ListPlot[MovingAverage[aimargbin,{(1/24),(3/24),(5/24),(6/24),(5/24),(3/24),(1/24)}],PlotRange->{{low,hi},{0,1.1*Max[yai]}},Frame->True,FrameTicks->{xaxisticks,Automatic,Automatic,Automatic}, Axes->False,PlotStyle->{Black,Thickness[0.007],PointSize[0.05]},Joined->True,AspectRatio->0.5,FrameLabel->{Style[label,{"Helvetica-Bold",12}],
		Style["PDF",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
],{i,1,np}
];
];

The line below gives the actual area within the credible region for each parameter achieved for the particular choice of the number of bins, nai.

If[ncsys<n&&lsel>0&&aidiff<1,N[areai],]

The 4th and 6th items of strtpbin2 are the lower boundaries of the phase windows for χ and ω 
that lead to the minimum credible regions for these two items..

strtpbin2

meanai

medianai

modeai

minai

maxai

Next step does the plotting.

(*If[ncsys<n&&lsel>0&&aidiff<1,plmarg1a=Show[GraphicsGrid[{{plpS[[1]],plpS[[2]],plpS[[5]]},{plpS[[7]],,}},Spacings->{Scaled[-0.0],Scaled[0.05]}]]]*)

(*Export["c:/Users/Philip/Documents/greg8/47UMaNew/fig/APT3planGene3C47UMaNewOffsetsCorrHETHJSM7Marg.eps", plmarg1a,ImageSize->72*9]*)

If[ncsys<n&&lsel>0&&aidiff<1,plmarg1allS=Show[GraphicsGrid[If[EvenQ[np],Table[{plpS[[ii]],plpS[[ii+1]]},{ii,1,np,2}],Join[Table[{plpS[[ii]],plpS[[ii+1]]},{ii,1,IntegerPart[(np)/2]*2,2}],{{plpS[[IntegerPart[(np)/2]*2+1]]}}]],Spacings->{Scaled[-0.1],Scaled[0.05]}]]]

Plot marginal  distribution for model parameters for most common peak. 
(moving average smoothing shorter shorter)

Extract parameter iterations beyond burn-in with period values between pLow and pHi.
By default pLow = avmin[[1]] and pHi = avMax[[1]], the search range.

trn=.;k=nβ;
If[ncsys<n,a1Low=-2.0;a1Hi=2.0;tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt,1}];trPsel=Select[tr,#[[2]]>a1Low&&#[[2]]<a1Hi&];lsel=Length[trPsel];If[lsel>0,trn=Transpose[trPsel];trnn=Drop[trn,{1}]];ng=Length[trnn[[1]]]]

Details behind generating the marginal distributions (using nai bins).

Specify the number of points nai to bin the continuous parameters and set up some tables we will need for some parameter distribution statistics.

nai=150;plpS=Table[1,{i,np}];meanai=Table[0,{j,1,np}];medianai=Table[0,{j,1,np}];modeai=Table[0,{j,1,np}];minai=Table[0,{j,1,np}];maxai=Table[0,{j,1,np}];areai=Table[0,{j,1,np}];

meanai

Next block computes the marginal PDF for the continuous parameters (the (Subscript[P, 1]) is discrete).
The text in red is devoted to handling the parameters χ and ω in a special way since
they are Modulo 1 and 2π, respectively. A problem can arise in calculating credible regions
if the PDF has a bump at both ends of the range and a gap in the middle. The two ends are 
actually close to one another in a wrap around sense because they are angular coordinates.
The aim of the red code is to put them together to achieve the minimum size credible region.

Recal we specified search ranges initially of
avmin[[4]]=0; avMax[[4]]=1.0;avmin[[6]]=0; avMax[[4]]=2π; 
The first step is to determine the repeat length of the angular parameters
If[avMax[[i]]>1.,useMax=2.0*π,useMax=1.0];
  Then add a second (repeat) cycle of each with
aimarg2=Join[aimarg1,(aimarg1+avMax[[i]]/2.)];
  Then we make a binned version (2 nai1 bins in total) of the double cycle.
  The idea is to Select one cycle (nai1 bins) of each varying the starting phase. 
  We compute the variance of each Select and note the start phase strtpbin2[[i]]
  of the first Select which has the minimum variance.
  Finally we Select one cycle from aimarg2 starting at a value strtpbin2[[i]].

strtpbin2=Table[0,{np+2}];swstrt=0;nai1=100;
If[ncsys<n&&lsel>0,
Do[label=StringForm[name[[i]]<>unit[[i]]];aidiff=0;aimarg=trn[[i+1]];swstrt=1;
If[type[[i]]==2||type[[i]]==4,aimarg1=Mod[aimarg,2*π];useMax=2π;
       aimarg2=Join[aimarg1,(aimarg1+useMax)];low1=0.0;hi1=N[2*useMax];dai1=(hi1-low1)/(2*nai1);
aibin1=BinCounts[aimarg2,{low1,hi1,dai1}];aimidpoints1=Table[N[k],{k,low1+dai1/2,hi1-dai1/2,dai1}];ddiff=Length[aimidpoints1]-Length[aibin1];If[ddiff>0,aidiff=ddiff,];
If[aidiff>0,,pbin1=Transpose[{aibin1,aimidpoints1}];varbin0=10000000.0;Do[pbin2=Select[pbin1,#[[2]]>((j-1)*dai1)&&#[[2]]<(((j-1)+nai1)*dai1)&];tpbin2=Transpose[pbin2];sumpbin2=Plus@@tpbin2[[1]];meanpbin2=tpbin2[[1]].tpbin2[[2]]/(sumpbin2);varpbin2=tpbin2[[1]].(((tpbin2[[2]]-meanpbin2))^2)/(sumpbin2);If[varpbin2<varbin0,strtpbin2[[i]]=tpbin2[[2,1]]-dai1/2;varbin0=varpbin2,],{j,1,nai1}];
aimarg=Select[aimarg2,#>strtpbin2[[i]]&&#<(strtpbin2[[i]]+nai1*dai1)&]]];
If[aidiff>0,,meanai[[i]]=Mean[aimarg];medianai[[i]]=Median[aimarg];low=Round[10000000*Min[aimarg]]/10000000;hi=Round[10000000*Max[aimarg]]/10000000;dai=(hi-low)/nai;aibin=BinCounts[aimarg,{low,hi,dai}];aimidpoints=Table[N[k],{k,low+dai/2,hi-dai/2,dai}];
aibintot=Plus@@aibin;yai=(aibin/(aibintot*dai));yaisum=Plus@@yai;pbin=Sort[Transpose[{yai/yaisum,aimidpoints}],OrderedQ[{#2,#1}]&];modeai[[i]]=pbin[[1,2]];minai[[i]]=1000000000;maxai[[i]]=-1000000000;area=0;j=1;
While[area≤0.683,area=area+ pbin[[j,1]];
If[pbin[[j,2]]<minai[[i]],minai[[i]]=pbin[[j,2]],Null];
If[pbin[[j,2]]>maxai[[i]],maxai[[i]]=pbin[[j,2]],Null];
j=j+1];areai[[i]]=area;
aimargbin=Transpose[{aimidpoints,yai}];xscale=1/N[10^Floor[Log[10,(Max[aimidpoints]-Min[aimidpoints])]]];If[xscale>0,xaxisticks=Round[2*xscale*Table[low+(hi-low)*(1+2*(j-1))/4.0,{j,1,2}]]/(2*xscale),xaxisticks=Round[xscale*Table[low+(hi-low)*(1+2*(j-1))/4.0,{j,1,2}]]/(xscale)];plpS[[i]]=ListPlot[MovingAverage[aimargbin,{(2/20),(5/20),(6/20),(5/20),(2/20)}],PlotRange->{{low,hi},{0,1.1*Max[yai]}},Frame->True,FrameTicks->{xaxisticks,Automatic,Automatic,Automatic}, Axes->False,PlotStyle->{Black,Thickness[0.007],PointSize[0.05]},Joined->True,AspectRatio->0.5,FrameLabel->{Style[label,{"Helvetica-Bold",16}],
		Style["PDF",{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]
],{i,1,np}
];
];

The line below gives the actual area within the credible region for each parameter achieved for the particular choice of the number of bins, nai.

If[ncsys<n&&lsel>0&&aidiff<1,N[areai],]

The 4th and 6th items of strtpbin2 are the lower boundaries of the phase windows for χ and ω 
that lead to the minimum credible regions for these two items..

strtpbin2

meanai

medianai

modeai

minai

maxai

Next step does the plotting.

(*If[ncsys<n&&lsel>0&&aidiff<1,plmarg1a=Show[GraphicsGrid[{{plpS[[1]],plpS[[2]],plpS[[5]]},{plpS[[7]],,}},Spacings->{Scaled[-0.0],Scaled[0.05]}]]]*)

(*Export["c:/Users/Philip/Documents/greg8/47UMaNew/fig/APT3planGene3C47UMaNewOffsetsCorrHETHJSM7Marg.eps", plmarg1a,ImageSize->72*9]*)

If[ncsys<n&&lsel>0&&aidiff<1,plmarg1allS=Show[GraphicsGrid[If[EvenQ[np],Table[{plpS[[ii]],plpS[[ii+1]]},{ii,1,np,2}],Join[Table[{plpS[[ii]],plpS[[ii+1]]},{ii,1,IntegerPart[(np)/2]*2,2}],{{plpS[[IntegerPart[(np)/2]*2+1]]}}]],Spacings->{Scaled[-0.1],Scaled[0.05]}]]]

Display smoothed marginal distributions together with true parameter values.

Comparison of smoothed and unsmoothed marginals.

Show[plp[[5]],plpS[[5]]]

trueIntercept

ToString[trueIntercept]

plalphaLineS=Show[plpS[[1]],Graphics[{Thickness[0.005],Line[{{trueIntercept,0},{trueIntercept,6}}],Text[Style[ToString[trueIntercept],FontSize->11,FontWeight->"Bold"],{-0.08,0.15},{-1,1},{0,1}]}]]

trueSlope

plbetaLineS=Show[plpS[[2]],Graphics[{Thickness[0.005],Line[{{trueSlope,0},{trueSlope,5}}],Text[Style[ToString[Round[1000.*trueSlope]/1000.],FontSize->11,FontWeight->"Bold"],{0.83,0.2},{-1,1},{0,1}]}]]

trueSigma

plsigmaLineS=Show[plpS[[3]],Graphics[{Thickness[0.005],Line[{{trueSigma,0},{trueSigma,6}}],Text[Style[ToString[Round[1000.*trueSigma]/1000.],FontSize->11,FontWeight->"Bold"],{0.22,0.1},{-1,1},{0,1}]}]]

trueMu

ToString[trueMu]

plmuLineS=Show[plpS[[4]],Graphics[{Thickness[0.005],Line[{{trueMu,0},{trueMu,3}}],Text[Style[ToString[trueMu],FontSize->11,FontWeight->"Bold"],{-0.15,0.2},{-1,1},{0,1}]}]]

trueTau

pltauLineS=Show[plpS[[5]],Graphics[{Thickness[0.005],Line[{{trueTau,0},{trueTau,3}}],Text[Style[ToString[Round[1000.*trueTau]/1000.],FontSize->11,FontWeight->"Bold"],{0.97,0.3},{-1,1},{0,1}]}]]

plmarg3S=Show[GraphicsGrid[{{plalphaLineS,plbetaLineS},{plsigmaLineS}},Spacings->{Scaled[-0.1],Scaled[0.05]}]]

plmarg5S=Show[GraphicsGrid[{{plalphaLineS,plbetaLineS},{plsigmaLineS,plmuLineS},{pltauLineS}},Spacings->{Scaled[-0.15],Scaled[0.05]}]]

Most probable set of parameters & 68.3% marginal credible regions for most common peak.

Extract parameter iterations beyond burn-in with period values between pLow and pHi. 
By default pLow = avmin[[1]] and pHi = avMax[[1]], the search range.
Order the list of parameter iterations according to the joint probability density in descending order.

lsel=0;a1Low=-2.0;a1Hi=2.0;tr=Take[Transpose[logpdfβi][[k]],{nburn,nnt,1}];trPsel=Select[tr,#[[2]]>a1Low&&#[[2]]<a1Hi&];lsel=Length[trPsel];lsel=Length[trPsel];If[lsel>0,trs=RotateLeft[Sort[RotateRight[trPsel,{0,1}],OrderedQ[{#2,#1}]&],{0,1}],];lsel

List the most probable iteration vector.

If[lsel>0,trs1=RotateLeft[trs[[1]],1];ntrs=Length[trs1];trs[[1]]=RotateRight[Join[Table[If[type[[i]]==2,Mod[trs1[[i]],2*π],trs1[[i]]],{i,1,np}],Take[trs1,-(ntrs-np)]],1]]

The next command extracts the peak parameter values.

type

If[lsel>0,peakai=Table[0,{i,1,np}];probMax=trs[[1,nl]];
Do[
peakai[[i]]=trs[[1,i+1]],{i,1,np}
];
]

peakai

Tabulate the parameter peak values and 68% credible regions (relative to peak).

swn=1;
If[ncsys<n&&lsel>0&&aidiff<1,
l1={{"Item","Peak","Mean","Median","Mode"},{"","","","",""}};
l2=
Table[{
If[i≤np,
label=StringForm[name[[i]]<>unit[[i]]],
If
[i>np&&i≤np,j=i-np;label=(Subscript[ToExpression["χ"], j]),
j=i-np;label=(Subscript[ToExpression["ω"], j])
]
],
peakai[[i]],meanai[[i]],medianai[[i]],modeai[[i]]},
{i,1,np}
];
l4={{"","","","",""},{"Item","68% Cred Reg","68% Cred Reg","Peak-Min","Max-Peak"},{"","Min","Max","68%","68%"}};
l5=
Table[{
If[i≤np,
label=StringForm[name[[i]]<>unit[[i]]],
If
[i>np&&i≤np+nplan,j=i-np;label=(Subscript[ToExpression["χ"], j]),
j=i-np-nplan;label=(Subscript[ToExpression["ω"], j])
]
],
minai[[i]],maxai[[i]],peakai[[i]]-minai[[i]],maxai[[i]]-peakai[[i]]},
{i,np}
];
smj=MatrixForm[Join[l1,l2,l4,l5]],]

Tabulate the parameter peak values and 68% credible regions (relative to median).

swn=1;
If[ncsys<n&&lsel>0&&aidiff<1,
l1={{"Item","Peak","Mean","Median","Mode"},{"","","","",""}};
l2=
Table[{
If[i≤np,
label=StringForm[name[[i]]<>unit[[i]]],
If
[i>np&&i≤np,j=i-np;label=(Subscript[ToExpression["χ"], j]),
j=i-np;label=(Subscript[ToExpression["ω"], j])
]
],
peakai[[i]],meanai[[i]],medianai[[i]],modeai[[i]]},
{i,1,np}
];
l4={{"","","","",""},{"Item","68% Cred Reg","68% Cred Reg","Median-Min","Max-Median"},{"","Min","Max","68%","68%"}};
l5=
Table[{
If[i≤np,
label=StringForm[name[[i]]<>unit[[i]]],
If
[i>np&&i≤np+nplan,j=i-np;label=(Subscript[ToExpression["χ"], j]),
j=i-np-nplan;label=(Subscript[ToExpression["ω"], j])
]
],
minai[[i]],maxai[[i]],medianai[[i]]-minai[[i]],maxai[[i]]-medianai[[i]]},
{i,np}
];
smj=MatrixForm[Join[l1,l2,l4,l5]],]

Note: the peak values correspond to the parameter values of the MCMC iteration having the largest 
joint posterior probability. The peak is the peak of the binned parameter marginals. It can happen
that one of the + or - error bar values can be a small negative quantity as the credible regions are 
based on the binned marginals.

Compute Gelman-Rubin statistic for β = 1 MCMC chain for probability peak of interest.
Only applicable to parallel tempering involving multiple chains.

Notes on burn-in period and tests for convergence:

Here we employ the Gelman-Rubin statistic to decide on convergence. Normally, the MCMC test of
    convergence is to run repeated simulations with different starting points and compute the 
    Gelman-Rubin statistic for each parameter from the repeated runs as described below.
    
         Let θ represent one of the model parameters.
         Let (Subscript[θ, j]^i) represent the (i^th) iteration of the (j^th) simulation.
         Extract the last η post burn-in iterations for each simulation.

Mean within chain variance  W=(1/(m(η-1)))(∑_(j=1)^m)(∑_(i=1)^η)((((Subscript[θ, j]^i)-(OverBar[(Subscript[θ, j])])))^2)
               Between chain variance  B=(η/(m-1))(∑_(j=1)^m)((((OverBar[(Subscript[θ, j])])-(OverBar[θ])))^2)
               Estimated variance   (OverHat[V])(θ)=(1-(1/η))W+(1/η)B
               Gelman-Rubin statistic =(Sqrt[(((OverHat[V])(θ))/W)])
                       The Gelman-Rubin statistic should be close to 1.0  (e.g. <1.05)  for all paramaters for convergence.
          Ref:Gelman,A.and D.B.Rubin (1992)'Inference from iterative simulations    
                      using multiple sequences (with discussion)',Statistical Science 7, pp. 457-511.

In our parallel tempering MCMC, roughly every 2*(nthin = 10)*(nβ=8) = 160 iterations the β=1 iterations accepts a swap from its neigboring independent simulation some of which are to locations that are remote in parameter space. The final β=1 simulation is thus a combination of a very large number of independent simulations. We divide the final β=1 post burn-in iterations into typically 10 equal intervals and intercompare the ten different and essentially independent distributions for each parameter using the Gelman-Rubin test.

The starting point for this section is the list of FMCMC samples, trPsel, for the probability peak of interest obtained in step (30).

If[ncsys<n&&lsel>0,trn=Transpose[trPsel];ng=Length[trn[[1]]]]

Compute Gelmen-Rubin statistic for each parameter for the post burn-in results.

AbsoluteTiming[nparts=10;seg=ng;gelmenrubin=Table[iter=Take[If[type[[i]]==2||type[[i]]==4,Mod[trn[[i+1]],2*π],trn[[i+1]]],-seg];deln=Floor[seg/nparts];w=Table[subiter=Take[iter,{1+(j-1)*deln,(j-1)*deln+deln}];{Mean[subiter],Total[(((subiter-Mean[subiter]))^2)]/(deln-1)},{j,1,nparts}];tw=Transpose[w];means=tw[[1]];W=Total[tw[[2]]]/nparts;B=Total[(((means-Mean[iter]))^2)]*deln/(nparts-1);(Sqrt[(((1-1/deln)*W+B/deln)/W)]),{i,1,np}]]

Examine the autocorrelation of MCMC samples

The autocorrelation formula has the form:

(∑_common)((Subscript[x, i])-(OverBar[x]))((Subscript[x', i+j])-(OverBar[x])')/((Sqrt[(∑_i)((((Subscript[x, i])-(OverBar[x])))^2)])(Sqrt[(∑_i)((((Subscript[x, i]^('))-(OverBar[x])'))^2)]))

Set nc, the number of correlation lags (different values of j in the above equation) and dj the increment in j.

k=nβ;a1Low=-2.0;a1Hi=2.0;tr=Take[Transpose[logpdfβi][[k]],-(nnt-nburn)];trPsel=Select[tr,#[[2]]>a1Low&&#[[2]]<a1Hi&];trn=Transpose[trPsel];a1good=trn[[2]];a2good=trn[[3]];a3good=trn[[4]];a4good=trn[[5]];a5good=trn[[6]];prob=trn[[nl]];ng=Length[trn[[1]]]

nc=1000;dj=2;   (* nc must be less than nnt-nburn *)

If[nc>ng,nc=ng,nc]

Timing[plcorr=Table[0,{i,np}];(* initiate table of graphs to be plotted *)
Do[h1={};label=StringForm[name[[i]]<>unit[[i]]<>"  ACF"];Do[fga=Take[trn[[i+1]],-(ng-j)];fa=fga-Mean[fga];fgb=Take[trn[[i+1]],(ng-j)];fb=fgb-Mean[fgb];h1=Append[h1,{j,fa.fb/((Sqrt[fa.fa])*(Sqrt[fb.fb]))}],{j,0,nc,dj}];plcorr[[i]]=ListPlot[h1,Frame->True,Joined->True,Axes->{True,False},AspectRatio->0.5,PlotStyle->{Black,Thick},PlotRange->{{-10,nc},{-0.1,1}},PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Lag",{"Helvetica-Bold",12}],
		Style[label,{"Helvetica-Bold",12}]},FrameStyle->Thickness[0.003],BaseStyle->{12,FontFamily->"Helvetica-Bold"}],{i,1,np}]]

plIter=Show[GraphicsGrid[If[EvenQ[np],Table[{plcorr[[ii]],plcorr[[ii+1]]},{ii,1,np,2}],Join[Table[{plcorr[[ii]],plcorr[[ii+1]]},{ii,1,IntegerPart[(np)/2]*2,2}],{{plcorr[[IntegerPart[(np)/2]*2+1]]}}]],Spacings->{Scaled[-0.25],Scaled[-0.02]}]]

Compute MAP parameters using Nelder-Mead and compare with maxlogpdfAll values.

fm3[a1_?NumericQ,a2_?NumericQ,a3_?NumericQ,a4_?NumericQ,a5_?NumericQ]:=
Module[{x,y},
ζ={a1+a2*a4,a4};yx-ζ;
vmat=Table[{{(a2^2)*(a5^2)+(a3^2)+((yerrvec[[i]])^2),a2*(a5^2)},{a2*(a5^2),(a5^2)+((xerrvec[[i]])^2)}},{i,1,nd}];
logprior=c -Log[a3+mJslim]-Log[a5+mJslim]; (* the prior for parameters *)
loglike=loglikeoffset+logfact-0.5*(∑_(i=1)^nd)(Log[Det[vmat[[i]]]])-0.5*(∑_(i=1)^nd)((Transpose[yx-ζ][[i]].Inverse[vmat[[i]]].Transpose[{Transpose[yx-ζ][[i]]}])[[1]]);
   (* likelihood fn for particular parameters *)  
minusloglike=-loglike-logprior
]

Evaluate fm for a particular set of variables.

Sequence@@Take[maxlogpdfAll,{2,np+1}]

fm3[Sequence@@Take[maxlogpdfAll,{2,np+1}]]

bound=Table[{ToExpression["vb"<>ToString[i]],minai[[i]],maxai[[i]]},{i,1,np}]

ans=Timing[NMinimize[{fm3[vb1,vb2,vb3,vb4,vb5],{vb3≥0.0,vb5≥0.0}},bound,Method->"NelderMead",MaxIterations->4000]]

vlist=ans[[2,2]]

ann=Table[ToExpression["vb"<>ToString[i]],{i,1,np}]

parMAP=ann/.vlist

Evaluate fm3 for the best parameter set found from the hybrid MCMC.

bestMCMC=fm3[Sequence@@Take[maxlogpdfAll,{2,np+1}]]

Compare this value of fm3 found by Nelder-Mead to that found above for the best parameter set from the hybrid MCMC.
Better means a smaller value, i.e., deeper minimum.

bestNedlerMead=ans[[2,1]]

If Nelder-Mead does not find a more probable set of parameters (i.e., smaller -Log[prior × likelihood]), then execute following command to use the hybrid MCMC best values after re-ordering the parameters so lowest period parameters occur first, etc..

If[bestMCMC<bestNedlerMead,parMap=Take[maxlogpdfAll,{2,np+1}]]

The next command extracts the peak parameter values and adds χ and ω parameter values.

type

lsel=2;

If[lsel>0,peakMAP=Table[0,{i,1,np}];
Do[
peakMAP[[i]]=parMAP[[i]],{i,1,np}
];
]

peakMAP

Examine the correlation coefficient of simulated true data

The autocorrelation formula has the form:

(∑_(i=1)^nd)((Subscript[x, i])-(OverBar[x]))((Subscript[x', i+j])-(OverBar[x])')/((Sqrt[(∑_i)((((Subscript[x, i])-(OverBar[x])))^2)])(Sqrt[(∑_i)((((Subscript[x, i]^('))-(OverBar[x])'))^2)]))

xx=xtvec;mnxvec=Mean[xx];yy=ytvec;mnyvec=Mean[yy];num=(xx-mnxvec).(yy-mnyvec);den1=(Sqrt[(xx-mnxvec).(xx-mnxvec)]);
den2=(Sqrt[(yy-mnyvec).(yy-mnyvec)]);num/(den1*den2)

acTrue=Covariance[(xx-mnxvec),(yy-mnyvec)]/((Sqrt[Variance[(xx-mnxvec)]])*(Sqrt[Variance[(yy-mnyvec)]]))

Correlation coefficient from multiple simulations of the straight line based on unscaled true parameters, errors in both coordinates and intrinsic scatter, followed by rescaling of the simulated data.

nd=15;

A=0.0;B=0.5;sigmay=1.0;sigmax=1.0;sig=1.0;tau=7.0;

corCoefAll={};
Do[xtvecRaw1=Table[RandomReal[NormalDistribution[0.0,tau]],{nd}];ya= A + B*xtvecRaw1;ytvecRaw1=ya+Table[RandomReal[NormalDistribution[0,sig]],{nd}];
y1vec1=ytvecRaw1+Table[Random[NormalDistribution[0,sigmay]],{nd}];
x1vec1=xtvecRaw1+Table[Random[NormalDistribution[0,sigmax]],{nd}];
xx=x1vec1;mnxvec=Mean[xx];yy=y1vec1;mnyvec=Mean[yy];num=(xx-mnxvec).(yy-mnyvec);
corCoef1=Covariance[(xx-mnxvec),(yy-mnyvec)]/((Sqrt[Variance[(xx-mnxvec)]])*(Sqrt[Variance[(yy-mnyvec)]]));
corCoefAll=Append[corCoefAll,corCoef1],
{i,1,10000}]

plCorCoef1a=Histogram[corCoefAll,PlotRange->{{0.4,1},All},Frame->True,Axes->False,AspectRatio->0.5,ChartStyle->{Gray},PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Correlation coefficient",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

plCorCoefLine=Show[plCorCoef1a,Graphics[{Thickness[0.005],Line[{{acTrue,0},{acTrue,1900}}],Text[Style[ToString[Round[1000*acTrue]/1000.],FontSize->11,FontWeight->"Bold"],{acTrue,1500},{-1,1},{0,1}]}]]

plCorCoef1=Show[plCorCoef1a,Graphics[{{Text[Style["Based on true parameters",FontSize->11,FontWeight->"Bold"],{0.635,1100}]},{Thickness[0.005],Line[{{acTrue,0},{acTrue,1200}}],Text[Style[ToString[Round[1000*acTrue]/1000.],FontSize->11,FontWeight->"Bold"],{acTrue,1000},{-1,1},{0,1}]}}]]

Correlation coefficient from multiple simulations for MAP estimates of intercept and slope, intrinsic scatter, mu and tau.

nd=15;

MAP values are given by peakMAP

peakMAP

σx
σy

A=peakMAP[[1]];B=peakMAP[[2]];sig=peakMAP[[3]];mu=peakMAP[[4]];tau=peakMAP[[5]]

corCoefMAPAll={};
Do[xtvecRaw1=Table[RandomReal[NormalDistribution[mu,tau]],{nd}];ya= A + B*xtvecRaw1;ytvecRaw1=ya+Table[RandomReal[NormalDistribution[0,sig]],{nd}];
y1vec1=ytvecRaw1+Table[Random[NormalDistribution[0,σy]],{nd}];
x1vec1=xtvecRaw1+Table[Random[NormalDistribution[0,σx]],{nd}];
xx=x1vec1;mnxvec=Mean[xx];yy=y1vec1;mnyvec=Mean[yy];num=(xx-mnxvec).(yy-mnyvec);
corCoefE=Covariance[(xx-mnxvec),(yy-mnyvec)]/((Sqrt[Variance[(xx-mnxvec)]])*(Sqrt[Variance[(yy-mnyvec)]]));
corCoefMAPAll=Append[corCoefMAPAll,corCoefE],
{i,1,10000}]

corCoefMAPAll1=Select[corCoefMAPAll,#>0.4&];Length[corCoefMAPAll1]

plCorCoefMAP=Histogram[corCoefMAPAll1,PlotRange->{{0.4,1},All},Frame->True,Axes->False,AspectRatio->0.5,ChartStyle->{Gray},PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["Correlation coefficient",{"Helvetica-Bold",14}],
		Style["Number",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

plCorCoefMAPa=Show[plCorCoefMAP,Graphics[{{Text[Style["Based on MAP parameters",FontSize->11,FontWeight->"Bold"],{0.6,650}]},{Thickness[0.005],Line[{{acTrue,0},{acTrue,700}}],Text[Style[ToString[Round[1000*acTrue]/1000.],FontSize->11,FontWeight->"Bold"],{acTrue,400},{-1,1},{0,1}]}}]]

plACsimMAP=Show[GraphicsGrid[{{plCorCoef1},{plCorCoefMAPa}},Spacings->{Scaled[-0.0],Scaled[0.0]}]]

Plot individual data points versus time with overlaid light curve for MAP solution for dominant peak.

(*peakMAP={5.368364175126877`,0.016748463147030497`,-0.44905284343398477`,6.4398672380644255`,12.780994160558485`,15.260854869968625`,3.3840208617178886`,0.7268885960785179`,8.155876798407073`};*)

peakai

maxlogpdfAll

peakMAP

xmin=1.2*Min[xvec];xmax=1.1*Max[xvec];
ymin=peakMAP[[1]]+peakMAP[[2]]*xmin;
ymax=peakMAP[[1]]+peakMAP[[2]]*xmax;
pnts={{xmin,ymin},{xmax,ymax}};
plfit=ListPlot[pnts,Joined->True,PlotStyle->{Black,Dashed,Thickness[0.003]},DisplayFunction->Identity];
plfitdat=Show[plxy1a,plfit,DisplayFunction->$DisplayFunction]

xmin=1.2*Min[xvec];xmax=1.1*Max[xvec];
ymin=peakMAP[[1]]+peakMAP[[2]]*xmin;
ymax=peakMAP[[1]]+peakMAP[[2]]*xmax;
pnts={{xmin,ymin},{xmax,ymax}};
plfitMAP=ListPlot[pnts,Joined->True,PlotStyle->{Black,DotDashed,Thickness[0.006]}];
plfitdata=Show[plxy1,plfitMAP,DisplayFunction->$DisplayFunction]

plfitdatpk=Show[plfitdata,Graphics[{Text["(b)",{0,-2.2}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,DisplayFunction->Identity]

Plot residuals for MAP solution for dominant peak.

ymvec=peakMAP[[1]]+peakMAP[[2]]*xvec

yvec

y=yvec-ymvec

yy=y/(Sqrt[(yerrvec^2)])

reducedchisqu=yy.yy/(nd-np)

yyL=y/(Sqrt[(yerrvec^2)+((peakMAP[[3]])^2)+((peakMAP[[2]])^2)*(xerrvec^2)])

yyL.yyL/(nd-np)

nd

rmsres=(Sqrt[((y.y)/nd)])

resid=Transpose[{xvec,yvec-ymvec}]

respk=StandardDeviation[resid]

Mean[yerrvec]

(Sqrt[((respk[[2]])^2)-((Mean[yerrvec])^2)])

(*residerr=Table[{resid[[i]],ErrorBar[(Sqrt[((yerrvec[[i]])^2)+((peakMAP[[3]])^2)])]},{i,nd}];*)

residerr=Table[{resid[[i]],ErrorBar[yerrvec[[i]]]},{i,nd}];
plshp2=ErrorListPlot[residerr,Frame->True,PlotRange->All,Axes->{True,False},PlotStyle->Black,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["x",{"Helvetica-Bold",14}],Style["y",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

plshp2b=Show[plshp2,Graphics[{Text["(c)",{0,-1.}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,DisplayFunction->Identity]

plboth=Show[GraphicsGrid[{{plxy1a},{plfitdatpk},{plshp2b}},Spacings->{Scaled[-0.0],Scaled[0.00]}]]

Plot data and fit residuals with overlaid mean light curve.

If[ncsys<n,k=nβ;a1Low=-2.0;a1Hi=2.0;r=20;tr=Take[Transpose[logpdfβi][[nβ]],{nburn,nnt,r}];xsel=Select[tr,#[[2]]>a1Low&&#[[2]]<a1Hi&];nsel=Length[xsel]
]

xsel[[1]]

modeai

xvec

mvel={};xlist=Table[Min[xvec]+(i-1)*(Max[xvec]-Min[xvec])/100,{i,1,101}];
Do[str=Take[xsel[[i]],{2,np+1}];
mveli=str[[1]]+str[[2]]*xlist;
mvel=Append[mvel,mveli]
,{i,1,nsel}
]

meany=Mean[mvel];

stdy=StandardDeviation[mvel];

Length[(meany+stdy)]

modlist=Transpose[{xlist,meany}];plfitmean=ListPlot[modlist,Joined->True,PlotStyle->{Black,Thickness[0.005]}];
plfitdatamean=Show[plxy1A,plfitmean]

modlist1=Transpose[{xlist,(meany+stdy)}];modlist2=Transpose[{xlist,(meany-stdy)}];plfitmeanr=ListPlot[{modlist1,modlist2},Joined->True,PlotStyle->{{Black,Opacity[0.7],Dashed,Thickness[0.005]},{Black,Opacity[0.7],Dashed,Thickness[0.005]}}];
plfitdatameanr=Show[plxy1A,plfitmeanr,plfitmean,plfitMAP]

(*Export["c:/Users/Philip/Documents/greg13/MCMCregression/fig/CorStLineFitIntrinsicScatterYesM7informativeGaussXtPriorAnalyticFitTrue.eps", plfitdatameanr,ImageSize->72*7]*)

Show[plxy1A,plfitmeanr,plfitmean,plfitMAP]

plmapfitmeana=Show[plfitdatameanr,Graphics[{Text["(a)",{0,-2.2}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,PlotRange->{All,{-2.6,2.6}}]

plmapfitmeanb=Show[plfitdatameanr,Graphics[{Text["(b)",{0,-2.2}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,PlotRange->{All,{-2.6,2.6}}]

(*mvec={};mres={};myerr={};
Do[str=Take[xsel[[i]],{2,np+1}];
mveci=str[[1]]+str[[2]]*xvec;
mresi=(mveci-yvec);
myerri=(Sqrt[(yerrvec^2)+((str[[3]])^2)]);
mvec=Append[mvec,mveci];
mres=Append[mres,mresi];
myerr=Append[myerr,myerri]
,{i,1,nsel}
]*)

mvec={};mres={};myerr={};
Do[str=Take[xsel[[i]],{2,np+1}];
mveci=str[[1]]+str[[2]]*xvec;
mresi=(yvec-mveci);
myerri=yerrvec;
mvec=Append[mvec,mveci];
mres=Append[mres,mresi];
myerr=Append[myerr,myerri]
,{i,1,nsel}
]

meanmvec=Mean[mvec]

meanres=Mean[mres]

meanyerr=Mean[myerr]

stdmvec=StandardDeviation[mvec]

y=meanres/meanyerr

reducedchi2=y.y/(nd-np)

StandardDeviation[mres]

meanres

residerr=Table[{{xvec[[i]],meanres[[i]]},ErrorBar[meanyerr[[i]]]},{i,nd}]

plshp2mean=ErrorListPlot[residerr,Frame->True,PlotRange->All,Axes->{True,False},PlotStyle->Black,PlotLabel->Style["",{"Helvetica-Bold",12}],FrameLabel->{Style["x",{"Helvetica-Bold",14}],Style["Residual",{"Helvetica-Bold",14}]},FrameStyle->Thickness[0.006],BaseStyle->{12,FontFamily->"Helvetica-Bold"}]

plshp2meanb=Show[plshp2mean,Graphics[{Text["(c)",{0,-1.1}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5]

plThree=Show[GraphicsGrid[{{plxy1a},{plmapfitmeanb},{plshp2meanb}},Spacings->{Scaled[-0.1],Scaled[-0.02]}]]

plboth2=Show[GraphicsGrid[{{plxy1a,plmapfitmeanb}},Spacings->{Scaled[-0.25],Scaled[-0.02]}]]

Figure caption: FMCMC analysis which includes the hidden variables as parameters. Panel (a) shows the true regression line (gray), the measured values with error bars. Panel (b) shows the true regression line (gray), the measured values with error bars, the mean fitted line (black) and MAP fit  (dot-dashed) derived from the FMCMC iterations, and the mean fitted line ±1 standard deviation in the FMCMC fit uncertainty (dashed curves).

plRawfit=Show[{plxtyt,plxyline,plfitMAP},Graphics[{Text["(d)",{0,-2.2}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,PlotRange->{All,{-2.6,2.6}}]

plRawfitb=Show[{plxtyt,plxyline,plfitMAP,plfitmeanr,plfitmean},Graphics[{Text["(b)",{0,-2.2}]}],TextStyle->{FontWeight->"Bold",FontSize->12},AspectRatio->0.5,PlotRange->{All,{-2.6,2.6}}]

plFour=Show[GraphicsGrid[{{plxy1a,plmapfitmeanb},{plshp2meanb,plRawfit}},Spacings->{Scaled[-0.2],Scaled[-0.02]}]]

Figure caption: FMCMC analysis which includes the hidden variables as parameters. Panel (a) shows the true regression line (gray), the measured values with error bars. Panel (b) shows the true regression line (gray), the measured values with error bars, the mean fitted line (black) and MAP fit  (dot-dashed) derived from the FMCMC iterations, and the mean fitted line ±1 standard deviation in the FMCMC fit uncertainty (dashed curves). Panel (c) shows the residuals from the mean fit. Panel (d) shows the true regression line (gray), the MAP fit  (dot-dashed) derived from the FMCMC iterations, and the true coordinates that include intrinsic scatter but not measurement uncertainties (black dots).

plTwo=Show[GraphicsGrid[{{plmapfitmeana,plRawfitb}},Spacings->{Scaled[-0.25],Scaled[-0.02]}]]

Figure caption: FMCMC analysis which includes the hidden variables as parameters. Panel (a) shows the true regression line (gray), the measured values with error bars, the mean fitted line (black) and MAP fit  (dot-dashed) derived from the FMCMC iterations, and the mean fitted line ±1 standard deviation in the FMCMC fit uncertainty (dashed curves). Panel (b) shows the true regression line (gray), the mean fitted line (black), the MAP fit  (dot-dashed) derived from the FMCMC iterations, the mean fitted line ±1 standard deviation in the FMCMC fit uncertainty (dashed curves), and the true coordinates that include intrinsic scatter but not measurement uncertainties (black dots).

NotebookSave[]