Commit ba110b4a authored by Médéric Boquien's avatar Médéric Boquien

Implement the possibility of carrying out the parameter estimation on log values

parent 24cd822a
......@@ -2,6 +2,8 @@
## Unreleased
### Added
- The evaluation of the parameters is always done linearly. This can be a problem when estimating the SFR or the stellar mass for instance as it is usual to estimate their log rather. Because the log is non-linear, the likelihood-weigthed mean of the log is not the log of the likelihood-weighted mean. Therefore the estimation of the log of these parameters has to be done during the analysis step. This is now possible. The variables to be analysed in log just need to be indicated with the suffix "_log", for instance "stellar.m_star_log". (Médéric Boquien, idea suggested by Samir Salim)
### Changed
### Fixed
- Running the scripts in parallel trigger a deadlock on OS X with python 3.5. A workaround has been implemented. (Médéric Boquien)
......
......@@ -126,6 +126,8 @@ class PdfAnalysis(AnalysisModule):
# Initalise variables from input arguments.
analysed_variables = config["analysed_variables"]
analysed_variables_nolog = [variable.rstrip('_log') for variable in
analysed_variables]
n_variables = len(analysed_variables)
save = {key: config["save_{}".format(key)].lower() == "true"
for key in ["best_sed", "chi2", "pdf"]}
......@@ -182,7 +184,7 @@ class PdfAnalysis(AnalysisModule):
n_params * n_variables),
(n_params, n_variables))
initargs = (params, filters, analysed_variables, model_redshifts,
initargs = (params, filters, analysed_variables_nolog, model_redshifts,
model_fluxes, model_variables, time.time(),
mp.Value('i', 0))
if cores == 1: # Do not create a new process
......
......@@ -265,14 +265,23 @@ def analysis(idx, obs):
# We compute the weighted average and standard deviation using the
# likelihood as weight.
for i, variable in enumerate(gbl_analysed_variables):
if variable.endswith('_log'):
variable = variable.strip('_log')
_ = np.log10
maxstd = lambda mean, std: max(0.02, std)
else:
_ = lambda x: x
maxstd = lambda mean, std: max(0.05 * mean, std)
if variable in sed.mass_proportional_info:
mean, std = weighted_param(gbl_model_variables[wz, i][wlikely]
* scaling[wlikely], likelihood)
mean, std = weighted_param(_(gbl_model_variables[wz, i][wlikely]
* scaling[wlikely]), likelihood)
else:
mean, std = weighted_param(gbl_model_variables[wz, i][wlikely],
mean, std = weighted_param(_(gbl_model_variables[wz, i][wlikely]),
likelihood)
gbl_analysed_averages[idx, i] = mean
gbl_analysed_std[idx, i] = max(0.05 * mean, std)
gbl_analysed_std[idx, i] = maxstd(mean, std)
gbl_best_fluxes[idx, :] = gbl_model_fluxes[best_index, :] \
* scaling[best_index_z]
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment