[Boost-commit] svn:boost r42757 - sandbox/random/boost/random

From: asutton_at_[hidden]
Date: 2008-01-14 08:26:57

Author: asutton
Date: 2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
New Revision: 42757
URL: http://svn.boost.org/trac/boost/changeset/42757

Log:
Added a bunch of new random distributions. Note that this commit also includes "fixes"
for older distributions that may not have supported enough parameters or computed values
in slightly "odd" ways (lognormal and gamma).

Added:
   sandbox/random/boost/random/beta_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/chi_squared_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/distributions.hpp   (contents, props changed)
   sandbox/random/boost/random/extreme_value_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/fisher_f_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/gamma_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/lognormal_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/pareto_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/rayleigh_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/skew_normal_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/students_t_distribution.hpp   (contents, props changed)
   sandbox/random/boost/random/weibull_distribution.hpp   (contents, props changed)

Added: sandbox/random/boost/random/beta_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/beta_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,119 @@
+// Copyright Jeremy Bruestle 2007
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_BETA_DISTRIBUTION_HPP
+#define BOOST_RANDOM_BETA_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+
+namespace boost {
+
+    // The Beta Distribution for random number generation as implemented by
+    // Jeremy Bruestle.
+    //
+    // This implementation is taken from the paper Cheng, "Generating Beta
+    // Variates with Nonintegral Shape Parameters", Communications of
+    // the ACM, 21(4), 1978.
+    template<class RealType = double>
+    class beta_distribution
+    {
+    public:
+        typedef RealType input_type;
+        typedef RealType result_type;
+
+        beta_distribution(const result_type& a = result_type(1),
+                        const result_type& b = result_type(1))
+            : _a(a)
+            , _b(b)
+            , _alpha()
+            , _beta()
+            , _theta()
+        {
+            BOOST_ASSERT(a > result_type(0));
+            BOOST_ASSERT(b > result_type(0));
+            init();
+        }
+
+        inline result_type a() const
+        { return _a; }
+
+        inline result_type b() const
+        { return _b; }
+
+        inline void reset()
+        { }
+
+        template<class Engine>
+        result_type operator ()(Engine& eng)
+        {
+            using std::log;
+            using std::exp;
+
+            // Try until we get a working case: expected number attempts varies
+            // based on a + b worst case number of expected attempts is 4.
+            result_type u1, u2, v, w, x;
+            do {
+                u1 = eng();
+                u2 = eng();
+                v = _beta * log (u1 / (result_type(1) - u1));
+                w = _a * exp(v);
+                x = _alpha * log(_alpha / (_b + w)) + _theta * v - log(result_type(4));
+            } while(x < log(u1 * u1 * u2));
+
+            return w / (_b + w);
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT, Traits>&
+        operator<<(std::basic_ostream<CharT, Traits>& os, const beta_distribution& dist)
+        {
+            os << dist.m_a << dist.m_b;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT, Traits>&
+        operator>>(std::basic_istream<CharT,Traits>& is, beta_distribution& dist)
+        {
+            is >> std::ws >> dist.m_a >> std::ws >> dist.m_b;
+            dist.init();
+            return is;
+        }
+#endif
+
+    private:
+        void init()
+        {
+            using std::exp;
+            using std::log;
+            using std::sqrt;
+            using std::max;
+
+            _alpha = _a + _b;
+            if (_a <= result_type(1) || _b <= result_type(1)) {
+                _beta = max(result_type(1) / _a, result_type(1) / _b);
+            }
+            else {
+                _beta = sqrt((_alpha - result_type(2)) / (result_type(2) * _a *_b - _alpha));
+            }
+            _theta = _a + result_type(1) / _beta;
+        }
+
+        result_type _a;
+        result_type _b;
+        result_type _alpha;
+        result_type _beta;
+        result_type _theta;
+    };
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_BETA_DISTRIBUTION_HPP

Added: sandbox/random/boost/random/chi_squared_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/chi_squared_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,75 @@
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_CHI_SQUARED_DISTRIBUTION_HPP
+#define BOOST_RANDOM_CHI_SQUARED_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+#include <boost/random/gamma_distribution.hpp>
+
+namespace boost {
+
+    // The Chi-Squared distribution...
+    //
+    // This implementation is drawn from the fact that the chi-squared is a
+    // special case of the gamma distribution with alpha being the degrees of
+    // freedom and theta, 2.
+    template<class RealType = double>
+    class chi_squared_distribution
+    {
+        typedef gamma_distribution<RealType> base_type;
+    public:
+        typedef RealType result_type;
+
+        explicit chi_squared_distribution(const result_type& df = result_type(1))
+            : m_gamma(df / result_type(2), result_type(2))
+        {
+            BOOST_ASSERT(degrees_of_freedom() > 0);
+        }
+
+        inline result_type degrees_of_freedom() const
+        { return m_gamma.alpha() * result_type(2); }
+
+        void reset()
+        { m_gamma.reset(); }
+
+        // Note that eng must provide values in the range [0, 1).
+        template<class Engine>
+        result_type operator ()(Engine& eng)
+        {
+            return m_gamma(eng);
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT, Traits>&
+        operator <<(std::basic_ostream<CharT, Traits>& os, const chi_squared_distribution& dist)
+        {
+            os << dist.m_gamma;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT, Traits>&
+        operator >>(std::basic_istream<CharT,Traits>& is, chi_squared_distribution& dist)
+        {
+            is >> std::ws >> dist.m_gamma;
+            return is;
+        }
+#endif
+
+    private:
+
+        result_type m_df;
+        base_type m_gamma;
+    };
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP

Added: sandbox/random/boost/random/distributions.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/distributions.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,29 @@
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_DISTRIBUTIONS_HPP
+#define BOOST_RANDOM_DISTRIBUTIONS_HPP
+
+#include <boost/random/bernoulli_distribution.hpp>
+#include <boost/random/binomial_distribution.hpp>
+#include <boost/random/poisson_distribution.hpp>
+#include <boost/random/geometric_distribution.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/lognormal_distribution.hpp>
+#include <boost/random/skew_normal_distribution.hpp>
+#include <boost/random/exponential_distribution.hpp>
+#include <boost/random/gamma_distribution.hpp>
+#include <boost/random/beta_distribution.hpp>
+#include <boost/random/chi_squared_distribution.hpp>
+#include <boost/random/fisher_f_distribution.hpp>
+#include <boost/random/students_t_distribution.hpp>
+#include <boost/random/cauchy_distribution.hpp>
+#include <boost/random/extreme_value_distribution.hpp>
+#include <boost/random/weibull_distribution.hpp>
+#include <boost/random/rayleigh_distribution.hpp>
+#include <boost/random/pareto_distribution.hpp>
+
+#endif // BOOST_RANDOM_DISTRIBUTIONS_HPP

Added: sandbox/random/boost/random/extreme_value_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/extreme_value_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,93 @@
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_EXTREME_VALUE_DISTRIBUTION_HPP
+#define BOOST_RANDOM_EXTREME_VALUE_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+
+namespace boost {
+
+    // The Gumbel distribution is one of the extreme value (I) distributions
+    // along with Weibull and Frechet. This is also known as the Fisher-Tippet
+    // distribution or the log-Weibull.
+    //
+    // This distribution actually implements a class of several different
+    // distributions. The extreme value distribution deals with the median
+    // of different distributions. There are three distinct flavors of extreme
+    // value distributions (I, II, and III). For type I distributions, there
+    // are two cases: minimum and maximum. Generally when authors refer to the
+    // type I distribution, they refer to the maximum case. The different cases
+    // only differ by a factor of -1 in most cases.
+    //
+    // Note that the type II and III are closely related to the minimum and
+    // maximum cases of the type I.
+    template<class RealType = double>
+    class extreme_value_distribution
+    {
+    public:
+        typedef RealType result_type;
+
+        explicit extreme_value_distribution(const result_type& mu = result_type(0),
+                                            const result_type& beta = result_type(1),
+                                            const result_type& sign = result_type(1))
+            : m_mu(mu)
+            , m_beta(beta)
+            , m_sign(sign)
+        {
+            BOOST_ASSERT(m_beta > result_type(0));
+            BOOST_ASSERT(m_sign == -1 || m_sign == 1);
+        }
+
+        result_type mu() const
+        { return m_mu; }
+
+        result_type beta() const
+        { return m_beta; }
+
+        result_type sign() const
+        { return m_sign; }
+
+        void reset()
+        { }
+
+        // Note that eng must provide values in the range (0, 1].
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            using std::log;
+            return m_mu - m_sign * m_beta * log(-log(eng()));
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT, Traits>&
+        operator <<(std::basic_ostream<CharT, Traits>& os, const extreme_value_distribution& dist)
+        {
+            os << dist.m_mu << " " << dist.m_beta << " " << dist.m_sign;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT,Traits>&
+        operator >>(std::basic_istream<CharT,Traits>& is, extreme_value_distribution& dist)
+        {
+            is >> std::ws >> dist.m_mu >> std::ws >> dist.m_beta >> std::ws >> dist.m_sign;
+            return is;
+        }
+#endif
+
+        result_type m_mu;
+        result_type m_beta;
+        result_type m_sign;
+    };
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP

Added: sandbox/random/boost/random/fisher_f_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/fisher_f_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,82 @@
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_F_DISTRIBUTION_HPP
+#define BOOST_RANDOM_F_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+#include <boost/random/chi_squared_distribution.hpp>
+
+namespace boost {
+
+    // The F distribution (a.k.a., Snedecor's F or Fisher-Snedecor) represents
+    // the ration of two chi-squared variates.
+    template<class RealType = double>
+    class fisher_f_distribution
+    {
+        typedef chi_squared_distribution<RealType> base_type;
+
+    public:
+        typedef RealType result_type;
+
+        explicit fisher_f_distribution(const result_type& d1 = result_type(1),
+                                       const result_type& d2 = result_type(1))
+            : m_first(d1)
+            , m_second(d2)
+        {
+            BOOST_ASSERT(first_degrees_of_freedom() > 0);
+            BOOST_ASSERT(second_degrees_of_freedom() > 0);
+        }
+
+        inline result_type first_degrees_of_freedom() const
+        { return m_first.degrees_of_freedom(); }
+
+        inline result_type second_degrees_of_freedom() const
+        { return m_second.degrees_of_freedom(); }
+
+        void reset()
+        {
+            m_first.reset();
+            m_second.reset();
+        }
+
+        // Note that eng must provide values in the range (0, 1].
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            result_type top = m_first(eng) / first_degrees_of_freedom();
+            result_type bot = m_second(eng) / second_degrees_of_freedom();
+            return top / bot;
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT, Traits>&
+        operator <<(std::basic_ostream<CharT, Traits>& os, const fisher_f_distribution& dist)
+        {
+            os << dist.m_first << " " << dist.m_second;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT,Traits>&
+        operator >>(std::basic_istream<CharT,Traits>& is, fisher_f_distribution& dist)
+        {
+            is >> std::ws >> dist.m_first << std::ws >> dist.m_second;
+            return is;
+        }
+#endif
+
+        base_type m_first;
+        base_type m_second;
+    };
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP

Added: sandbox/random/boost/random/gamma_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/gamma_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,153 @@
+// Copyright Jens Maurer 2000-2001
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
+#define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+
+namespace boost {
+
+    // A rewrite of the old gamma distribution, this adds support for the
+    // two-parameter version described almost everywhere else (Wikipedia,
+    // MathWorld, etc). The first parameter is named alpha (or k), the second
+    // is named theta (or beta).
+    //
+    // The gamma distribution represents the summation of alpha exponentially
+    // distributed random variables,each with the parameter theta.
+    template<class RealType = double>
+    class gamma_distribution
+    {
+    public:
+        typedef RealType result_type;
+
+        explicit gamma_distribution(const result_type& alpha = result_type(1),
+                                    const result_type& theta = result_type(1))
+            : m_alpha(alpha)
+            , m_theta(theta)
+        {
+            BOOST_ASSERT(m_alpha > result_type(0));
+            BOOST_ASSERT(m_theta > result_type(0));
+        }
+
+        result_type alpha() const
+        { return m_alpha; }
+
+        result_type theta() const
+        { return m_theta; }
+
+        void reset()
+        { }
+
+        // Note that eng must provide values in the range [0, 1).
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            using std::sqrt;
+            using std::pow;
+            using std::log;
+            using std::exp;
+            using std::tan;
+
+            if(m_alpha == result_type(1)) {
+                // Pick a number, not 0, and return -log(u) * theta
+                result_type u;
+                for(u = eng(); u == 0; u = eng()) ;
+                return -log(u) * m_theta;
+            }
+            else if(m_alpha > result_type(1)) {
+                // This part of the algorithm appears to come from Cheng'1977,
+                // "The generation of Gamma variables with non-integral shape
+                // parameters", in Applied Statistices, 26(1), pp. 71-74.
+                // At least that's what the Python code sites :)
+
+                // Preliminaries
+                result_type a = sqrt(result_type(2) * m_alpha - result_type(1));
+                result_type b = m_alpha - log(result_type(4));
+                result_type c = m_alpha + a;
+
+                while(1) {
+                    // Pick two numbers s.t. 0 != u1 != 1. Note that the Python
+                    // implementation constrains these with a precision of 1e-7
+                    // away from the values.
+                    result_type u1 = eng();
+                    if(u1 == 0 || u1 == 1) {
+                        continue;
+                    }
+                    result_type u2 = result_type(1) - eng();
+
+                    // A bunch of computing. I don't know what this does.
+                    result_type v = log(u1 / (result_type(1) - u1)) / a;
+                    result_type x = m_alpha * exp(v);
+                    result_type z = u1 * u1 * u2;
+                    result_type r = b + c * v - x;
+                    result_type magic = result_type(1) + log(result_type(4.5));
+
+                    if((r + magic - result_type(4.5) * z >= 0.0) ||
+                       (r >= log(z)))
+                    {
+                        return x * m_theta;
+                    }
+                }
+            }
+            else {
+                // For lack of a constant...
+                result_type e = exp(1);
+
+                // For alpha < 0, we use the ALGORITHM GS of Statistical Computing
+                // by Kennedy and Gentle. I didn't read this either, but did basically
+                // reimplement the Python version.
+                result_type x;
+                while(1) {
+                    result_type u = eng();
+                    result_type b = (m_alpha + e) / e;
+                    result_type p = b * u;
+                    x = p < result_type(1)
+                            ? pow(p, result_type(1) / m_alpha)
+                            : -log((b - p) / m_alpha);
+                    result_type u1 = eng();
+                    if(p > result_type(1)) {
+                        if(u1 <= pow(x, m_alpha - result_type(1))) {
+                            break;
+                        }
+                    }
+                    else if(u1 <= exp(-x)) {
+                        break;
+                    }
+                }
+                return x * m_theta;
+            }
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT, Traits>&
+        operator<<(std::basic_ostream<CharT, Traits>& os, const gamma_distribution& dist)
+        {
+            os << dist.m_alpha << " " << dist.m_theta;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT,Traits>&
+        operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& dist)
+        {
+            is >> std::ws >> dist.m_alpha << std::ws >> dist.m_theta;
+            return is;
+        }
+#endif
+
+        result_type m_alpha;
+        result_type m_theta;
+    };
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP

Added: sandbox/random/boost/random/lognormal_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/lognormal_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,84 @@
+// Copyright Jens Maurer 2000-2001
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP
+#define BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+#include <boost/random/normal_distribution.hpp>
+
+namespace boost {
+    // A rewrite of the original lognormal dsitribution. This removes the work
+    // done to maintain an adjusted mean and sigma and simply exponentiates
+    // values of a normal distribution. This implementation completely removes
+    // the artificial constraints on the mean (couldn't be 0). Also, it's a
+    // bit cleaner now.
+    template<typename RealType = double>
+    class lognormal_distribution
+    {
+        typedef normal_distribution<RealType> base_type;
+    public:
+        typedef RealType result_type;
+
+        explicit lognormal_distribution(result_type mean = result_type(1),
+                                         result_type sigma = result_type(1))
+            : m_norm(mean, sigma)
+        {
+            BOOST_ASSERT(sigma > 0);
+        }
+
+        inline result_type mean() const
+        { return m_norm.mean(); }
+
+        inline result_type sigma() const
+        { return m_norm.sigma(); }
+
+        result_type reset()
+        { m_norm.reset(); }
+
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            using std::exp;
+
+            // The previous implementation scales and translates by some
+            // kind of a factor. Why don't we just exponentiate a normal
+            // variate with the given mean and sigma.
+            //
+            // return exp(_norm(eng) * _nsigma + _nmean);
+
+            return exp(m_norm(eng));
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT, Traits>&
+        operator<<(std::basic_ostream<CharT, Traits>& os, const lognormal_distribution& dist)
+        {
+            os << dist.m_normal;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT, Traits>&
+        operator>>(std::basic_istream<CharT, Traits>& is, lognormal_distribution& dist)
+        {
+            is >> dist.m_normal;
+            return is;
+        }
+#endif
+
+    private:
+        base_type m_norm;
+    };
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP

Added: sandbox/random/boost/random/pareto_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/pareto_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,91 @@
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_PARETO_DISTRIBUTION_HPP
+#define BOOST_RANDOM_PARETO_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+
+namespace boost {
+
+    // The Pareto distribution can be used to Pareto-distributed random numbers.
+    // This is basically the set of functions similar to 1/x. The Pareto
+    // distribution is given in terms of a location, xi and a shape, alpha. The
+    // implementation of this class comes from information gleaned from
+    // Wikipedia. Note that MathWorld uses k as the location, whereas Wikipedia
+    // uses k as the shape. We prefer to use k to denote location.
+    //
+    // Note that there is also a generalized Pareto distribution. We should
+    // probably implement this in terms of that, but we don't need to. Moreover,
+    // this matches the Pareto distribution given in the Math toolkit.
+    template <typename ResultType = double>
+    class pareto_distribution
+    {
+    public:
+        typedef ResultType result_type;
+
+        explicit pareto_distribution(const result_type& location = result_type(1),
+                                     const result_type& shape = result_type(1))
+            : m_location(location)
+            , m_shape(shape)
+        {
+            BOOST_ASSERT(m_location > 0);
+            BOOST_ASSERT(m_shape > 0);
+        }
+
+        // Returns the minimum value of the distribution.
+        result_type location() const
+        { return m_location; }
+
+        // Returns the shape of the distribution.
+        result_type shape() const
+        { return m_shape; }
+
+        void reset()
+        { }
+
+        // Note that Engine must generate uniformly distributed random numbers
+        // in the range [0, 1).
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            using std::pow;
+
+            // This is actually taken from the Wikipedia section on sampling
+            // rather than the generalized Pareto generator.
+            return m_location / pow(eng(), result_type(1) / m_shape);
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT,Traits>&
+        operator<<(std::basic_ostream<CharT,Traits>& os, const pareto_distribution& dist)
+        {
+            os << dist.m_location << " " << dist.m_shape;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT, Traits>&
+        operator>>(std::basic_istream<CharT, Traits>& is, pareto_distribution& dist)
+        {
+            is >> std::ws >> dist.m_location >> std::ws >> dist.m_shape;
+            return is;
+        }
+#endif
+
+    private:
+        result_type m_location;
+        result_type m_shape;
+    };
+
+} // namespace boost
+
+#endif
+

Added: sandbox/random/boost/random/rayleigh_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/rayleigh_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,102 @@
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_RAYLEIGH_DISTRIBUTION_HPP
+#define BOOST_RANDOM_RAYLEIGH_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+#include <boost/random/normal_distribution.hpp>
+
+namespace boost {
+
+    // They Rayleigh distribution generates random numbers distributed along
+    // a Rayleigh distribution with the parameter sigma. The implentation of
+    // this generator is taken from Wikipedia. Specifically R ~ Rayleigh(sigma)
+    // if R = sqrt(X^2 + Y^2) where X, Y ~ N(0, sigma^2).
+    //
+    // Implementation-wise, we actually use two different distribution objects
+    // to represent the random vectors upon which the distribution is based.
+    // I'm not sure if we actually need two since they're both have the same
+    // distribution. Just something to check.
+    template <typename ResultType = double>
+    class rayleigh_distribution
+    {
+        typedef normal_distribution<ResultType> base_type;
+    public:
+        typedef ResultType result_type;
+
+        // Default sigma is 1 - which looks like the maximum of the curve
+        // defined by the distribution should also be 1.
+        explicit rayleigh_distribution(const result_type& sigma = 1)
+            : m_sigma(sigma)
+            , m_x(0, sigma)
+            , m_y(0, sigma)
+        {
+            assert(sigma > result_type(0));
+        }
+
+        rayleigh_distribution(const rayleigh_distribution& dist)
+            : m_sigma(dist.m_sigma)
+            , m_x(0, dist.m_sigma * dist.m_sigma)
+            , m_y(0, dist.m_sigma * dist.m_sigma)
+        { }
+
+        result_type sigma() const
+        { return m_sigma; }
+
+        void reset()
+        {
+            m_x.reset();
+            m_y.reset();
+        }
+
+        // Note that Engine must generate uniformly distributed random numbers
+        // in the range (0, 1).
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            using std::sqrt;
+
+            // For this algorithm, draw a number and square it. We can't use
+            // two random values for each x and y since we can forseeably (and
+            // it will happen) end up with negative numbers under the root.
+            // Using only one random number will prevent that.
+            result_type x = m_x(eng);
+            result_type y = m_y(eng);
+            return sqrt(x * x + y * y);
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT,Traits>&
+        operator<<(std::basic_ostream<CharT,Traits>& os, const rayleigh_distribution& wd)
+        {
+            os << wd._sigma;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT,Traits>&
+        operator>>(std::basic_istream<CharT,Traits>& is, rayleigh_distribution& wd)
+        {
+            is >> std::ws >> wd._sigma;
+            return is;
+        }
+#endif
+
+    private:
+        result_type m_sigma;
+        base_type m_x;
+        base_type m_y;
+    };
+
+} // namespace boost
+
+#endif
+

Added: sandbox/random/boost/random/skew_normal_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/skew_normal_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,130 @@
+// Copyright Edward M. Morrison 2007
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_SKEW_NORMAL_DISTRIBUTION_HPP
+#define BOOST_RANDOM_SKEW_NORMAL_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+#include <boost/random/normal_distribution.hpp>
+
+namespace boost
+{
+    // Skew normal distribution based on descriptions of Adelchi Azzalini given
+    // here. http://azzalini.stat.unipd.it/SN.
+    //
+    // This distribution is constructed over a standard normal distribution object
+    // that provides a pair of random numbers for skewing. The inputs are as
+    // follows:
+    //  - mean - The center of the distribution.
+    //  - sigma - The standard deviation of the distribution.
+    //  - alpha - Shape parameter. Negative values skew to the left, positive
+    //            to the right. Note that a shape of 0 is just a normal
+    //            distribution.
+    //
+    // Note that the alpha parameter is used to compute a delta value that
+    // controls the correlation between two random numbers. This value is
+    // cached inside the distribution to avoid redundant computation.
+    template <typename RealType = double>
+    class skew_normal_distribution
+    {
+        typedef normal_distribution<RealType> base_type;
+    public:
+        typedef RealType result_type;
+
+        skew_normal_distribution(const result_type& mean = result_type(0),
+                                 const result_type& sigma = result_type(1),
+                                 const result_type& alpha = result_type(0))
+            : m_mean(mean)
+            , m_sigma(sigma)
+            , m_alpha(alpha)
+            , m_delta(make_delta(alpha))
+            , m_norm()
+        {
+            BOOST_ASSERT(m_sigma >= result_type(0));
+            BOOST_ASSERT(m_delta >= result_type(-1) && m_delta <= result_type(1));
+        }
+
+        result_type mean() const
+        { return m_mean; }
+
+        result_type sigma() const
+        { return m_sigma; }
+
+        result_type delta() const
+        { return m_delta; }
+
+        result_type alpha() const
+        { return m_alpha; }
+
+
+        void reset()
+        { m_norm.reset(); }
+
+        // Note that the engine must generate uniform random numbers in the
+        // range of (0, 1).
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            using std::sqrt;
+
+            // Generate a couple of base random numbers.
+            result_type u0 = m_norm(eng);
+            result_type v = m_norm(eng);
+
+            // Compute u1 such that u0, u1 have a correlation of delta.
+            result_type one = result_type(1);
+            result_type u1 =  u0 * m_delta + v * sqrt(one - m_delta * m_delta);
+            result_type z = (u0 >= 0 ) ? u1 : -u1;
+
+            // Scale and tranlate the result so that it conforms to the mean
+            // and standard deviation of the distribution.
+            return m_mean + m_sigma * z;
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT,Traits>&
+        operator <<(std::basic_ostream<CharT,Traits>& os, const skew_normal_distribution& dist)
+        {
+            os << dist.m_mean << " " << dist.m_sigma << " "
+               << dist.m_delta << " " << dist.m_delta;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT,Traits>&
+        operator >>(std::basic_istream<CharT,Traits>& is, skew_normal_distribution& dist)
+        {
+            is >> std::ws >> dist.m_mean >> std::ws >> dist.m_sigma
+               >> std::ws >> dist.m_alpha >> std::ws >> dist.m_delta;
+            return is;
+        }
+#endif
+
+    private:
+        result_type make_delta(result_type alpha)
+        {
+            using std::sqrt;
+            result_type one(1);
+            return alpha / sqrt(one + alpha * alpha);
+        }
+
+    private:
+        result_type m_mean;
+        result_type m_sigma;
+        result_type m_alpha;
+        result_type m_delta;
+        base_type m_norm;
+    };
+
+} // namespace boost
+
+#endif
+

Added: sandbox/random/boost/random/students_t_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/students_t_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,81 @@
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_STUDENTS_T_DISTRIBUTION_HPP
+#define BOOST_RANDOM_STUDENTS_T_DISTRIBUTION_HPP
+
+#include <cmath>
+#include <iostream>
+
+#include <boost/assert.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/chi_squared_distribution.hpp>
+
+namespace boost {
+
+    // The Student's T distribution...
+    //
+    // The implementation of this distribution is taken from Wikipedia, which
+    // describes the ratio of two random variates (a Normal and Chi-Square)
+    // being t-distributed over v diegrees of freedom. The normal distribution
+    // is the standard normal distributuion, and the Chi-Square distribution
+    // has v degrees of freedom.
+    template<class RealType = double>
+    class students_t_distribution
+    {
+        typedef normal_distribution<RealType> normal_type;
+        typedef chi_squared_distribution<RealType> chi_square_type;
+    public:
+        typedef RealType result_type;
+
+        explicit students_t_distribution(const result_type& df = result_type(1))
+            : m_norm()
+            , m_chi(df)
+        {
+            BOOST_ASSERT(df > 0);
+        }
+
+        inline result_type degrees_of_freedom() const
+        { return m_chi.degrees_of_freedom(); }
+
+        void reset()
+        {
+            m_norm.reset();
+            m_chi.reset();
+        }
+
+        // Note that eng must provide values in the range (0, 1].
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            return m_norm(eng) / sqrt(m_chi(eng) / degrees_of_freedom());
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT, Traits>&
+        operator <<(std::basic_ostream<CharT, Traits>& os, const students_t_distribution& dist)
+        {
+            os << dist.m_norm << " " << dist.m_chi;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT,Traits>&
+        operator >>(std::basic_istream<CharT,Traits>& is, students_t_distribution& dist)
+        {
+            is >> std::ws >> dist.m_norm << std::ws >> dist.m_chi;
+            return is;
+        }
+#endif
+
+        normal_type m_norm;
+        chi_square_type m_chi;
+    };
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP

Added: sandbox/random/boost/random/weibull_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/random/boost/random/weibull_distribution.hpp	2008-01-14 08:26:55 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,96 @@
+// Copyright Andrew Sutton 2007
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_RANDOM_WEIBULL_DISTRIBUTION_HPP
+#define BOOST_RANDOM_WEIBULL_DISTRIBUTION_HPP
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/assert.hpp>
+
+namespace boost {
+
+    // A Weibull distribution for random number generation. The implementation
+    // is based off the description found in Wikipedia with alpha as lambda and
+    // beta as k.
+    //
+    // Note that the weibull distribution as a special case of the exponential
+    // distribution s.t. Y ~ Weibull(a, b) if X ~ Exponential(b^-a) and Y = X^-a.
+    // This implementation is not built as such.
+    //
+    // Note that MathWorld uses alpha/beta rather than lambda and k.
+    template <typename ResultType = double>
+    class weibull_distribution
+    {
+    public:
+        typedef ResultType result_type;
+
+        explicit weibull_distribution(const result_type& alpha = result_type(1),
+                                      const result_type& beta = result_type(1))
+            : m_alpha(alpha)
+            , m_beta(beta)
+        {
+            BOOST_ASSERT(alpha > result_type(0));
+            BOOST_ASSERT(beta > result_type(0));
+        }
+
+        result_type alpha() const
+        { return m_alpha; }
+
+        result_type beta() const
+        { return m_beta; }
+
+
+        void reset()
+        { }
+
+        // Note that Engine must generate uniformly distributed random numbers
+        // in the range (0, 1).
+        template<class Engine>
+        result_type operator()(Engine& eng)
+        {
+            using std::log;
+            using std::pow;
+
+            // Generate a random number in [0, 1) that is _not_ 0. Taking the
+            // log of 0 is bad!
+            result_type x;
+            for(x = eng(); x == 0; x = eng()) ;
+
+            result_type one = result_type(1);
+            result_type y = m_alpha * (pow(-log(x), one / m_beta));
+
+            return y;
+        }
+
+#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)
+        template<class CharT, class Traits>
+        friend std::basic_ostream<CharT,Traits>&
+        operator<<(std::basic_ostream<CharT,Traits>& os, const weibull_distribution& dist)
+        {
+            os << dist.m_alpha << " " << dist.m_beta;
+            return os;
+        }
+
+        template<class CharT, class Traits>
+        friend std::basic_istream<CharT,Traits>&
+        operator>>(std::basic_istream<CharT,Traits>& is, weibull_distribution& dist)
+        {
+            is >> std::ws >> dist.m_alpha >> std::ws >> dist.m_beta;
+            return is;
+        }
+#endif
+
+    private:
+        result_type m_alpha;
+        result_type m_beta;
+    };
+
+} // namespace boost
+
+#endif
+