CRejectionSamplingCapable.h

Go to the documentation of this file.

/* +------------------------------------------------------------------------+
   |                     Mobile Robot Programming Toolkit (MRPT)            |
   |                          http://www.mrpt.org/                          |
   |                                                                        |
   | Copyright (c) 2005-2017, Individual contributors, see AUTHORS file     |
   | See: http://www.mrpt.org/Authors - All rights reserved.                |
   | Released under BSD License. See details in http://www.mrpt.org/License |
   +------------------------------------------------------------------------+ */
#ifndef CRejectionSamplingCapable_H
#define CRejectionSamplingCapable_H

#include <mrpt/utils/utils_defs.h>
#include <mrpt/bayes/CProbabilityParticle.h>
#include <mrpt/random.h>

namespace mrpt
{
/// \ingroup mrpt_bayes_grp
namespace bayes
{
/** A base class for implementing rejection sampling in a generic state space.
 *   See the main method CRejectionSamplingCapable::rejectionSampling
 *  To use this class, create your own class as a child of this one and
 * implement the desired
 *   virtual methods, and add any required internal data.
 * \ingroup mrpt_bayes_grp
 */
template <class TStateSpace>
class CRejectionSamplingCapable
{
   public:
        typedef CProbabilityParticle<TStateSpace> TParticle;

        /** Virtual destructor
          */
        virtual ~CRejectionSamplingCapable() {}
        /** Generates a set of N independent samples via rejection sampling.
          * \param desiredSamples The number of desired samples to generate
          * \param outSamples The output samples.
          * \param timeoutTrials The maximum number of rejection trials for each
         * generated sample (i.e. the maximum number of iterations). This can be
         * used to set a limit to the time complexity of the algorithm for difficult
         * probability densities.
          *  All will have equal importance weights (a property of rejection
         * sampling), although those samples
          *   generated at timeout will have a different importance weights.
          */
        void rejectionSampling(
                size_t desiredSamples, std::vector<TParticle>& outSamples,
                size_t timeoutTrials = 1000)
        {
                MRPT_START

                TStateSpace x;
                typename std::vector<TParticle>::iterator it;

                // Set output size:
                if (outSamples.size() != desiredSamples)
                {
                        // Free old memory:
                        outSamples.clear();

                        // Reserve new memory:
                        outSamples.resize(desiredSamples);
                        for (it = outSamples.begin(); it != outSamples.end(); ++it)
                                it->d.reset(new TStateSpace);
                }

                // Rejection sampling loop:
                double acceptanceProb;
                for (it = outSamples.begin(); it != outSamples.end(); ++it)
                {
                        size_t timeoutCount = 0;
                        double bestLik = -1e250;
                        TStateSpace bestVal;
                        do
                        {
                                RS_drawFromProposal(*it->d);
                                acceptanceProb = RS_observationLikelihood(*it->d);
                                ASSERT_(acceptanceProb >= 0 && acceptanceProb <= 1);
                                if (acceptanceProb > bestLik)
                                {
                                        bestLik = acceptanceProb;
                                        bestVal = *it->d;
                                }
                        } while (acceptanceProb < mrpt::random::getRandomGenerator().drawUniform(
                                                                                  0.0, 0.999) &&
                                         (++timeoutCount) < timeoutTrials);

                        // Save weights:
                        if (timeoutCount >= timeoutTrials)
                        {
                                it->log_w = log(bestLik);
                                *it->d = bestVal;
                        }
                        else
                        {
                                it->log_w = 0;  // log(1.0);
                        }
                }  // end for it

                MRPT_END
        }

   protected:
        /** Generates one sample, drawing from some proposal distribution.
          */
        virtual void RS_drawFromProposal(TStateSpace& outSample) = 0;

        /** Returns the NORMALIZED observation likelihood (linear, not
         * exponential!!!) at a given point of the state space (values in the range
         * [0,1]).
          */
        virtual double RS_observationLikelihood(const TStateSpace& x) = 0;

};  // End of class def.

}  // End of namespace
}  // End of namespace

#endif