📄 random.h
字号:
// -*- c++ -*-///////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2001 Oh-Wook Kwon, all rights reserved. ohwook@yahoo.com
//
// Easy Matrix Template Library
//
// This Easy Matrix Template Library is provided "as is" without any express
// or implied warranty of any kind with respect to this software.
// In particular the authors shall not be liable for any direct,
// indirect, special, incidental or consequential damages arising
// in any way from use of the software.
//
// Everyone is granted permission to copy, modify and redistribute this
// Easy Matrix Template Library, provided:
// 1. All copies contain this copyright notice.
// 2. All modified copies shall carry a notice stating who
// made the last modification and the date of such modification.
// 3. No charge is made for this software or works derived from it.
// This clause shall not be construed as constraining other software
// distributed on the same medium as this software, nor is a
// distribution fee considered a charge.
//
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// Filename: Random.h
//
// Note:
// 1. Distributions other than uniform and Gaussian were not tested.
// You need to test the resulting distribution before use.
// In particular, the Gamma distribution has a poor pdf when a<=1.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef _RANDOM_H_#define _RANDOM_H_
// Undefine rand and srand to avoid name collisions
// Use RANDF and SRAND instead.
#ifdef USE_LOCAL_RAND
#ifdef rand
#undef rand
#endif
#ifdef srand
#undef srand
#endif
#else
#include <stdlib.h>
#endif
// RANDF() returns a uniform random number [0,1).
#if defined(USE_LOCAL_RAND)
#define RANDF() UniformRandom()
#define SRAND(x) SRand(x)
#define SRANDN(x) SRandn(x)
double Rand();
void SRand(unsigned int seed);
#elif defined(__unix)// Prototype for C Library functions deand48 and srand48.double drand48();void srand48(long);
#define RANDF() drand48()#define SRAND(x) srand48(x)
#define SRANDN(x) SRandn(x)#elif defined(WIN32)// If not unix use ANSI C defaults.#define __rand() ((double)rand()/RAND_MAX)#define RANDF() __rand()
#define SRAND(x) srand(x)
#define SRANDN(x) SRandn(x)
#else
#pragma message ( "Your system does not have random number generators " __FILE__ )
#pragma message ( "Define __unix or WIN32 " )
!!error!! // to stop compiling#endif
double UniformRandom();
double GaussianRandom();
void SRandn(unsigned int seed);
template <typename T> class UniformDist;
template <typename T> class GaussianDist;
template <typename T> class MiscellaneousDist;
///////////////////////////////////////////////////////////////////////////////
// Exported routines
///////////////////////////////////////////////////////////////////////////////
// continuous distributions
#ifdef USE_LOCAL_RAND
inline double UniformRandom();
inline void SRand(unsigned int seed);
#endif
inline double GaussianRandom();
inline void SRandn(unsigned int seed);
inline double ExponentialRandom(double lambda);
inline double TRandom(double n); // Student's t-distribution with n degree of freedom
inline double GammaRandom(double a, double b);
inline double BetaRandom(double a, double b);
inline double FisherRandom(double m, double n); // m degrees of freedom in the numerator and n degrees of freedom in the denominator
inline double ChiSquareRandom(double n); // n degree of freedom
inline double CauchyRandom();
// discrete distributions
inline int BernoulliRandom(double p);
inline int BinomialRandom(int n,double p);
inline int PossionRandom(double lambda);
///////////////////////////////////////////////////////////////////////////////
// Implementation of random number generation
///////////////////////////////////////////////////////////////////////////////
#ifdef USE_LOCAL_RAND
struct Int32__ {
unsigned short high;
unsigned short low;
};
///////////////////////////////////////////////////////////////////////////////
// There must exist only one uniform random number generator in the whole system.
// I had to implement a singleton object using a static member function and
// a static object defined within the static member function.
///////////////////////////////////////////////////////////////////////////////
template <typename T>
class UniformDist {
public:
UniformDist() {}
virtual ~UniformDist() {}
static UniformDist<double>* GetUniformDist();
static void SRand_(unsigned int seed);
static double Rand_();
private:
unsigned int rand_seed;
Int32__ rand_value;
};
template <typename T> UniformDist<double>* UniformDist<T>::GetUniformDist(){
static UniformDist<double> uniformDist;
return &uniformDist;
}
template <typename T> void UniformDist<T>::SRand_(unsigned int seed){
UniformDist* p=UniformDist::GetUniformDist();
p->rand_seed=seed;
p->rand_value.high=seed/65536;
p->rand_value.low=seed%65536;
}
///////////////////////////////////////////////////////////////////////////////
// Uniform distribution [0,1)
// x = (69069*x+1)%pow(2,32);
// u = x/pow(2,32)
// Reference: C.P. Robert, The Bayesian Choice, 2nd ed., Springer, 2001.
///////////////////////////////////////////////////////////////////////////////
template <typename T> double UniformDist<T>::Rand_()
{
UniformDist* p=UniformDist::GetUniformDist();
static Int32__ c={1,3533}; // 69069
static double rand_max=pow((double)2.0,(double)32);
unsigned int x1=c.low*p->rand_value.low+1;
unsigned int x2=c.low*p->rand_value.high;
unsigned int x3=c.high*p->rand_value.low;
unsigned int x4=c.high*p->rand_value.high;
p->rand_value.low=x1 & 0xffff;
unsigned int x1h=(unsigned int)(0xffff & ((x1&0xffff0000)>>16));
unsigned int x2l=(unsigned int)(x2&0xffff);
unsigned int x3l=(unsigned int)(x3&0xffff);
p->rand_value.high=(unsigned int)(x1h+x2l+x3l);
double u=(unsigned int)(65536*p->rand_value.high+p->rand_value.low)/rand_max;
return u;
}
#endif
///////////////////////////////////////////////////////////////////////////////
// Gaussian distribution
///////////////////////////////////////////////////////////////////////////////
template <typename T>
class GaussianDist {
public:
GaussianDist() { gaussSaved=0; }
virtual ~GaussianDist() {}
static GaussianDist<double>* GetGaussianDist();
static void SRandn_(unsigned int seed);
static double Randn_();
private:
int gaussSaved; // GaussDeviate generates numbers in pairs
double gaussSave; // 2nd of pair is remembered here
#ifdef USE_LOCAL_RAND
unsigned int rand_seed;
Int32__ rand_value;
#endif
};
template <typename T> GaussianDist<double>* GaussianDist<T>::GetGaussianDist(){
static GaussianDist<double> gaussianDist;
return &gaussianDist;
}
template <typename T> void GaussianDist<T>::SRandn_(unsigned int seed){
GaussianDist* p=GaussianDist<double>::GetGaussianDist();
#ifdef USE_LOCAL_RAND
p->rand_seed=seed;
p->rand_value.high=seed/65536;
p->rand_value.low=seed%65536;
#endif
p->gaussSaved=0;
p->gaussSave=0;
}
// random number with a N(0,1) distribution
template <typename T> double GaussianDist<T>::Randn_()
{
GaussianDist* p=GaussianDist<double>::GetGaussianDist();
double fac,r,v1,v2,x;
if (p->gaussSaved) {
x = p->gaussSave; p->gaussSaved = 0;
}
else {
do {
#ifdef USE_LOCAL_RAND
v1 = 2.0*UniformRandom() - 1.0;
v2 = 2.0*UniformRandom() - 1.0;
#else
v1 = 2.0*RANDF() - 1.0;
v2 = 2.0*RANDF() - 1.0;
#endif
r = v1*v1 + v2*v2;
}
while (r>=1.0);
fac = sqrt(-2.0*log(r)/r);
p->gaussSaved = 1;
p->gaussSave = v1*fac;
x = v2*fac;
}
return x;
}
///////////////////////////////////////////////////////////////////////////////
// Other distributions
///////////////////////////////////////////////////////////////////////////////
template <typename T>
class MiscellaneousDist {
public:
// continuous distributions
static double ExponentialRandom_(double lambda=1);
static double TRandom_(double n);
static double GammaRandom_(double a, double b=1);
static double BetaRandom_(double a, double b);
static double FisherRandom_(double m, double n);
static double ChiSquareRandom_(double lambda);
static double CauchyRandom_();
// discrete distributions
static int BernoulliRandom_(double p);
static int BinomialRandom_(int n,double p);
static int NegativeBinomialRandom_(int n,double p);
static int PossionRandom_(double lambda);
static int HyperGeometricalbRandom_(int N,int n,double p);
};
template <typename T> double MiscellaneousDist<T>::ExponentialRandom_(double lambda){
double u;
while(1){
u=UniformRandom();
if(u!=0) break;
}
return -log(u)/lambda;
}
// n degree of freedom
template <typename T> double MiscellaneousDist<T>::TRandom_(double n){
while(1){
double u1=UniformRandom();
double u2=UniformRandom();
double x,v;
if(u1 < 0.5){
x=1/(4*u1-1);
v=u2/(x*x);
}
else{
x=4*u1-3;
v=u2;
}
if( v < 1-(fabs(x)/2) || v < pow(1+(x*x/n), -(n+1)/2.0) )
return x;
}
}
template <typename T> double MiscellaneousDist<T>::GammaRandom_(double a, double b){
assert(a>0 && b>0);
double u1,u2;
if(a > 50){
return GaussianRandom()*sqrt(a)/b + (a/b);
}
else if(a > 1){
double c1=a-1, c2=(a-(1/(6*a)))/c1, c3=2/c1, c4=1+c3, c5=1/sqrt(a);
while(1){
while(1){
u1=UniformRandom();
u2=UniformRandom();
if(a > 2.5){
u1 = u2 + c5 * (1-1.86*u1);
}
if(u1 > 0 && u1 < 1) break;
}
double w=c2*u2/u1;
if(c3*u1+w+1/w <= c4 || c3*log(u1)-log(w)+w <= 1)
return (b*c1*w);
}
}
else{
double e=exp(1);
while(1){
u1=UniformRandom();
u2=UniformRandom();
double x,y;
if(u1 > e/(e+a)){
x=-log((a+e)*(1-u1)/(a*e));
y=pow(x,a-1);
}
else{
x=pow((a+e)*u1*e,1/a);
y=exp(-a);
}
if(u2<y)
return (b*x);
}
}
}
template <typename T> double MiscellaneousDist<T>::BetaRandom_(double a, double b){
double y1=GammaRandom(a,1);
double y2=GammaRandom(b,1);
return y1/(y1+y2);
}
// m degrees of freedom in the numerator and n degrees of freedom in the denominator
template <typename T> double MiscellaneousDist<T>::FisherRandom_(double m, double n){
double xnom=ChiSquareRandom(m)/m;
double xdenom;
while(1){
xdenom=ChiSquareRandom(n)/n;
if(xdenom > xnom/mtl_numeric_limits<double>::max()) break;
}
return xnom/xdenom;
}
// n degree of freedom
template <typename T> double MiscellaneousDist<T>::ChiSquareRandom_(double n){
return GammaRandom(n/2.0,2.0);
}
template <typename T> double MiscellaneousDist<T>::CauchyRandom_(){
static double pi=4*atan(1);
//f(x)=1/(pi*(1+x^2))
return tan(pi*(UniformRandom()-0.5));
}
template <typename T> int MiscellaneousDist<T>::BernoulliRandom_(double p){
if(UniformRandom()<=p) return 1;
else return 0;
}
template <typename T> int MiscellaneousDist<T>::BinomialRandom_(int n,double p){
int K=n;
if(n <= 30){
int count = 0;
for(int i=1;i<=n;i++){
double u=UniformRandom();
if(u < p) count++;
}
return count;
}
else{
int i, k=n, x=0;
double l=p;
while(1){
i=(int)(1+k*l);
double v=BetaRandom(i,k+1-i);
if(l>v){
l=l/v;
k=i-1;
}
else{
x=x+i;
l=(l-v)/(1-v);
k=k-i;
}
if(k<=K) break;
}
for(i=1;i<=k;i++){
double u=UniformRandom();
if(u<p) x++;
}
return x;
}
}
template <typename T> int MiscellaneousDist<T>::PossionRandom_(double lambda){
if(lambda<30){
double p=1;
int N=0;
double c=exp(-lambda);
while(1){
double u=UniformRandom();
p=p*u;
N=N+1;
if(p<c) break;
}
return N-1;
}
else{
static double pi=4*atan(1);
double c=0.767-3.36/lambda;
double beta=pi*pow((2*lambda),-0.5);
double alpha=beta*lambda;
double k=log(c)-lambda-log(beta);
double x;
while(1){
while(1){
double u1=UniformRandom();
x=(alpha-log((1-u1)/u1))/beta;
if(x>-0.5) break;
}
double u2=UniformRandom();
int N=(int)(x+0.5);
if( alpha-beta*x+log(u2/(1+exp(2*(alpha-beta*x)))) <= k+N*log(lambda)-log(N) )
return N;
}
}
}
////////////////////////////////////////////////////////////////////////////
// Implementation of the exported routines
////////////////////////////////////////////////////////////////////////////
#ifdef USE_LOCAL_RAND
inline double UniformRandom(){
return UniformDist<double>::Rand_();
}
inline void SRand(unsigned int seed){
UniformDist<double>::SRand_(seed);
}
#endif
inline double GaussianRandom(){
return GaussianDist<double>::Randn_();
}
inline void SRandn(unsigned int seed){
#ifdef USE_LOCAL_RAND
GaussianDist<double>::SRandn_(seed);
#endif
}
inline double ExponentialRandom(double lambda){
return MiscellaneousDist<double>::ExponentialRandom_(lambda);
}
// n degree of freedom
inline double TRandom(double n){
return MiscellaneousDist<double>::TRandom_(n);
}
inline double GammaRandom(double a, double b){
return MiscellaneousDist<double>::GammaRandom_(a,b);
}
inline double BetaRandom(double a, double b){
return MiscellaneousDist<double>::BetaRandom_(a,b);
}
// m degrees of freedom in the numerator and n degrees of freedom in the denominator
inline double FisherRandom(double m, double n){
return MiscellaneousDist<double>::FisherRandom_(m,n);
}
// n degree of freedom
inline double ChiSquareRandom(double n){
return MiscellaneousDist<double>::ChiSquareRandom_(n);
}
inline double CauchyRandom(){
return MiscellaneousDist<double>::CauchyRandom_();
}
inline int BernoulliRandom(double p){
return MiscellaneousDist<double>::BernoulliRandom_(p);
}
inline int BinomialRandom(int n,double p){
return MiscellaneousDist<double>::BinomialRandom_(n,p);
}
inline int PossionRandom(double lambda){
return MiscellaneousDist<double>::PossionRandom_(lambda);
}
#endif⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -