45 vector<double> pp(npoints), log_ll(npoints);
46 double binSize = (parameter_limits[0][1]-parameter_limits[0][0])/(npoints-1);
47 for (
int i=0; i<npoints; i++) {
48 pp[i] = parameter_limits[0][0]+binSize*i;
49 vector<double> par = {pp[i]};
50 log_ll[i] = this->log(par);
53 shared_ptr<statistics::STR_likelihood_inputs> likelihood_inputs = static_pointer_cast<statistics::STR_likelihood_inputs>(m_likelihood_inputs);
54 likelihood_inputs->interp_function1D = make_shared<glob::FuncGrid>(
glob::FuncGrid(pp, log_ll,
"Spline"));
57 m_likelihood_function_grid = [&] (vector<double> &par,
const shared_ptr<void> input) {
return exp(m_log_likelihood_function_grid(par, input));};
60 ofstream fout(output_file);
checkIO(fout, output_file);
61 for (
int i=0; i<npoints; i++)
62 fout << pp[i] <<
" " << log_ll[i] << endl;
74 vector<double> pp1(npoints), pp2(npoints);
75 vector<vector<double>> log_ll(npoints, vector<double>(npoints, 0));
76 double binSize1 = (parameter_limits[0][1]-parameter_limits[0][0])/(npoints-1);
77 double binSize2 = (parameter_limits[1][1]-parameter_limits[1][0])/(npoints-1);
78 for (
int i=0; i<npoints; i++) {
79 pp1[i] = parameter_limits[0][0]+binSize1*i;
80 for (
int j=0; j<npoints; j++) {
81 pp2[j] = parameter_limits[1][0]+binSize2*j;
82 vector<double> par = {pp1[i], pp2[j]};
83 log_ll[i][j] = this->log(par);
87 shared_ptr<statistics::STR_likelihood_inputs> likelihood_inputs = static_pointer_cast<statistics::STR_likelihood_inputs>(m_likelihood_inputs);
88 likelihood_inputs->interp_function2D = make_shared<glob::FuncGrid2D>(
glob::FuncGrid2D(pp1, pp2, log_ll,
"Spline"));
91 m_likelihood_function_grid = [&] (vector<double> &par,
const shared_ptr<void> input) {
return exp(m_log_likelihood_function_grid(par, input));};
94 ofstream fout(output_file);
checkIO(fout, output_file);
95 for (
int i=0; i<npoints; i++) {
96 for (
int j=0; j<npoints; j++)
97 fout << pp1[i] <<
" " << pp2[j] <<
" " << log_ll[i][j] << endl;
112 vector<double> pp, log_ll;
116 shared_ptr<statistics::STR_likelihood_inputs> likelihood_inputs = static_pointer_cast<statistics::STR_likelihood_inputs>(m_likelihood_inputs);
117 likelihood_inputs->interp_function1D = make_shared<glob::FuncGrid>(
glob::FuncGrid(pp, log_ll,
"Spline"));
120 m_likelihood_function_grid = [&] (vector<double> &par,
const shared_ptr<void> input) {
return exp(m_log_likelihood_function_grid(par, input));};
130 vector<double> pp1, pp2;
131 vector<vector<double>> log_ll;
135 shared_ptr<statistics::STR_likelihood_inputs> likelihood_inputs = static_pointer_cast<statistics::STR_likelihood_inputs>(m_likelihood_inputs);
136 likelihood_inputs->interp_function2D = make_shared<glob::FuncGrid2D>(
glob::FuncGrid2D(pp1, pp2, log_ll,
"Spline"));
139 m_likelihood_function_grid = [&] (vector<double> &par,
const shared_ptr<void> input) {
return exp(m_log_likelihood_function_grid(par, input));};
146 cbl::statistics::Likelihood::Likelihood (
const std::shared_ptr<data::Data> data,
const std::shared_ptr<Model> model,
const LikelihoodType likelihood_type,
const std::vector<size_t> x_index,
const int w_index,
const std::shared_ptr<ModelParameters> model_parameters,
const double prec,
const int Nres)
149 set_model(model, model_parameters);
150 set_function(likelihood_type, x_index, w_index, prec, Nres);
160 set_model(model, model_parameters);
161 set_log_function(log_likelihood_function);
170 return exp(this->log(pp));
179 return (m_use_grid) ? m_log_likelihood_function_grid(parameter, m_likelihood_inputs) : m_log_likelihood_function(parameter, m_likelihood_inputs);
197 switch (model->dimension()) {
199 m_model = make_shared<Model1D>(*static_pointer_cast<Model1D>(model));
202 m_model = make_shared<Model2D>(*static_pointer_cast<Model2D>(model));
205 ErrorCBL(
"the model dimension must be Dim::_1D_ or Dim::_2D_ !",
"set_model",
"Likelihood.cpp");
208 m_model_parameters = (model_parameters == NULL) ? make_shared<LikelihoodParameters>(
LikelihoodParameters(m_model->parameters()->nparameters(), m_model->parameters()->type(), m_model->parameters()->name())) : model_parameters;
209 m_model->set_parameters(m_model_parameters);
227 if (m_likelihood_type==statistics::LikelihoodType::_NotSet_)
228 ErrorCBL(
"the Likelihood function is not set!",
"set_grid",
"Likelihood.cpp");
230 unsigned int npar = m_model->parameters()->nparameters_free();
233 ErrorCBL(
"there is no parameter free to vary!",
"set_grid",
"Likelihood.cpp");
235 ErrorCBL(
"wrong size for the vector of starting parameters!",
"set_grid",
"Likelihood.cpp");
236 if (parameter_limits.size()!=npar)
237 ErrorCBL(
"wrong size for the vector of parameter limits!",
"set_grid",
"Likelihood.cpp");
239 coutCBL <<
"Computing tabulated likelihood!" << endl;
243 m_set_grid_likelihood_1D(file);
245 m_set_grid_likelihood_1D(npoints, parameter_limits, file);
249 m_set_grid_likelihood_2D(file);
251 m_set_grid_likelihood_2D(npoints, parameter_limits, file);
266 m_likelihood_type = likelihood_type;
268 if (m_data->dataType()==data::DataType::_1D_ || m_data->dataType()==data::DataType::_1D_extra_) {
270 switch (m_likelihood_type)
273 case (LikelihoodType::_Gaussian_Error_):
274 for (
int i=0; i<m_data->ndata(); i++) {
275 if (std::isinf(pow(1./m_data->error(i), 2)))
284 case (LikelihoodType::_Gaussian_Covariance_):
285 m_data->invert_covariance(prec, Nres);
289 case (LikelihoodType::_Poissonian_):
294 ErrorCBL(
"type of likelihood not recognized or not yet implemented!",
"set_function",
"Likelihood.cpp");
299 else if (m_data->dataType()==data::DataType::_2D_ || m_data->dataType()==data::DataType::_2D_extra_) {
300 switch (m_likelihood_type)
303 case (LikelihoodType::_Gaussian_Error_):
304 for (
int i=0; i<m_data->xsize(); i++) {
305 for (
int j=0; j<m_data->ysize(); j++) {
306 if (std::isinf(pow(1./m_data->error(i, j), 2)))
315 case (LikelihoodType::_Poissonian_):
320 ErrorCBL(
"type of likelihood not recognized or not yet implemented!",
"set_function",
"Likelihood.cpp");
325 ErrorCBL(
"data type not recognized or not yet implemented!",
"set_function",
"Likelihood.cpp");
327 m_likelihood_function = [&] (vector<double> &par,
const shared_ptr<void> input) {
return exp(m_log_likelihood_function(par, input)); };
336 m_likelihood_type = LikelihoodType::_UserDefined_;
337 m_likelihood_function = likelihood_function;
338 m_log_likelihood_function = [&] (vector<double> &par,
const shared_ptr<void> input) {
return std::log(m_likelihood_function(par, input)); };
347 m_likelihood_type = LikelihoodType::_UserDefined_;
348 m_log_likelihood_function = log_likelihood_function;
349 m_likelihood_function = [&] (vector<double> &par,
const shared_ptr<void> input) {
return exp(m_log_likelihood_function(par, input)); };
356 void cbl::statistics::Likelihood::maximize (
const std::vector<double> start,
const std::vector<std::vector<double>> parameter_limits,
const unsigned int max_iter,
const double tol,
const double epsilon)
358 if (m_likelihood_type==statistics::LikelihoodType::_NotSet_)
359 ErrorCBL(
"the Likelihood function is not set!",
"maximize",
"Likelihood.cpp");
361 vector<double> starting_par;
362 vector<vector<double>> limits_par;
364 unsigned int npar_free = m_model->parameters()->nparameters_free();
365 unsigned int npar = m_model->parameters()->nparameters();
367 if (start.size()==npar_free && parameter_limits.size()==npar_free)
369 starting_par = start;
370 limits_par = parameter_limits;
372 else if (start.size()==npar && parameter_limits.size()==npar)
374 for (
size_t i=0; i<npar_free; i++) {
375 starting_par.push_back(start[m_model->parameters()->free_parameter()[i]]);
376 limits_par.push_back(parameter_limits[m_model->parameters()->free_parameter()[i]]);
380 ErrorCBL(
"check your inputs!",
"maximize",
"Likelihood.cpp");
382 m_likelihood_inputs = make_shared<STR_likelihood_inputs>(
STR_likelihood_inputs(m_data, m_model, m_x_index, m_w_index));
384 function<double(vector<double> &)> ll = [
this] (vector<double> &pp)
386 return -m_log_likelihood_function(pp, m_likelihood_inputs);
389 coutCBL <<
"Maximizing the likelihood..." << endl;
391 coutCBL <<
"Done!" << endl << endl;
393 m_model_parameters->set_bestfit_values(result);
394 m_model_parameters->write_bestfit_info();
395 coutCBL <<
"log(Likelihood) = " << this->log(result) << endl << endl;
404 coutCBL <<
"Writing results of Likelihood minimization on " << dir_output+file << endl;
405 vector<double> bestFitValue = m_model_parameters->bestfit_value();
407 double likelihoodValue = this->log(bestFitValue);
409 string mkdir =
"mkdir -p "+dir_output;
410 if (system(mkdir.c_str())) {}
412 ofstream fout(dir_output+file);
414 fout <<
"#Parameters information" << endl;
415 fout <<
"nParameters = " << bestFitValue.size() << endl;
417 for (
size_t i=0; i<bestFitValue.size(); i++) {
418 fout <<
"par" << i+1 <<
"_name = " << m_model_parameters->name(i) << endl;
419 fout <<
"par" << i+1 <<
"_status = " << m_model_parameters->status(i) << endl;
420 fout <<
"par" << i+1 <<
"_bestfit_value = " << bestFitValue[i] << endl;
423 fout <<
"#Likelihood information" << endl;
424 fout <<
"likelihoodType = " << name << endl;
425 fout <<
"logLikelihoodValue = " << likelihoodValue << endl;
427 fout.clear(); fout.close();
428 coutCBL <<
"I wrote the file " << dir_output+file << endl;
437 switch (m_model->dimension()) {
441 vector<double> xvec = xx;
445 m_model->write(output_dir, output_file, xvec, parameters);
451 vector<double> xvec = xx, yvec = yy;
457 m_model->write(output_dir, output_file, xvec, yvec, parameters);
462 ErrorCBL(
"the input dimension must be Dim::_1D_ or Dim::_2D_!",
"write_model",
"Likelihood.cpp");
#define coutCBL
CBL print message.
The class LikelihoodParameters.
The class LikelihoodParameters.
void write_model(const std::string output_dir, const std::string output_file, const std::vector< double > parameters, const std::vector< double > xx={}, const std::vector< double > yy={})
write the model at xx at bestfit
void write_results(const std::string dir_output, const std::string file)
write best-fit results on a file
void set_model(std::shared_ptr< Model > model=NULL, const std::shared_ptr< ModelParameters > model_parameters=NULL)
set the model for the likelihood analysis
void maximize(const std::vector< double > start, const std::vector< std::vector< double >> parameter_limits, const unsigned int max_iter=10000, const double tol=1.e-6, const double epsilon=1.e-3)
function that maximizes the likelihood, finds the best-fit parameters and stores them in the model
void m_set_grid_likelihood_1D(const int npoints, const std::vector< std::vector< double >> parameter_limits, const std::string output_file)
set the likelihood grid and write the grid on a file
double log(std::vector< double > ¶meter) const
evaluate the logarithm of the likelihood for the input parameters
void set_grid(const int npoints, const std::vector< std::vector< double >> parameter_limits, const std::string file, const bool read=false)
set the likelihood function as a grid, to speed up computation: this only works for one or two free p...
double operator()(std::vector< double > &pp) const
evaluate the likelihood
void unset_grid()
set the likelihood function with internal values of LikelihoodType
void set_log_function(const LogLikelihood_function loglikelihood_function)
set the natural logarithm of the likelihood function
void m_set_grid_likelihood_2D(const int npoints, const std::vector< std::vector< double >> parameter_limits, const std::string output_file)
set the likelihood grid and write the grid on a file
void set_function(const LikelihoodType likelihood_type, const std::vector< size_t > x_index={0, 2}, const int w_index=-1, const double prec=1.e-10, const int Nres=-1)
set the likelihood type using the LikelihoodType object
void set_data(std::shared_ptr< data::Data > data)
set the data for the likelihood analysis
Likelihood()
default constructor
static const char fINT[]
conversion symbol for: int -> std::string
static const std::string defaultString
default std::string value
double LogLikelihood_Gaussian_1D_covariance(std::vector< double > &likelihood_parameter, const std::shared_ptr< void > input)
function to compute the gaussian loglikelihood
double LogLikelihood_Gaussian_1D_error(std::vector< double > &likelihood_parameter, const std::shared_ptr< void > input)
function to compute the Gaussian log-likelihood
double LogLikelihood_Poissonian_1D_(std::vector< double > &likelihood_parameter, const std::shared_ptr< void > input)
function to compute the poissonian loglikelihood
LikelihoodType
the type of likelihood function
std::vector< std::string > LikelihoodTypeNames()
return a vector containing the LikelihoodType names
double LogLikelihood_1D_interpolated(std::vector< double > &likelihood_parameter, const std::shared_ptr< void > input)
function to compute the loglikelihood on a grid
std::function< double(std::vector< double > &, const std::shared_ptr< void >)> Likelihood_function
definition of a function for computation of the Likelihood
double LogLikelihood_Poissonian_2D_(std::vector< double > &likelihood_parameter, const std::shared_ptr< void > input)
function to compute the poissonian loglikelihood
double LogLikelihood_2D_interpolated(std::vector< double > &likelihood_parameter, const std::shared_ptr< void > input)
function to compute the loglikelihood on a grid
double LogLikelihood_Gaussian_2D_error(std::vector< double > &likelihood_parameter, const std::shared_ptr< void > input)
function to compute the gaussian loglikelihood model with one parameter
std::vector< double > GSL_minimize_nD(FunctionDoubleVector func, const std::vector< double > start, const std::vector< std::vector< double >> ranges, const unsigned int max_iter=1000, const double tol=1.e-6, const double epsilon=0.1)
minimize the provided function using GSL procedure
The global namespace of the CosmoBolognaLib
std::string conv(const T val, const char *fact)
convert a number to a std::string
void read_matrix(const std::string file_matrix, std::vector< double > &xx, std::vector< double > &yy, std::vector< std::vector< double >> &matrix, const std::vector< int > col={})
read a matrix from file
void checkIO(const std::ifstream &fin, const std::string file="NULL")
check if an input file can be opened
int ErrorCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_)
throw an exception: it is used for handling exceptions inside the CosmoBolognaLib
void read_vector(const std::string file_vector, std::vector< double > &xx, std::vector< double > &vector, const std::vector< int > col={})
read a vector from file