Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Using prospector #166

Open
steveschulze opened this issue Mar 2, 2020 · 18 comments
Open

Using prospector #166

steveschulze opened this issue Mar 2, 2020 · 18 comments
Labels

Comments

@steveschulze
Copy link

Hi,

I have a few questions in using prospector.

  1. In the file demo_params.py in the directory demo is written that prospector uses the WMAP9 cosmology. Can the user change that?

  2. Galaxy mass: I would like to extract the mass currently locked in stars. There are instructions here, but I don't understand them. Could you explain in more detail how I can extract the mass locked in stars for parametric and non-parametric SFHs after the SED fitted?

  3. SFR: Is it possible to extract the SFR averaged over the past 100 years from the SED model or have the SFR as a free parameter in the fit?

  4. Attenuation: The Kriek & Conroy (2013) attenuation curve is determined by dust_index. However, the fsps manual says that dust_index is only used when dust_type = 0. If I would like to compute an EBV value from this model how would I have to proceed?

  5. Confidence intervals for best SEDs? The figure here shows contours of fits. How does one extract those from the final result file?

  6. Does the writer module have the option to overwrite hdf5 files?

  7. I want to change the labels of the corner plots. Is it possible to change the labels in reader.subcorner routine?

@bd-j
Copy link
Owner

bd-j commented Mar 5, 2020

Hi @steveschulze,

  1. No, sorry. You could create a branch and change https://github.com/bd-j/prospector/blob/master/prospect/sources/constants.py#L4 in your own repo and reinstall. For most cases it should make very very little difference.

  2. can you give more info about what you find confusing? The model and the obs object in the instructions at the link come from prospect.io.read_results.results_from; the sps object comes from prospect.io.read_results.get_sps.

  3. It is only possible to have SFR be a free parameter through some more or less complicated parameter transformations. In general we find it easier to work out the recent SFR after the fact based on the model parameters, but each SFH parameterization has a different conversion between model parameters and a recent SFR. Sorry there are not yet convenience functions for this.

  4. That part of the FSPS manual is out of date. Note however that dust_index has a different meaning for dust_type=4 than dust_type=0. Getting precise E(B-V) values for arbitrary dust curves is not trivial; how badly do you need them? One way would be to generate a spectrum as normal and also generate a spectrum with dust2=dust1=0 and compute the difference in the B-V colors of these spectra.

  5. Construct an SED for each of some number of fair draws from the posterior (i.e., using the weights if it's a dynesty run to get fair draws bi something like np.random.choice) and then compute the percentiles of the flux at each wavelength. The more draws you make, the more accurate your percentiles will be.

  6. By default a timestamp is added to each output file. But you could pass a filename without a timestamp to the writer method and I think that would overwrite any existing file of the same name

  7. No. But you could make your own version of the method (it is quite short) that takes a dictionary to remap the parameter names to something prettier.

@steveschulze
Copy link
Author

Hi Ben,

Thanks for your detailed answer.

  1. I read your post again. It's clearer now. I should probably just jump into the water.

  2. Do you know papers that show how to extract the SFRs for different SFHs? In your opinion, if I would use the attenuation-corrected UV luminosity and the standard Kennicutt relation, would this give a good estimate of the SFR averaged over the past 100 years?

  3. Having a precise value would not be a top priority. It would be just useful to have an estimate of how much a galaxy is reddened.

Cheers,

Steve

@bd-j
Copy link
Owner

bd-j commented Mar 27, 2020

Hi Steve, sorry for the delay but (gestures broadly at everything)

  1. For the nonparametric sfhs you can use prospect.models.transforms.logsfr_ratios_to_sfrs and then just read off the SFR of the most recent bin. For tau and delay-tau you just compute it (tau is e^{tage/tau}, delay-tau would be tage/tau * e^{tage/tau}) but you need to scale by the total mass formed (basically the integral of the previous functions, which is an incomplete gamma function)

  2. you could get the function for A_lambda(lambda, A_V, delta) from Kriek & Conroy, then replace A_V with 1.086 * dust2 and delta with dust_index, then compute the values at 4500AA (~B band) and 5500AA (~V band). Or if delta ~ 0, then E(B-V) = A_V/ 4.05

@astro-vic
Copy link

Hi Ben,
when you say "you need to scale by the total mass formed", how should I go about doing that? I fitted a galaxy that has a log(mass) = 10.34 and I found the SFR to be 0.02159 using Tage = 2.83 and tau = 0.51.

@bd-j
Copy link
Owner

bd-j commented Apr 2, 2020

you could do something like

import numpy as np
from scipy.special import gamma, gammainc

tage, tau, mass = 2.83, 0.51, 10**10.34
# for delay tau this function gives the (unnormalized) SFR 
# for any t, tau combo in M_sun/Gyr
sfr = lambda t,tau: return (t/tau) * np.exp(-t/tau)
# now we numerically integrate this SFH from 0 to tage to get the mass formed
times = np.linspace(0, tage, 1000)
A = np.trapz(sfr(times, tau), times)
# But this could also be done using an incomplete gamma function (integral of xe^{-x})
A = tau * gamma(2) * gammainc(2, tage/tau)
# and now we renormalize the formed mass to the actual mass value 
# to get the the SFR in M_sun per Gyr 
psi = mass * sfr(tage, tau) / A
# if we want SFR in Msun/year
psi /= 1e9

or, in one line

psi = mass * (tage/tau**2) * np.exp(-tage/tau) / (gamma(2) * gammainc(2, tage/tau)) * 1e-9

@steveschulze
Copy link
Author

Hi Ben,

I have a few more questions.

  1. In the interactive demo, you use the high-resolution templates to fit a broadband SED. Is it possible to load either low-res templates or rebin them to increase the speed?

  2. The fit in the interactive demo shows PAH features between 5 and 10 µm. Is it possible to load templates with only a stellar component?

  3. The IMF in the interactive demo was set to Chabrier (imf_type = 2). According to the FSPS manual, the Chabrier IMF is '1'. Could you confirm that Chabrier is imf_type = 2?

  4. In the discussion about how to compute the SFR, you mentioned that there is an exponentially declining SFH. The FSPS manual doesn't mention an exponentially declining SFH (http://dfm.io/python-fsps/current/stellarpop_api/?highlight=sfh).

Thanks a lot. :)

