rv_samp.h

Go to the documentation of this file.

#ifndef RV_SAMP_H
#define RV_SAMP_H
 
#include <chrono>
 
#ifdef DROPPINGTHISINRPACKAGE
    #include <RcppEigen.h>
    // [[Rcpp::depends(RcppEigen)]]
#else
    #include <Eigen/Dense>
#endif
 
#include <random>
 
 
namespace pf {
 
namespace rvsamp{
 
 
 
class rvsamp_base
{
public:
 
    inline rvsamp_base() : 
        m_rng{static_cast<std::uint32_t>(std::chrono::high_resolution_clock::now().time_since_epoch().count())} 
    {}
 
protected:
 
    std::mt19937 m_rng;
 
};
 
 
 
template<typename float_t>
class UnivNormSampler : public rvsamp_base
{
    
public:
 
 
    UnivNormSampler();
 
 
    UnivNormSampler(float_t mu, float_t sigma);
 
 
    void setStdDev(float_t sigma);
    
    
    void setMean(float_t mu);
    
        
    float_t sample();    
    
 
private:
 
    std::normal_distribution<float_t> m_z_gen;
    
    float_t m_mu;
    
    float_t m_sigma;
 
};
 
 
template<typename float_t>
UnivNormSampler<float_t>::UnivNormSampler()
    : rvsamp_base()
    , m_z_gen(0.0, 1.0)
{
    setMean(0.0);
    setStdDev(1.0);
}
 
 
template<typename float_t>
UnivNormSampler<float_t>::UnivNormSampler(float_t mu, float_t sigma)
    : rvsamp_base()
    , m_z_gen(0.0, 1.0)
{
    setMean(mu); 
    setStdDev(sigma);
}
 
 
template<typename float_t>
void UnivNormSampler<float_t>::setMean(float_t mu)
{
    m_mu = mu;
}
 
 
template<typename float_t>
void UnivNormSampler<float_t>::setStdDev(float_t sigma)
{
    m_sigma = sigma;
}
 
 
template<typename float_t>
float_t UnivNormSampler<float_t>::sample()
{
    return m_mu + m_sigma * m_z_gen(m_rng);
}
 
 
 
template<typename float_t>
class UnivStudTSampler : public rvsamp_base
{
    
public:
 
 
    UnivStudTSampler() = delete;
 
 
    UnivStudTSampler(float_t dof);
 
 
    float_t sample();    
    
 
private:
 
    std::student_t_distribution<float_t> m_t_gen;
    
};
 
 
template<typename float_t>
UnivStudTSampler<float_t>::UnivStudTSampler(float_t dof)
    : rvsamp_base()
    , m_t_gen(dof)
{
}
 
 
template<typename float_t>
float_t UnivStudTSampler<float_t>::sample()
{
    return m_t_gen(m_rng);
}
 
 
 
template<typename float_t>
class UnivLogNormSampler : public rvsamp_base
{
    
public:
 
 
    UnivLogNormSampler();
 
 
    UnivLogNormSampler(float_t mu, float_t sigma);
 
 
    void setSigma(float_t sigma);
    
    
    void setMu(float_t mu);
    
        
    float_t sample();    
    
 
private:
 
    std::normal_distribution<float_t> m_z_gen;
    
    float_t m_mu;
    
    float_t m_sigma;
 
};
 
 
template<typename float_t>
UnivLogNormSampler<float_t>::UnivLogNormSampler()
    : rvsamp_base()
    , m_z_gen(0.0, 1.0)
{
    setMu(0.0);
    setSigma(1.0);
}
 
 
template<typename float_t>
UnivLogNormSampler<float_t>::UnivLogNormSampler(float_t mu, float_t sigma)
    : rvsamp_base()
    , m_z_gen(0.0, 1.0)
{
    setMu(mu); 
    setSigma(sigma);
}
 
 
template<typename float_t>
void UnivLogNormSampler<float_t>::setMu(float_t mu)
{
    m_mu = mu;
}
 
 
template<typename float_t>
void UnivLogNormSampler<float_t>::setSigma(float_t sigma)
{
    m_sigma = sigma;
}
 
 
template<typename float_t>
float_t UnivLogNormSampler<float_t>::sample()
{
    return std::exp(m_mu + m_sigma * m_z_gen(m_rng));
}
 
 
 
template<typename float_t>
class UnivGammaSampler : public rvsamp_base
{
    
public:
 
 
    UnivGammaSampler() = delete;
 
 
    UnivGammaSampler(float_t alpha, float_t beta);
 
 
    float_t sample();    
    
 
private:
 
    std::gamma_distribution<float_t> m_gamma_gen;
    
    float_t m_alpha;
    
    float_t m_beta;
 
};
 
 
template<typename float_t>
UnivGammaSampler<float_t>::UnivGammaSampler(float_t alpha, float_t beta)
    : rvsamp_base()
    , m_gamma_gen(alpha, beta)
{
}
 
 
template<typename float_t>
float_t UnivGammaSampler<float_t>::sample()
{
    return m_gamma_gen(m_rng);
}
 
 
 
template<typename float_t>
class UnivInvGammaSampler : public rvsamp_base
{
    
public:
 
 
    UnivInvGammaSampler();
 
 
    UnivInvGammaSampler(float_t alpha, float_t beta);
 
 
    float_t sample();    
    
 
private:
 
    std::gamma_distribution<float_t> m_gamma_gen;
    
    float_t m_alpha;
    
    float_t m_beta;
 
};
 
 
template<typename float_t>
UnivInvGammaSampler<float_t>::UnivInvGammaSampler()
    : rvsamp_base()
    , m_gamma_gen(1.0, 1.0)
{
}
 
 
template<typename float_t>
UnivInvGammaSampler<float_t>::UnivInvGammaSampler(float_t alpha, float_t beta)
    : rvsamp_base()
    , m_gamma_gen(alpha, beta)
{
}
 
 
template<typename float_t>
float_t UnivInvGammaSampler<float_t>::sample()
{
    return 1.0/m_gamma_gen(m_rng);
}
 
 
 
template<typename float_t>
class TruncUnivNormSampler : public rvsamp_base
{
    
public:
 
 
    TruncUnivNormSampler(float_t mu, float_t sigma, float_t lower, float_t upper);
 
        
    float_t sample();    
    
 
private:
 
    std::normal_distribution<float_t> m_z_gen;
 
    float_t m_mu;
    
    float_t m_sigma;
 
    float_t m_lower;
 
    float_t m_upper;
};
 
 
template<typename float_t>
TruncUnivNormSampler<float_t>::TruncUnivNormSampler(float_t mu, 
                                                    float_t sigma, 
                                                    float_t lower, 
                                                    float_t upper)
    : rvsamp_base()
    , m_z_gen(0.0, 1.0)
    , m_mu(mu)
    , m_sigma(sigma)
    , m_lower(lower)
    , m_upper(upper)
{
}
 
 
template<typename float_t>
float_t TruncUnivNormSampler<float_t>::sample()
{
    float_t proposal;
    bool accepted = false;
    while(!accepted)
    {
        proposal = m_mu + m_sigma*m_z_gen(this->m_rng);
        if((m_lower <= proposal) & (proposal <= m_upper))
            accepted = true;
    }
    return proposal;
}
 
 
 
template<typename float_t, typename int_t>
class PoissonSampler : public rvsamp_base
{
    
public:
 
 
    PoissonSampler();
 
 
    PoissonSampler(float_t lambda);
 
 
    void setLambda(float_t lambda);
    
 
    int_t sample();    
    
 
private:
 
    std::poisson_distribution<int_t> m_p_gen;
};
 
 
template<typename float_t, typename int_t>
PoissonSampler<float_t, int_t>::PoissonSampler() 
    : rvsamp_base(), m_p_gen(float_t(1.0))
{
}
 
 
template<typename float_t, typename int_t>
PoissonSampler<float_t, int_t>::PoissonSampler(float_t lambda) 
    : rvsamp_base(), m_p_gen(lambda)
{
}
 
 
template<typename float_t, typename int_t>
void PoissonSampler<float_t, int_t>::setLambda(float_t lambda)
{
   m_p_gen.param(typename decltype(m_p_gen)::param_type(lambda)); 
}
 
 
template<typename float_t, typename int_t>
int_t PoissonSampler<float_t, int_t>::sample()
{
    return m_p_gen(m_rng);
}
 
 
 
template<typename float_t, typename int_t>
class BernSampler : public rvsamp_base
{
    
public:
 
 
    BernSampler();
 
 
    BernSampler(float_t p);
 
 
    void setP(float_t p);
    
 
    int_t sample();    
    
 
private:
 
    std::bernoulli_distribution m_B_gen;
    
    float_t m_p;
};
 
 
template<typename float_t, typename int_t>
BernSampler<float_t, int_t>::BernSampler() 
    : rvsamp_base(), m_B_gen(.5)
{
}
 
 
template<typename float_t, typename int_t>
BernSampler<float_t, int_t>::BernSampler(float_t p) 
    : rvsamp_base(), m_B_gen(p)
{
}
 
 
template<typename float_t, typename int_t>
void BernSampler<float_t, int_t>::setP(float_t p)
{
    m_p = p;
}
 
 
template<typename float_t, typename int_t>
int_t BernSampler<float_t, int_t>::sample()
{
    return (m_B_gen(m_rng)) ? 1 : 0;
}
 
 
 
 
template<size_t dim, typename float_t>
class MVNSampler : public rvsamp_base
{
public:
 
    using Vec = Eigen::Matrix<float_t,dim,1>;
    using Mat = Eigen::Matrix<float_t,dim,dim>;
    
     MVNSampler();
 
 
    MVNSampler(const Vec &meanVec, const Mat &covMat);
 
 
    void setCovar(const Mat &covMat);
    
    
    void setMean(const Vec &meanVec);
    
        
    auto sample() -> Vec;    
    
private:
 
    std::normal_distribution<float_t> m_z_gen;
    
    Mat m_scale_mat;
    
    Vec m_mean;
    
};
 
 
template<size_t dim, typename float_t>
MVNSampler<dim, float_t>::MVNSampler()
        : rvsamp_base()
        , m_z_gen(0.0, 1.0)
{
    setMean(Vec::Zero());
    setCovar(Mat::Identity());  
}
 
 
template<size_t dim, typename float_t>
MVNSampler<dim, float_t>::MVNSampler(const Vec &meanVec, const Mat &covMat)
    : rvsamp_base()
    , m_z_gen(0.0, 1.0)
{
    setCovar(covMat);
    setMean(meanVec);
}
 
 
template<size_t dim, typename float_t>
void MVNSampler<dim, float_t>::setCovar(const Mat &covMat)
{
    Eigen::SelfAdjointEigenSolver<Mat> eigenSolver(covMat);
    m_scale_mat = eigenSolver.eigenvectors() * eigenSolver.eigenvalues().cwiseMax(0).cwiseSqrt().asDiagonal();
}
 
 
template<size_t dim, typename float_t>
void MVNSampler<dim, float_t>::setMean(const Vec &meanVec)
{
    m_mean = meanVec;
}
 
 
template<size_t dim, typename float_t>
auto MVNSampler<dim, float_t>::sample() -> Vec
{
    Vec Z;
    for (size_t i=0; i< dim; ++i) 
    {
        Z(i) = m_z_gen(this->m_rng);
    }
    return m_mean + m_scale_mat * Z;
}
 
 
 
 
template<typename float_t>
class UniformSampler : public rvsamp_base
{
public:
 
    UniformSampler();
    
    
    UniformSampler(float_t lower, float_t upper);
    
    
    float_t sample();
 
private:
 
    std::uniform_real_distribution<float_t> m_unif_gen;
 
};
 
 
template<typename float_t>
UniformSampler<float_t>::UniformSampler() 
        : rvsamp_base()
        , m_unif_gen(0.0, 1.0)
{
}
 
 
template<typename float_t>
UniformSampler<float_t>::UniformSampler(float_t lower, float_t upper) 
        : rvsamp_base()
        , m_unif_gen(lower, upper)
{    
}
 
 
template<typename float_t>
float_t UniformSampler<float_t>::sample()
{
    return m_unif_gen(m_rng);
}
 
 
 
 
template<size_t N, typename float_t>
class k_gen : public rvsamp_base
{
public:
    k_gen();
 
    
    std::array<unsigned int, N> sample(const std::array<float_t, N> &logWts);     
};
 
 
template<size_t N, typename float_t>
k_gen<N, float_t>::k_gen() : rvsamp_base() {}
 
 
template<size_t N, typename float_t>
std::array<unsigned int, N> k_gen<N, float_t>::sample(const std::array<float_t, N> &logWts)
{
    // these log weights may be very negative. If that's the case, exponentiating them may cause underflow
    // so we use the "log-exp-sum" trick
    // actually not quite...we just shift the log-weights because after they're exponentiated
    // they have the same normalized probabilities
    
   // Create the distribution with exponentiated log-weights
   // subtract the max first to prevent underflow
   // normalization is taken care of by std::discrete_distribution
    std::array<float_t, N> w;
    float_t m = *std::max_element(logWts.begin(), logWts.end());
    std::transform(logWts.begin(), logWts.end(), w.begin(), 
                   [&m](float_t d) -> float_t { return std::exp(d-m); } );
    std::discrete_distribution<> kGen(w.begin(), w.end());
    
    // sample and return ks
    std::array<unsigned int, N> ks;
    for(size_t i = 0; i < N; ++i){
        ks[i] = kGen(this->m_rng);
    }
    return ks;
}
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
template<typename float_t>
class BetaSampler : public rvsamp_base
{
    
public:
 
 
    BetaSampler() = delete;
 
 
    BetaSampler(float_t alpha, float_t beta);
 
 
    float_t sample();    
    
 
private:
 
    std::gamma_distribution<float_t> m_first_gamma_gen;
    
    std::gamma_distribution<float_t> m_second_gamma_gen;
  
    float_t m_alpha;
    
    float_t m_beta;
 
};
 
 
template<typename float_t>
BetaSampler<float_t>::BetaSampler(float_t alpha, float_t beta)
    : rvsamp_base()
    , m_first_gamma_gen(alpha, 1.0)
    , m_second_gamma_gen(beta, 1.0)
{
}
 
 
template<typename float_t>
float_t BetaSampler<float_t>::sample()
{
    float_t first = m_first_gamma_gen(m_rng);
    float_t second = m_second_gamma_gen(m_rng);
    return first / (first + second);
}
 
 
} // namespace rv_samp
    
} // namespace pf
    
    
#endif // RV_SAMP_H