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

#include "stdafx.h" 
 
#define _CRT_RAND_S 
#include  
#define _USE_MATH_DEFINES 
#include  
 
#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(&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=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