Posterior.h

Go to the documentation of this file.

/********************************************************************
 *  Copyright (C) 2014 by Federico Marulli and Alfonso Veropalumbo  *
 *  federico.marulli3@unibo.it                                      *
 *                                                                  *
 *  This program is free software; you can redistribute it and/or   *
 *  modify it under the terms of the GNU General Public License as  *
 *  published by the Free Software Foundation; either version 2 of  *
 *  the License, or (at your option) any later version.             *
 *                                                                  *
 *  This program is distributed in the hope that it will be useful, *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of  *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the   *
 *  GNU General Public License for more details.                    *
 *                                                                  *
 *  You should have received a copy of the GNU General Public       *
 *  License along with this program; if not, write to the Free      *
 *  Software Foundation, Inc.,                                      *
 *  59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.       *
 ********************************************************************/
 
#ifndef __POSTERIOR__
#define __POSTERIOR__
 
#include "RandomNumbers.h"
#include "Prior.h"
#include "Likelihood.h"
 
 
// ===================================================================================================
 
 
namespace cbl {
 
  namespace statistics {
 
    class Posterior : public Likelihood {
 
    protected:
 
      std::shared_ptr<cbl::statistics::Prior> m_prior;
      
      std::shared_ptr<statistics::Model> m_response_func = NULL;
 
      std::vector<double> m_log_likelihood;
 
      std::vector<double> m_acceptance;
 
      int m_seed;
 
      std::shared_ptr<cbl::random::UniformRandomNumbers_Int> m_seed_generator;
 
      std::vector<double> m_log_posterior;
 
      int m_Nparameters;
 
      std::vector<std::vector<double>> m_parameters;
 
      std::vector<double> m_weight;
 
      void m_set_seed (const int seed);
 
      int m_generate_seed () { return m_seed_generator->operator()(); }
      
      void m_isSet_response () {if (m_response_func == NULL) ErrorCBL("the response function for the SSC is not set!", "m_isSet_response", "Posterior.h");};
 
    public:
 
 
      Posterior () = default;
 
      Posterior (const std::vector<std::shared_ptr<PriorDistribution>> prior_distributions, const Likelihood &likelihood, const int seed=5341);
 
      Posterior (const std::vector<std::shared_ptr<PriorDistribution>> prior_distributions, 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 int seed=5341);
 
      Posterior (const std::string file_name, const std::string path_file, const std::vector<int> usecols,  const std::vector<std::string> parameter_names, const int skip_header);
 
      ~Posterior () = default;
 
    
    
      std::shared_ptr<ModelParameters> parameters () const { return m_model_parameters; }
 
      double operator () (std::vector<double> &pp) const;
 
      std::vector<double> weight (const int start=0, const int thin=1)  const;
 
      double log (std::vector<double> &pp) const;
 
      double chi2 (const std::vector<double> parameter={}) const;
 
      void set_model (std::shared_ptr<Model> model=NULL, const std::shared_ptr<ModelParameters> model_parameters=NULL);
 
      void set (const std::vector<std::shared_ptr<PriorDistribution>> prior_distributions, 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 int seed);
 
      int get_Nparameters () { return m_Nparameters; }
 
      int m_get_seed () { return m_seed; }
 
      bool get_m_use_grid() { return m_use_grid; }
 
      std::vector<std::vector<double>> get_parameters () { return m_parameters; }
 
      std::vector<double> get_log_posterior () { return m_log_posterior; }
 
      std::shared_ptr<void> get_m_likelihood_inputs() { return m_likelihood_inputs; }
 
      std::shared_ptr<cbl::statistics::Prior> get_m_prior() { return m_prior; }
 
      LogLikelihood_function get_m_log_likelihood_function () { return m_log_likelihood_function; };
 
      Likelihood_function get_m_likelihood_function () { return m_likelihood_function; };
 
      Likelihood_function get_m_likelihood_function_grid () { return m_likelihood_function_grid; };
 
      LogLikelihood_function get_m_log_likelihood_function_grid () { return m_log_likelihood_function_grid; };
 
      void set_response_function (std::shared_ptr<statistics::Model> response) {m_response_func = response;};
 
      std::shared_ptr<statistics::Model> get_response_function () {m_isSet_response(); return m_response_func;};
 
 
      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)
      { (void)start; (void)parameter_limits; (void)max_iter; (void)tol; (void)epsilon; ErrorCBL("the method is used without parameter_ranges!", "maximize", "Posterior.h"); }
 
      void maximize (const std::vector<double> start, const unsigned int max_iter=10000, const double tol=1.e-6, const double epsilon=1.e-4);
 
      void initialize_chains (const int chain_size, const int n_walkers);
 
      void initialize_chains (const int chain_size, const int n_walkers, const double radius, const std::vector<double> start, const unsigned int max_iter=10000, const double tol=1.e-6, const double epsilon=1.e-3);
 
      void initialize_chains (const int chain_size, const int n_walkers, std::vector<double> &value, const double radius);
 
      void initialize_chains (const int chain_size, const std::vector<std::vector<double>> chain_value);
 
      void initialize_chains (const int chain_size, const int n_walkers, const std::string input_dir, const std::string input_file);
 
      void sample_stretch_move (const double aa=2, const bool parallel=true, const std::string outputFile=par::defaultString, const int start=0, const int thin=1, const int nbins=50);
 
      void importance_sampling (const std::string input_dir, const std::string input_file, const std::vector<size_t> column={}, const int header_lines_to_skip=1, const bool is_FITS_format=false, const bool apply_to_likelihood=false, const int n_walkers=100);
 
      void read_chain (const std::string input_dir, const std::string input_file, const int n_walkers, const std::vector<size_t> column={}, const int header_lines_to_skip=1, const bool is_FITS_format=false);
    
      void read_chain_ascii (const std::string input_dir, const std::string input_file, const int n_walkers, const std::vector<size_t> column={}, const int header_lines_to_skip=1);
 
      void read_chain_fits (const std::string input_dir, const std::string input_file, const int n_walkers, const std::vector<size_t> column);
 
      void write_chain (const std::string output_dir, const std::string output_file, const int start=0, const int thin=1, const bool is_FITS_format=false, const int prec=5, const int ww=14);
    
      void write_chain_ascii (const std::string output_dir, const std::string output_file, const int start=0, const int thin=1, const int prec=5, const int ww=14);
 
      void write_chain_fits (const std::string output_dir, const std::string output_file, const int start=0, const int thin=1);
 
    
      void write_maximization_results (const std::string output_dir, const std::string root_file);
 
      void show_results (const int start, const int thin, const int nbins=50, const bool show_mode=false, const int ns=-1, const int nb=-1);
 
      void write_results (const std::string output_dir, const std::string root_file, const int start=0, const int thin=1, const int nbins=50, const bool fits=false, const bool compute_mode=false, const int ns=-1, const int nb=-1);
 
      void write_model_from_chain (const std::string output_dir, const std::string output_file, const std::vector<double> xx={}, const std::vector<double> yy={}, const int start=0, const int thin=1);
    };
  }
}
 
#endif