Subversion Repositories public iLand

#include <cmath>
#include <stdexcept>
#include <sstream>
#include <iostream>
#include "SimpleRNG.h"

#define PI 3.1415926535897932384626433832795

SimpleRNG::SimpleRNG()
{
    // These values are not magical, just the default values Marsaglia used.
    // Any unit should work.
    m_u = 521288629;
    m_v = 362436069;
}

void SimpleRNG::SetState(unsigned int u, unsigned int v)
{
    m_u = u;
    m_v = v;
}

void SimpleRNG::GetState(unsigned int& u, unsigned int& v)
{
    u = m_u;
    v = m_v;
}

double SimpleRNG::GetUniform(unsigned int& u, unsigned int& v)
{
    // 0 <= u <= 2^32
    unsigned int z = GetUint(u, v);
    // The magic number is 1/(2^32 + 1) and so result is positive and less than 1.
    return z*2.328306435996595e-10;
}

unsigned int SimpleRNG::GetUint(unsigned int& u, unsigned int& v)
{
    v = 36969*(v & 65535) + (v >> 16);
    u = 18000*(u & 65535) + (u >> 16);
    return (v << 16) + u;
}

double SimpleRNG::GetUniform()
{
    return GetUniform(m_u, m_v);
}

unsigned int SimpleRNG::GetUint()
{
    return GetUint(m_u, m_v);
}
   
// Get normal (Gaussian) random sample with specified mean and standard deviation
double SimpleRNG::GetNormal(double mean = 0.0, double standardDeviation = 1.0)
{
    if (standardDeviation <= 0.0)
    {
                std::stringstream os;
        os << "Standard deviation must be positive." << "\n"
           << "Received standard deviation " << standardDeviation;
        throw std::invalid_argument( os.str() );
    }
            // Use Box-Muller algorithm
    double u1 = GetUniform();
    double u2 = GetUniform();
    double r = sqrt( -2.0*log(u1) );
    double theta = 2.0*PI*u2;
    return mean + standardDeviation*r*sin(theta);
}
       
// Get exponential random sample with specified mean
double SimpleRNG::GetExponential(double mean = 1.0)
{
    if (mean <= 0.0)
    {
                std::stringstream os;
        os << "Exponential mean must be positive." << "\n"
           << "Received standard deviation " << mean;
        throw std::invalid_argument( os.str() );
    }
    return -mean*log( GetUniform() );
}

double SimpleRNG::GetGamma(double shape, double scale)
{
    // Implementation based on "A Simple Method for Generating Gamma Variables"
    // by George Marsaglia and Wai Wan Tsang.  ACM Transactions on Mathematical Software
    // Vol 26, No 3, September 2000, pages 363-372.

    double d, c, x, xsquared, v, u;

    if (shape >= 1.0)
    {
        d = shape - 1.0/3.0;
        c = 1.0/sqrt(9.0*d);
        for (;;)
        {
            do
            {
                x = GetNormal();
                v = 1.0 + c*x;
            }
            while (v <= 0.0);
            v = v*v*v;
            u = GetUniform();
            xsquared = x*x;
            if (u < 1.0 -.0331*xsquared*xsquared || log(u) < 0.5*xsquared + d*(1.0 - v + log(v)))
                return scale*d*v;
        }
    }
    else if (shape <= 0.0)
    {
                std::stringstream os;
        os << "Shape parameter must be positive." << "\n"
           << "Received shape parameter " << shape;
        throw std::invalid_argument( os.str() );
    }
    else
    {
        double g = GetGamma(shape+1.0, 1.0);
        double w = GetUniform();
        return scale*g*pow(w, 1.0/shape);
    }
}

double SimpleRNG::GetChiSquare(double degreesOfFreedom)
{
    // A chi squared distribution with n degrees of freedom
    // is a gamma distribution with shape n/2 and scale 2.
    return GetGamma(0.5 * degreesOfFreedom, 2.0);
}

