Line data    Source code

       1             : /**
       2             :    Bojan Nikolic <bojan@bnikolic.co.uk> 
       3             :    Initial version 2008
       4             : 
       5             :    This file is part of BNMin1 and is licensed under GNU General
       6             :    Public License version 2
       7             : 
       8             :    \file metropolis.hxx
       9             :    Renamed to metropolis.h 2023
      10             : 
      11             :    Metropolis sampling of the likelihood 
      12             : */
      13             : #ifndef _BNMIN1_METROPOLIS_HXX__
      14             : #define _BNMIN1_METROPOLIS_HXX__
      15             : 
      16             : #include <stdexcept>
      17             : #include <list>
      18             : #include <vector>
      19             : 
      20             : #include <memory> // was <boost/scoped_ptr.hpp>
      21             : #include <functional> // was <boost/function.hpp>
      22             : 
      23             : #include "minim.h"
      24             : #include "minimmodel.h"
      25             : #include "mcpoint.h"
      26             : 
      27             : 
      28             : 
      29             : namespace Minim {
      30             : 
      31             :   class MLikelihood;
      32             :   class MetroPropose;
      33             : 
      34             :   /** \brief Markov Chain Monte Carlo with Metropolis sampling
      35             :    */
      36             :   class MetropolisMCMC:
      37             :     public ModelDesc 
      38             :   {
      39             : 
      40             :     MLikelihood & ml;
      41             :     
      42             :     std::unique_ptr<MetroPropose> prop;
      43             : 
      44             :   public:
      45             : 
      46             :     class Error :
      47             :       public std::runtime_error
      48             :     {
      49             :     public:
      50           0 :       Error(const std::string &msg):
      51           0 :         std::runtime_error(msg)
      52             :       {
      53           0 :       }
      54             :     };
      55             :     
      56             :     /// Options for the Metropolis MCMC sampling algorithm
      57             :     enum Options {
      58             :       /// No special options, default
      59             :       NoOp=0,
      60             :       /** Normally proposal points are generated by perturbing the
      61             :       whole parameter vector simultanousely. If this option is set,
      62             :       each parameter is perturbed in turn. This reduces the impact of
      63             :       sigmas which are too large by allowing other parameters to
      64             :       evolve. However, if sigmas are too low and there are many
      65             :       parameters, time to converge increases.
      66             :       */
      67             :       Sequence=1};
      68             : 
      69             :     // ---------- Public data -----------------------------
      70             : 
      71             :     /** \brief A function to record at each accepted distribution
      72             :         point
      73             :         
      74             :         If this pointer is set, this function will be called at each
      75             :         point in the distribution and values stored in parameter x
      76             :         will be saved in the field fval.
      77             :     */
      78             :     std::function< void (std::vector<double> &x) >  f;    
      79             : 
      80             :     // ---------- Construction / Destruction --------------
      81             : 
      82             :     /**
      83             :        \brief Initialise the class
      84             : 
      85             :        \param sigmas The standard deviations of the parameter
      86             :        perturbation used in the Metropolis algorithm.
      87             :        
      88             :        \param seed The seed for the random number generator
      89             : 
      90             :        \param opt Options
      91             :      */
      92             :     MetropolisMCMC(MLikelihood & ml,
      93             :                    const std::vector<double> & sigmas,
      94             :                    unsigned seed=0,
      95             :                    Options opt=NoOp
      96             :                    );
      97             :     
      98             :     virtual ~MetropolisMCMC();
      99             : 
     100             :     // ---------- Public interface --------------------------
     101             : 
     102             :     /**
     103             :        \brief Create a sample chain.
     104             : 
     105             :        \param npropose number of proposal samples to run
     106             : 
     107             :      */
     108             :     std::list<MCPoint> *
     109             :     sample(size_t npropose);
     110             :     
     111             : 
     112             : 
     113             :   };
     114             : 
     115             : 
     116             : }
     117             : #endif