1
$\begingroup$

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

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

Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ 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. $\endgroup$ Commented May 29, 2017 at 3:25
  • $\begingroup$ I added more to my question. $\endgroup$ Commented May 29, 2017 at 4:19
  • $\begingroup$ I'm asking how I would be able to get the 68% confidence interval for my parameters. This is the info that I know. $\endgroup$ Commented May 29, 2017 at 4:20
  • 2
    $\begingroup$ 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. $\endgroup$ Commented May 29, 2017 at 5:08
  • $\begingroup$ Alright I will try to post the code. $\endgroup$ Commented May 29, 2017 at 16:38

1 Answer 1

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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"? $\endgroup$ Commented Jul 7, 2018 at 22:47

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.