double SimpleRNG::GetInverseGamma(double shape, double scale)
{
    // If X is gamma(shape, scale) then
    // 1/Y is inverse gamma(shape, 1/scale)
    return 1.0 / GetGamma(shape, 1.0 / scale);
}

double SimpleRNG::GetWeibull(double shape, double scale)
{
    if (shape <= 0.0 || scale <= 0.0)
    {
                std::stringstream os;
        os << "Shape and scale parameters must be positive." << "\n"
           << "Received shape " << shape << " and scale " << scale;
        throw std::invalid_argument( os.str() );
    }
    return scale * pow(-log(GetUniform()), 1.0 / shape);
}

double SimpleRNG::GetCauchy(double median, double scale)
{
    if (scale <= 0)
    {
                std::stringstream os;
                os << "Shape parameter must be positive." << "\n"
           << "Received shape parameter " << scale;
        throw std::invalid_argument( os.str() );
    }

    double p = GetUniform();

    // Apply inverse of the Cauchy distribution function to a uniform
    return median + scale*tan(PI*(p - 0.5));
}

double SimpleRNG::GetStudentT(double degreesOfFreedom)
{
    if (degreesOfFreedom <= 0)
    {
                std::stringstream os;
                os << "Degrees of freedom must be positive." << "\n"
           << "Received shape parameter " << degreesOfFreedom;
        throw std::invalid_argument( os.str() );
    }

    // See Seminumerical Algorithms by Knuth
    double y1 = GetNormal();
    double y2 = GetChiSquare(degreesOfFreedom);
    return y1 / sqrt(y2 / degreesOfFreedom);
}

// The Laplace distribution is also known as the double exponential distribution.
double SimpleRNG::GetLaplace(double mean, double scale)
{
    double u = GetUniform();
    return (u < 0.5) ?
        mean + scale*log(2.0*u) :
        mean - scale*log(2*(1-u));
}

double SimpleRNG::GetLogNormal(double mu, double sigma)
{
    return exp(GetNormal(mu, sigma));
}

double SimpleRNG::GetBeta(double a, double b)
{
    if (a <= 0.0 || b <= 0.0)
    {
                std::stringstream os;
                os << "Beta parameters must be positive." << "\n"
           << "Received parameters " << a << " and " << b;
        throw std::invalid_argument( os.str() );
    }

    // There are more efficient methods for generating beta samples.
    // However such methods are a little more efficient and much more complicated.
    // For an explanation of why the following method works, see
    // http://www.johndcook.com/distribution_chart.html#gamma_beta

    double u = GetGamma(a, 1.0);
    double v = GetGamma(b, 1.0);
    return u / (u + v);
}

int SimpleRNG::GetPoisson(double lambda)
{
    return (lambda < 30.0) ? PoissonSmall(lambda) : PoissonLarge(lambda);
}

int SimpleRNG::PoissonSmall(double lambda)
{
    // Algorithm due to Donald Knuth, 1969.
    double p = 1.0, L = exp(-lambda);
        int k = 0;
        do
        {
                k++;
                p *= GetUniform();
        }
        while (p > L);
    return k - 1;
}

int SimpleRNG::PoissonLarge(double lambda)
{
        // "Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson
        // Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979)
        // The article is on pages 29-35. The algorithm given here is on page 32.

        double c = 0.767 - 3.36/lambda;
        double beta = PI/sqrt(3.0*lambda);
        double alpha = beta*lambda;
        double k = log(c) - lambda - log(beta);

        for(;;)
        {
                double u = GetUniform();
                double x = (alpha - log((1.0 - u)/u))/beta;
                int n = (int) floor(x + 0.5);
                if (n < 0)
                        continue;
                double v = GetUniform();
                double y = alpha - beta*x;
        double temp = 1.0 + exp(y);
                double lhs = y + log(v/(temp*temp));
                double rhs = k + n*log(lambda) - LogFactorial(n);
                if (lhs <= rhs)
                        return n;
        }
}

double SimpleRNG::LogFactorial(int n)
{

    if (n < 0)
    {
        std::stringstream os;
        os << "Invalid input argument (" << n
        << "); may not be negative"