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