Switch to EMCEE for determining inital DPI1 and eps_g

[nielsenmb]

In the l=1 models that include mixed modes, PBjam uses Dynesty to sample the prior conditional on a few input parameters in order to get initial guesses for DPi1 and eps_g for the n_g list.

For SG stars this is reasonably quick (~10s-1min), but for RG stars this can for some reason be pretty slow. It is unclear exactly what makes the process so much slower. However, using a nested sampler here is probably overkill, and EMCEE may be able to make this step somewhat faster.

[redeaglekr]

Hey, I tried using samplers.EmceeSampling on l20 models from Asyl20model(). It spits out nan values in log_likelihood values, which arises from the nan values of parameter values passed to likelihood fuction, I guess.

[nielsenmb]

Hi @redeaglekr can you show an example of the function calls that produce this error? If this is the peakbag stage then it's a bit strange that it produces NaN values if the inputs to the class are correct.

[redeaglekr]

Hi,

I first ran the l20model on with nested sampling, it went fine. Then, I changed the argument to samplers.EmceeSampling Asyl20model() and replaced line 105 of runl20model to self.l20Samples = self.l20model.runSampler(nsamples = 5000). The result of the run is

Burning in sampler
Parameters: [  0.06371292  -0.18220788  -0.07618848   0.20307707  -0.03551503
  -0.03079507 -15.94589277   3.80291599   0.04679628   3.14416928
   2.16587044   1.07442912]
lnlike: nan
lnprior: -inf
Traceback (most recent call last):
  File "/home/PBjam_files/for_keplerstar.py", line 52, in <module>
    self.modeID.runl20model()
  File "/home/pbjam_env/lib/python3.12/site-packages/pbjam/modeID.py", line 112, in runl20model
    self.l20Samples = self.l20model.runSampler(nsamples = 5000)
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/pbjam_env/lib/python3.12/site-packages/pbjam/samplers.py", line 546, in runSampler
    burnSampler, ctrlChain = self.burnInSampler(nchains, DEfrac, earlyStop, walltime, checkEvery, t0, ncontrol, DtauLimit)
                             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/pbjam_env/lib/python3.12/site-packages/pbjam/samplers.py", line 246, in burnInSampler
    p0 = self.initSamples(nchains)
         ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/homepbjam_env/lib/python3.12/site-packages/pbjam/samplers.py", line 115, in initSamples
    print('add terms:', self.AddLikeTerms(thetaU))
                        ^^^^^^^^^^^^^^^^^
AttributeError: 'Asyl20model' object has no attribute 'AddLikeTerms'

First, I couldnt find the function AddLikeTerms in jar. So I changed AddLikeTerms function to addAddObsLike from jar. I am not sure btw..

the result was

Burning in sampler
Parameters: [ -0.22196413  -0.17355245  -0.06878502   0.21866154  -0.14983814
   0.13674455 -15.92010177   3.75139231   0.02953804   3.16675585
   2.18725809   1.09632213]
lnlike: nan
lnprior: -inf
add terms: -18.241743232572368
Parameters: [-4.11537835e-01 -1.65858730e-01 -6.83726373e-02  6.51680865e-02
 -5.75043670e-03  1.84253456e-01 -1.64034910e+01  3.84220724e+00
 -5.95988954e-02  3.09818312e+00  2.18780170e+00  1.07835030e+00]
lnlike: nan
lnprior: -inf
add terms: -18.224371000340444
Parameters: [-1.48241017e-01 -3.36082237e-01 -9.10665050e-02 -1.15635700e-02
 -9.09854142e-02  1.57849232e-01 -1.58774074e+01  3.87055857e+00
  8.57542928e-02  3.40460764e+00  2.16392947e+00  1.01914599e+00]
lnlike: nan
lnprior: -inf
add terms: -17.991404041882642

and this goes on

/home/pbjam_env/lib/python3.12/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract
  lnpdiff = f + nlp - state.log_prob[j]
/home/pbjam_env/lib/python3.12/site-packages/emcee/autocorr.py:38: RuntimeWarning: invalid value encountered in divide
  acf /= acf[0]
Convergence >> 1: nan
and this goes on forever.. 

That being said, my whole purpose to switch to emcee was to make the samplinng quicker. ie, instead of sampling the full prior space, I was looking for a sort of initial best fit and then move on to dynesty, once we are sure the fit looks better.

[nielsenmb]

The missing AddLikeTerms is because the emcee sampling class we set up was only meant to be run within the peakbagging module. In the modeID module this attribute isn't set. I guess replacing it with the observational parameter likelihoods would be OK but it's not something I've thought about. I'm not entirely sure why it would lead to NaN values though :(

Dynesty was chosen for the modeID stage because it's where multi-modal posteriors tend to appear due to the prior volume potentially encompassing multiple mode peaks. emcee doesn't tend to be very good with these situations and can easily get stuck in one local maximum of the posterior. Once near the maximum however things tend to be much more uni-modal in which case emcee would probably be fine.

In my experience Dynesty tends to be pretty good at reducing the prior volume very quickly initially. You'll notice that the Delta_logZ value in the progress bar will be large initially indicating the prior volume is being reduce very rapidly. Once it becomes small however you'll see it takes a long time to converge to the target threshold.

So I think if you want to speed up the process it might be worth investigating doing it the other way around, i.e., use Dynesty to get close to the maximum and then use emcee the rest of the way when the posterior region that needs to be explored probably won't be multi-modal. This, at least the philosophy we adopted when splitting PBjam into the modeID and peakbagging stages. But it could be worth trying a different target Delta_logz in Dynesty to stop the sampling earlier, and then inspecting the results. It's not something I've looked into that much, but it could be it's already a very good solution that peakbagging can work with when it reaches Delta_logZ <= 10 or something like that.