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"