SpectrumModel log-likelihood not showing

[mileslucas]

Hi, quick follow up question. I wonder if the same thing might be happening when I train the SpectrumModel object? For instance, when I now train my model (with a converged/trained emulator) to do a MAP fit following the example, I still get log likelihood as None even after training.

Originally posted by @iyera22 in https://github.com/iancze/Starfish/issues/134#issuecomment-790234684

[mileslucas]

@iyera22 would you mind posting a snippet of the code you're working with? You don't need to share any data, yet, I just want to get an idea of what we're working with.

[aishaiyer]

Sure, first I construct the Spectrum Model object:

model = SpectrumModel(
    "mygrid_emu.hdf5",
    'data.hdf5',
     grid_params=[2700, 5.0, 0.0, 0.5],
)

then, i define a set of priors:

priors = {
    "Teff": st.uniform(2000, 2000),
    "logg": st.uniform(4.0, 1.25),
    "logZ": st.uniform(-1., 2.),
    "CtoO": st.uniform(0.3,0.6),
}

Then, I train the model object:

model.train(priors)

Now when I print the model, I get log likelihood: None

SpectrumModel
-------------
Data: data
Emulator: mygrid_emu
Log Likelihood: None

Parameters
  Teff: 2597.239758212062
  logg: 5.249999999925803
  logZ: 0.0005916055749307812
  CtoO: 0.37563832498034355

[mileslucas]

Okay, can you tell me what outputs you get from the following:

model._lnprob
model.log_likelihood()
model._lnprob

you can also look inside model.residuals after calling train, it'll be a list of spectrum residuals, you can check if they look valid.

[aishaiyer]

ok so here are the following outputs:

(Pdb) model._lnprob
(Pdb) model.log_likelihood()
2152.916379999185
(Pdb) model._lnprob
2152.916379999185

I do get reasonable residuals however this brings me back to another concern that the emcee run doesn't seem to converge and gives flat posteriors for all fitted parameters using such a model object.

[aishaiyer]

Basically Im wondering if this has any effect on the emcee run if it's only a display problem and its actually behind the scene calculating the log likelihood?

[mileslucas]

Okay, the output there gives me confidence everything is being computed alright. For your info, when you instantiate a model and print it, it will say "None" for the log-likelihood, because we don't want to spend time calculating the log-likelihood just for printing. The first time you call log_likelihood, though, it will update and print appropriately, like you can see in your debugging.

If you're having trouble with convergence with emcee, that's an issue related with your data, priors, and your sampling parameters (number of walkers, number samples, etc.). There's no "one size fits all" fix for that, however the emcee documentation has some good FAQs as a start.

[aishaiyer]

Great, thank you! This helps!

[mileslucas]

Okay, I'm going to go ahead and close this issue- if something seems to be broken please don't hesitate to open a new one!

[aishaiyer]

Thank you!

[aishaiyer]

Hi! a quick update on the emcee convergence issue. I have the following situation: 1) I am working with flux-calibrated data 2) I have a stellar model with flux computed at the surface of the star 3) I am trying to fit for the log_scale parameter and use that to calculate the radius

The problem however is that within the spectrum model code, the log_scale is scaling the normalized emulator interpolated spectrum (mean flux) to match the flux calibrated data flux, and so the best log_scale doesn't fall within physically realistic radius values.

        # Only rescale flux_mean and flux_std
        if "log_scale" in self.params:
            scale = np.exp(self.params["log_scale"])
            fluxes[-2:] = rescale(fluxes[-2:], scale)

Additionally, say for example if my raw model flux (computed at the surface of the star) is about 10-20 dex higher than the flux calibrated data, and I use an optimum log_scale that visually works well to match the two, then the emcee posteriors for all best fit parameters appeared flat. The same happens when I do not fit for log_scale and simply let it autoscale (the following code). I understand this huge scaling offset between data and model is contributing massively to the covariance matrix.

        # Renorm to data flux if no "log_scale" provided
        if "log_scale" not in self.params:
            factor = _get_renorm_factor(self.data.wave, flux, self.data.flux)
            flux = rescale(flux, factor)
            X = rescale(X, factor)

The only situation where I get well constrained posteriors for all parameters is when I scale up the flux calibrated data to roughly match the level of the stellar model flux at the surface of the star before performing the fit.

This is an example of the retrieved Teff, both when I fit for a visually good log_scale with appropriate priors and also when I do not fit for log_scale at all and let it autoscale.

This is an example of the same retrieved Teff, where I again fit for a visually good log_scale and also when I do not fit for log_scale at all, only here I've modified the input flux calibrated data to roughly match the level of the model flux at the surface of the star. Of course this is not the desired way to do this if I'm looking to properly solve for the radius of the object.

I'm wondering if this version of the code stores the flux scaling value somewhere before everything gets normalized inside the emulator? Maybe I could call that value and use it within the spectrum model code to specifically fit for a log((R/d)^2)? Thanks!