chi square value for each chain

[YuichiHarikane]

Hi,

I have two questions.

1. Is there any way to obtain chi-square values of all chains in the fitting? Are the chi-square values already stored in the output h5 file?

2. I'm using the parametric sfh model and getting the star formation rate with

   res, obs, model = reader.results_from('my_h5_file')
   sps = reader.get_sps(res)
   model.predict(best["parameter"], obs=obs, sps=sps)
   theta_best_dict=model.params
   sfr = prospect.plotting.sfh.parametric_sfr(times=0.05, **theta_best_dict)

   Can I also get 16%, 50%, and 84% percentiles of the SFR?

Thank you in advance,

Yuichi

[bd-j]

Hi Yuichi,

The result dictionary has the keys "lnprobability" and, if using a nested sampler, "lnlikelihood" which are arrays with the these quantities for each sample. For the default noise models lnlikelihood is equal to -0.5 * chi_square - 0.5 * np.log(np.product(2 * pi * sigma**2)) where sigma is the uncertainty vector for your data, while lnprobability has an additional term from the ln of the prior probability.

To compute percentiles of the SFR I would probably do something like

from prospect.io.write_results import chain_to_struct
# get all the posterior samples in a structured array
chaincat = chain_to_struct(res["chain"].reshape(-1, model.ndim), model)
# make it a dictionary 
chaindict = {k:chaincat[k] for k in chaincat.dtype.names}
# get the sfr of every sample
sfrs = prospect.plotting.sfh.parametric_sfr(times=0.05, **chaindict)
# compute (weighted) quantile from all the sfrs
from prospect.plotting.corner import quantile
sfr_ptiles = quantile(sfrs[None, :], [0.16, 0.5, 0.84], weights=res.get("weights", None))[0]

[YuichiHarikane]

Hi Ben,

Thanks for your answers!

As for the percentiles of the SFR, I can successfully get them by changing the line 263 of prospect/plotting/sfh.py if tmass < 0: to if tmass.any() < 0: because I got the following error:

File "make_summary_plot.py", line 143, in main
    sfrs = prospect.plotting.sfh.parametric_sfr(times=0.05, **chaindict)
  File "/Users/hari/opt/anaconda3/envs/astroconda/lib/python3.7/site-packages/prospect/plotting/sfh.py", line 128, in parametric_sfr
    sfr, mass = compute_mass_formed(tage - times, pset)
  File "/Users/hari/opt/anaconda3/envs/astroconda/lib/python3.7/site-packages/prospect/plotting/sfh.py", line 263, in compute_mass_formed
    if tmass < 0:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

As for the chi-square value, I'm using emcee so I only have lnprobability. Is there any way to get the chi-square values for all chains? I could calculate model fluxes in photometric bands by running model.predict for each chain but it may take a long time. I found model.predict_phot but I'm not sure how to use this...

I have another question. In the parametric sfh fitting, I want to keep tage less than the cosmic age. I found prospect.models.transforms.tage_from_tuniv() written in the FAQ, but could you let me know how to use this? Does the following setting make sense? How can I add the tage_tuniv parameter?

model_params = TemplateLibrary["parametric_sfh"]
model_params["tage"]["init"] = 0.1
model_params["tage"]["init_disp"] = 1
model_params["tage"]["isFree"] = False
model_params["tage"]["depend"] = transforms.tage_from_tuniv

Thank you in advance.

Yuichi

[bd-j]

I pushed a fix to the main branch for the bug you noticed (though one would want to do (tmass < 0).any() or, more generally, np.any(tmass < 0) to make sure the inequality is evaluated before the call to any)

Re the chi-square, you could make new predicted SEDs via model.predict() and then use prospect.fitting.lnprobfn with appropriate flags, but it's probably faster to compute the prior probabilities and subtract them from the lnprobability vector. E.g.

lnp_prior = model.prior_product(res["chain"])
chisquare = res["lnprobability"] - lnp_prior

Also just FYI note that predict_phot() would not be faster than predict(); the spectrum generation that is the time consuming part is required anyway before the filter projections for photometric prediction.

[bd-j]

and re ages, that's almost right. But you need to introduce the new free parameter "tage_tuniv" which is the ratio tage/t_univ. Also dependencies are specified with "depends_on" key, not "depends"

model_params = TemplateLibrary["parametric_sfh"]
model_params["tage"]["isFree"] = False
model_params["tage"]["depends_on"] = transforms.tage_from_tuniv
# add new parameter
model_params["tage_tuniv"] = dict(N=1, isFree=True, prior=priors.Uniform(mini=0, maxi=1), init=0.5)

