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.cxx
       9             :    Renamed to metropolis.cc 2023
      10             : 
      11             : */
      12             : 
      13             : #include <list>
      14             : #include <memory>
      15             : #include <cmath>
      16             : #include <sstream>
      17             : 
      18             : #include "metropolis.h"
      19             : #include "metro_propose.h"
      20             : 
      21             : 
      22             : namespace Minim {
      23             : 
      24           0 :   MetropolisMCMC::MetropolisMCMC(MLikelihood & ml,
      25             :                                  const std::vector<double> & sigmas,
      26             :                                  unsigned seed,
      27           0 :                                  Options opt):
      28             :     ModelDesc(ml),
      29           0 :     ml(ml),
      30           0 :     f(NULL)
      31             :   {
      32           0 :     if (opt & Sequence)
      33           0 :       prop.reset(new MetroProposeSeq(sigmas, seed));
      34             :     else
      35           0 :       prop.reset(new MetroPropose(sigmas, seed));
      36           0 :   }
      37             : 
      38           0 :   MetropolisMCMC::~MetropolisMCMC()
      39             :   {
      40           0 :   }
      41             : 
      42             :   std::list<MCPoint> *
      43           0 :   MetropolisMCMC::sample(size_t npropose)
      44             :   {
      45           0 :     if (prop->nPars() != NParam() )
      46             :     {
      47           0 :       std::ostringstream buf;
      48           0 :       buf << "Number of sigmas not consistent with number of pars set to fit: expected " << NParam( ) << ", got " << prop->nPars( );
      49           0 :       throw Error(buf.str( ));
      50           0 :     }
      51             : 
      52             :     std::unique_ptr< std::list<MCPoint> > res
      53           0 :       (new std::list<MCPoint> ) ;
      54             : 
      55             :     // Current point in chain
      56           0 :     std::vector<double> cpoint( prop->nPars());
      57           0 :     get(cpoint);
      58             : 
      59             :     // Current likelihood
      60           0 :     double clogl = - ml.lLikely();
      61             :     
      62           0 :     std::vector<double> propose( prop->nPars() );
      63             :     
      64           0 :     for(size_t sn =0 ; sn < npropose; ++sn)
      65             :     {
      66           0 :       propose=cpoint;
      67           0 :       prop->displace(propose);      
      68           0 :       put(propose);
      69             : 
      70           0 :       double proplogl = - ml.lLikely();
      71           0 :       double    ldiff = proplogl - clogl;
      72             : 
      73             :       // If new point more likely than old, definitly accept (ldiff
      74             :       // >0). If new point less likely accept with probability
      75             :       // exp(ldiff)
      76             : 
      77           0 :       if (ldiff > 0 or 
      78           0 :           prop->raccept() < exp(ldiff) )
      79             :       {
      80           0 :         MCPoint mcp(propose);
      81           0 :         mcp.ll=proplogl;
      82           0 :         if (f)
      83           0 :           f(mcp.fval);
      84           0 :         res->push_back(mcp);
      85           0 :         cpoint=propose;
      86           0 :         clogl =proplogl;      
      87           0 :       }
      88             :     }
      89             : 
      90           0 :     return res.release();
      91             : 
      92           0 :   }
      93             : 
      94             :   
      95             : 
      96             : 
      97             : }
      98             :