Undergrad Efficient Monte Carlo Fitting with Python Package Emcee and MCMC Analysis

Click For Summary
The discussion focuses on using the Python package emcee for fitting functions to data points through MCMC analysis. The user experiences issues with parameter plots showing a rapid convergence and an unusual initial shape, suspecting it may be due to random walks affecting parameter errors. Suggestions include discarding the initial samples, known as the burn-in period, to improve results. It is noted that the emcee package may not automatically account for burn-in, leading to the observed issues. Properly configuring the sampler and resetting it after initial runs can help mitigate these problems.
Silviu
Messages
612
Reaction score
11
Hello I use a python package called emcee to fit a function to some data points. The fit looks great, but when I want to plot the value of each parameter at each step I get the plot.png below. In their example (with a different function and data points) they get the plot1.png image. Why is my function converging so fast, and why does it have that weird shape in the beginning. I apply MCMC using likelihood and posterior probability. And even if the fit looks very good, the error on parameters of function are very small (10^10 smaller than the actual value) and I think it is because of the random walks. Any idea how to fix it?
 

Attachments

  • plot.png
    plot.png
    9.1 KB · Views: 538
  • plot1.png
    plot1.png
    56 KB · Views: 564
Physics news on Phys.org
That looks like it is just the burn in period. Usually people just discard the first few thousand points in order to avoid that.

I have done all of my MCMC work in R, so I am not familiar with the Python package you mention, but it looks like it is either set with a burn in of 0 or it plots the burn in samples.
 
If you're using emcee something like this will help overcome the burn in:
sampler = EnsembleSampler(...)

pos, prob, state = sampler.run_mcmc(...)

sampler.reset()

sampler.run_mcmc(pos, 1000, rstate0=state)
 
The standard _A " operator" maps a Null Hypothesis Ho into a decision set { Do not reject:=1 and reject :=0}. In this sense ( HA)_A , makes no sense. Since H0, HA aren't exhaustive, can we find an alternative operator, _A' , so that ( H_A)_A' makes sense? Isn't Pearson Neyman related to this? Hope I'm making sense. Edit: I was motivated by a superficial similarity of the idea with double transposition of matrices M, with ## (M^{T})^{T}=M##, and just wanted to see if it made sense to talk...

Similar threads

  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 12 ·
Replies
12
Views
3K
  • · Replies 2 ·
Replies
2
Views
1K
Replies
8
Views
2K
  • · Replies 16 ·
Replies
16
Views
3K
  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 1 ·
Replies
1
Views
1K
Replies
28
Views
4K
  • · Replies 1 ·
Replies
1
Views
2K
Replies
1
Views
2K