Note the initial value and dispersion of "tage" is now irrelevant, as it will be determined from the redshift (i.e. t_univ) and the new parameter "tage_tuniv"

[YuichiHarikane]

Hi Ben,

Thank you for your replies! I can successfully got the chi-square values for all chains.

I tried the fitting using your suggestions of tage_from_tuniv and tage_tuniv, and got the following outputs:

<class 'prospect.models.sedmodel.SpecModel'>

Free Parameters: (name: prior) 
-----------
  zred: <class 'prospect.models.priors.Uniform'>(mini=0.0,maxi=20.0)
  mass: <class 'prospect.models.priors.LogUniform'>(mini=100000000.0,maxi=100000000000.0)
  logzsol: <class 'prospect.models.priors.TopHat'>(mini=-2,maxi=0.4)
  dust2: <class 'prospect.models.priors.TopHat'>(mini=0.0,maxi=2.0)
  tage: <class 'prospect.models.priors.TopHat'>(mini=0.001,maxi=13.8)
  tau: <class 'prospect.models.priors.LogUniform'>(mini=0.1,maxi=10)

Fixed Parameters: (name: value [, depends_on]) 
-----------
  sfh: [4] 
  imf_type: [1] 
  dust_type: [2] 
  tage_tuniv: [0.5] 
  mass_units: ['mstar'] 
  add_igm_absorption: [ True] 
  add_neb_emission: [ True] 
  add_neb_continuum: [ True] 
  nebemlineinspec: [False] 
  gas_logz: [0.] <function stellar_logzsol at 0x7fcb7a48b560>
  gas_logu: [-2.] 
  elines_to_ignore: ['Ly alpha 1216'] 
number of walkers=118

It seems that the tage is still a free parameter and tage_univ is a fixed parameter (and the initial value is set to 0.5), but did I do something wrong? The obtained chain also looks strange. Although I'm fitting the redshift of zred~8, the tage ranges up to 5 (it should be smaller than t_universe=~0.6 Gyr at z~8), as shown in the attached plot (x-axis is zred and y-axis is tage).

Please let me know if you find something wrong.

Best wishes,

Yuichi

[bd-j]

Hi Yuichi,

There was a typo; "isfree" should be all lowercase. For me

from prospect.models import SpecModel, transforms, templates, priors

model_params = templates.TemplateLibrary["parametric_sfh"]
model_params["tage"]["isfree"] = False
model_params["tage"]["depends_on"] = transforms.tage_from_tuniv
# add new parameter
model_params["tage_tuniv"] = dict(N=1, isfree=True, prior=priors.Uniform(mini=0, maxi=1), init=0.5)

model = SpecModel(model_params)
print(model)

yields

:::::::
<class 'prospect.models.sedmodel.SpecModel'>

Free Parameters: (name: prior) 
-----------
  mass: <class 'prospect.models.priors.LogUniform'>(mini=100000000.0,maxi=1000000000000.0)
  logzsol: <class 'prospect.models.priors.TopHat'>(mini=-2,maxi=0.19)
  dust2: <class 'prospect.models.priors.TopHat'>(mini=0.0,maxi=2.0)
  tau: <class 'prospect.models.priors.LogUniform'>(mini=0.1,maxi=30)
  tage_tuniv: <class 'prospect.models.priors.Uniform'>(mini=0,maxi=1)

Fixed Parameters: (name: value [, depends_on]) 
-----------
  zred: [0.1] 
  sfh: [4] 
  tage: [1] <function tage_from_tuniv at 0x1aecd4310>
  imf_type: [2] 
  dust_type: [0] 

[bd-j]

note that because tau is no longer a free parameter it will not be logged in the output chain the same way. Thus, to make prospect.plotting.sfh.parametric_sfr() work you will have to construct the value of tau from the sampled parameters and zred, and then add that as a keyword argument to parametric_sfr()

e.g.

from prospect.models import transforms
tau = transforms.tage_from_tuniv(zred=model.params["zred"], tage_tuniv=chaindict["tage_tuniv"])
sfrs = prospect.plotting.sfh.parametric_sfr(times=0.05, tau=tau, **chaindict)

[YuichiHarikane]

Hi Ben,

Thanks. I can successfully run with tage_from_tuniv and tage_tuniv after correcting the typo.

Just to confirm, not tau but tage is no longer a free parameter so

from prospect.models import transforms
tage = transforms.tage_from_tuniv(zred=chaindict["zred"], tage_tuniv=chaindict["tage_tuniv"])
sfrs = prospect.plotting.sfh.parametric_sfr(times=0.05, tage=tage, **chaindict)

is that correct?

Yuichi

[bd-j]

yes, of course, tage has to be reconstructed, not tau 🤦

[YuichiHarikane]

Thanks!