www.pudn.com > NoiseEstimate.rar > random_variate.cpp


#include "stdafx.h"

#define _CRT_RAND_S
#include <stdlib.h>
#define _USE_MATH_DEFINES
#include <math.h>

#include "random_variate.h"


double gen_uniform(double min, double max)
{
//return (min + (double) rand() / (double) RAND_MAX * (max - min));
unsigned int number;
errno_t err;
err = rand_s(&amt;number);
return min+(double)number/(double)UINT_MAX*(max-min);
}

double gen_normal(double mu, double sigma)
{
/* The ratio-of-uniforms method */
double b, a_pos, a_neg;
b=1; a_pos=sqrt(2/M_E); a_neg=-a_pos;
double u, v;
double x;
bool accept;
do
{
u=gen_uniform(0, b); v=gen_uniform(a_neg, a_pos);
x=v/u;
if (x*x<=4-4*u)
{
accept=true;
} else if (x*x>=4/u-4)
{
accept=false;
} else
{
accept=(x*x<=-4*log(u));
}
} while(!accept);

return mu + sigma * x;
}

double gen_exponential(double lambda)
{
return -1/lambda * log(gen_uniform(0, 1));
}

double gen_rayleigh(double sigma)
{
return sigma * sqrt(-2 * log(gen_uniform(0, 1)));
}

double gen_gamma(double alpha, double lambda)
{
double u, v;
double x, y;
double b, c;
double w, z;
bool accept;
double t;

if (alpha > 1.0) {
/* Best's rejection algorithm XG for gamma random variates (Best, 1978) */
b = alpha - 1;
c = 3*alpha-0.75;
do {
u = gen_uniform();
v = gen_uniform();

w = u*(1-u);
y = sqrt(c/w)*(u-0.5);
x = b+y;
if (x>=0) {
z = 64*w*w*w*v*v;
accept=(z <= 1-2*y*y/x);
if (!accept) {
accept = (log(z) <= 2 * (b*log(x/b)-y));
}
}
} while(!accept);
return lambda*x;
} else if (alpha == 1.0) {
return gen_exponential(1/lambda);
} else if (alpha < 1.0) {
/* Algorithm RGS for gamma variates (Best, 1983) */
t=0.07+0.75*sqrt(1-alpha);
b=1+alpha*exp(-t)/t;
c=1/alpha;
do
{
u=gen_uniform(0, 1);
w=gen_uniform(0, 1);
v=b*u;
if (v<=1)
{
x=t*pow(v, c);
accept=(w<=(2-x)/(2+x));
if (!accept)
{
accept=(w<=exp(-x));
}
}
else
{
x=-log(c*t*(b-v));
y=x/t;
accept=(w*(alpha+y-alpha*y)<=1);
if (!accept)
{
accept=(w<=pow(y, alpha-1));
}
}
} while(!accept);
return lambda * x;
}
return -1;
}

double gen_beta(double alpha, double beta)
{
/* Johnk's beta generator */
double u, v;
double x, y;
do
{
u=gen_uniform(0, 1);
v=gen_uniform(0, 1);
x=pow(u, 1/alpha);
y=pow(v, 1/beta);
} while(x+y>1);
return x/(x+y);
}

double Uniform(double x1, double x2)
{
return gen_uniform(0,1) * (x2 - x1) + x1;
}
//配极法产生正态分布, 一次产生两个独立变量, 只取一个即可
//Knuth: Art of Computer Programming, Vol2, Page 107
static double Normal_01_Peiji()
{
double v1,v2,s;
double x1,x2;
do {
v1=gen_uniform(0,1)*2.0-1.0;
v2=gen_uniform(0,1)*2.0-1.0;
s= v1*v1+v2*v2;
} while(s>=1.0);
if (s==0)
x1=x2=0.0;
else{
x1= v1 * sqrt(-2.0 * log(s)/s);
x2= v2 * sqrt(-2.0 * log(s)/s);
}
return x1;
}
//比例法产生正态分布
//Knuth: Art of Computer Programming, Vol2, Page 114
static double Normal_01_Bili()
{
double u,v,x,x2;
const double c1=sqrt(8.0/M_E);
const double c2=exp(0.25)*4;
const double c3=exp(-1.35)*4;
while(true){
do{
u=gen_uniform(0,1);
}while(u==0.0);
v=gen_uniform(0,1);
x=c1*(v-0.5)/u;
x2=x*x;
if (x2<=5.0-c2*u)
break;
if (x2>=c3/u+1.4)
continue;
if (x2<=-4.8*log(u))
break;
}
return x;
}

//正态分布, 目前使用比例方法
double Normal(double avg, double sigma)
{
return Normal_01_Bili() * sigma + avg;
}

//Possion分布
int Poisson(double lamda)
{
static double s_lamda=1.0;
static double expl= exp(-s_lamda);
if (lamda!=s_lamda){
s_lamda= lamda;
expl= exp(-s_lamda);
}
double s=1.0;
int n=0;
do {
s*=gen_uniform(0,1);
++n;
} while(s>expl);
return n-1;
}
//! 三角分布, 密度函数为三角形, 左下x坐标为a, 右下x坐标为b, 且假设a<0<b;
// 高度为 2/(a-b).
// 分布函数为 d(x)= ((1/a)*x*x-2*x+a)/(a-b);when x<=0; ((1/b)*x*x-2*x+a)/(a-b); when x>=0
// d(a)=0; d(0)=a/(a-b); d(b)=1;
double Triangle(double a, double b)
{
double y=gen_uniform(0,1), x=0;
if (y<a/(a-b)){
x=a*(1-sqrt(y*(1-b/a)));
}
else{
x=b*(1-sqrt((1-y)*(1-a/b)));
}
return x;
}