# Deriving 68% confidence level for each parameter after MCMC

I used the MCMC package https://github.com/joshburkart/mathematica-mcmc to find the best parameters for model of Massive Gravity. Now I need to find the confidence 68% for each of those parameters and have no idea how to proceed. This is my first time doing anything like this. https://arxiv.org/pdf/1205.1613v1.pdf I'm trying to reproduce the confidence levels in this paper which are in results section. Any help will be much appreciated. The current output of my MCMC

The current output of my MCMC MCMC[(1.725- RR[M, W, w, h, m, 1091.3])) (-1/2), {{M, .28, .001, Real}, {w, -1, .001, Real}, {W, 0, .001, Real}, {h, .7, .001, Real}, {m, .135, .001, Real}}, 39000]

"MCMCResult"[{M -> 0.292338, w -> -1.03469, W -> 0.0506029,h -> 0.680667,m -> 0.140849}, "\[LeftSkeleton]39000\[RightSkeleton]"]

"ParameterErrors" -> {M -> 0.0136284, w -> 0.0564952, W -> 0.0418342,

h -> 0.0131136, m -> 0.00542007}, "AverageAcceptance" -> 0.789823, "TimeSpent" -> 149855.631056 Second, "NumSteps" -> 39000, "Parameters" -> {M, w, W, h, m}, "ProposalSpreads" -> {0.001, 0.001, 0.001, 0.001, 0.001}, "ParameterDomains" -> {Real, Real, Real, Real, Real}, "BurnFraction" -> 0.1, "BurnEnd" -> 3900

• If you include your code you used to estimate the parameters (or at least a minimal working example), you'll likely get better and quicker help.
– JimB
May 29 '17 at 3:25
• I added more to my question. May 29 '17 at 4:19
• I'm asking how I would be able to get the 68% confidence interval for my parameters. This is the info that I know. May 29 '17 at 4:20
• Please add code as text, not (only) as image. You need to add all necessary input to reproduce your resuts, otherwise this question will be closed. May 29 '17 at 5:08
• Alright I will try to post the code. May 29 '17 at 16:38

This is an extended comment rather than an answer.

Without the definition of the function RR, one can't give you very specific advice.

However, if the output of MCMC is placed in mcmc, then you can get a list of all of the MCMC runs for all of the parameters using

results = mcmc["ParameterRun"];

This array will include all of the burn-in values. By looking at the code and trying a few other examples it appears that the runs that you want to keep are given by the following:

numsteps = 3900 (* The number of runs you requested *)
k = Min[1000, Floor[nnumsteps/2]];
results = mcmc["ParameterRun"][[k;;]]

From results you can duplicate mcmc["BestFitParameters"] and mcmc["ParameterErrors"]:

mcmc["BestFitParameters"]
{Mean[results[[All, 1]]], Mean[results[[All, 2]]]}

mcmc["ParameterErrors"]
{StandardDeviation[results[[All, 1]]], StandardDeviation[results[[All, 2]]]}

If you are able to assume (rather than just willing to assume) that the resulting distributions are normal, then parameter estimates plus-or-minus the associated standard deviation gets the "68% confidence intervals". (If you're really making Bayesian inferences, you should really call these "credible intervals" and explicitly state your prior distributions.)

If you can't assume normal distributions, then one can use the Quantile function to get the intervals:

Table[{Quantile[results[[All, i]], (100 - 68)/200],
Quantile[results[[All, i]], 1 - (100 - 68)/200]}, {i, Dimensions[results][[2]]}]

This finds the middle 68% of the results for each parameter.

• this last expression gives you an interval at whose extrema the distribution takes different values. How would you find the so-called "Highest Density Interval"? Jul 7 '18 at 22:47