Steve

@bd-j
Copy link
Owner

bd-j commented Apr 8, 2020

  1. In principle yes, but the FSPS fortran source code has fixed wavelength gridding; the change has to be made to FSPS stellar libraries and source code.

  2. set model_params["add_dust_emission"]["init"] = False

  3. No, imf_type=2 is Kroupa; the interactive demo is wrong (thanks for finding this)

  4. A "tau model" (sfh=1) is exponentially declining.

@steveschulze
Copy link
Author

Thanks a lot Ben.

@steveschulze
Copy link
Author

Going back to the templates. Are there low-res templates is FSPS?

@bd-j
Copy link
Owner

bd-j commented Apr 8, 2020

you could use the older BaSeL theoretical stellar library; this is something like 1000 wavelength points instead of 7000 for MILES. See the FSPS documentation for how to do this (once you recompile FSPS you have to reinstall python-FSPS)

@steveschulze
Copy link
Author

Thanks.

@gkhullar
Copy link

gkhullar commented Apr 9, 2020

I would like to piggyback on Steve's question regarding SFR calculation (I am also happy to move this to a new issue/question). I am a little confused in the usage of tage vs sf_start for the tau/delayed-tau models here.

a) While fitting SFHs to a galaxy spectrum to calculate galaxy ages and/or star formation start times, should tage be a free parameter, or sf_start?

e.g. For calculations of the SFR in the time-step here (#166 (comment)):

times = np.linspace(0, tage, 1000)

Is this indicating that sf_start=0 (as is the case by default)? Should sf_start always be 0 while using prospector's parametric models?

b) If tage is fixed to the age of the universe at the epoch of observation and sf_start is a free parameter, wouldn't tage - sf_start technically be the age of the galaxy at the epoch of observation?

@bd-j
Copy link
Owner

bd-j commented Apr 9, 2020

It's (almost?) never necessary to set sf_start != 0- it's simple translation of the time axis. But yes, your logic in b) would hold.

@steveschulze
Copy link
Author

Is it possible to output the total stellar-mass formed and the mass locked into the living stars, simultaneously? Do any of the mass estimate take into account how the mass returned to the ISM?

@bd-j
Copy link
Owner

bd-j commented Apr 17, 2020

The formed stellar mass is an input parameter.

The ratio of the mass locked into stars (and remnants) to the formed mass is returned by SedModel.mean_model. It's sometimes referenced as mfrac in code or demos. I'm not sure what you mean by 'take into account', but the difference between the formed mass and the surviving stellar mass is presumably returned to the ISM; FSPS does not compute self-consistent chemical evolution.

@steveschulze
Copy link
Author

Thanks. :) For instance, BEAGLE includes a correction term for material returned to the ISM, in addition to computing the total mass formed and mass locked in living stars. These are only tiny details but important for comparing results of different SED fitting codes.

@idoirani
Copy link

Hi, I'd like to join in on the discussion here :)
In the fsps documentation I see that the number 4 sfh "delayed tau" model used in the parametrized sfh prospector model has 6 parameters. It has a constant component and an additional starburst component with an arbitrary fraction of the total mass formed. I'm adding the relevant paragraph bellow describing the parameters from the documentation. My question is are these (besides tau) set a fixed parameters anywhere? Shouldn't these parts of the model be taken into account when extracting the SFR, if they are non 0?

Thanks!
Ido

tau – (default: 1.0) Defines e-folding time for the SFH, in Gyr. Only used if sfh=1 or sfh=4. The range is 0.1<τ<1020.1<τ<102.
const – (default: 0.0) Defines the constant component of the SFH. This quantity is defined as the fraction of mass formed in a constant mode of SF; the range is therefore 0≤C≤10≤C≤1. Only used if sfh=1 or sfh=4.
sf_start – (default: 0.0) Start time of the SFH, in Gyr. Only used if sfh=1 or sfh=4 or sfh=5.
sf_trunc – (default: 0.0) Truncation time of the SFH, in Gyr. If set to 0.0, there is no trunction. Only used if sfh=1 or sfh=4 or sfh=5.
fburst – (default: 0.0) Defines the fraction of mass formed in an instantaneous burst of star formation. Only used if sfh=1 or sfh=4.
tburst – (default: 11.0) Defines the age of the Universe when the burst occurs. If tburst > tage then there is no burst. Only used if sfh=1 or sfh=4.

@jrleja
Copy link
Collaborator

jrleja commented May 14, 2020

Parameters which aren't specified in Prospector will be set automatically to their default python-fsps / FSPS values (I believe in that order, though I'm not sure of any meaningful conflicts). For the SFH=4 option, most of the 6 components have defaults which mean they do not affect the SFH. If the user does not modify sf_start, sf_trunc, fburst, and tburst, then SFH=4 will operate like a standard delayed-tau model.

If you do modify any of these parameters, you'll have to adjust the calculation of SFR accordingly.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

6 participants