steve_bank
Diabetic retinopathy and poor eyesight. Typos ...
For a distribution such as normal the practical problem is generating random samples for arbitrary means and standard deviations.
That is the quantile function. For a normal distribution the inverse error function is part of the normal quantile. Random normally distrusted numbers have a range of uses for example Monte Carlo methods. Same with the exponential distribution.
Using a factorial table improves speed.
Part of my small math library I have been working on.
I tried numeral integration for the erf() and it was not as acurate as the math tools. The approximation closely matches libraries.
That is the quantile function. For a normal distribution the inverse error function is part of the normal quantile. Random normally distrusted numbers have a range of uses for example Monte Carlo methods. Same with the exponential distribution.
Using a factorial table improves speed.
Part of my small math library I have been working on.
I tried numeral integration for the erf() and it was not as acurate as the math tools. The approximation closely matches libraries.
Code:
class RAND_DIST{
private:
//URAND random number algorithm
//URAND A Universal Random Number Generator by Malcom and Moler
unsigned long long randll = 1; //random number
unsigned long long m = 0x7fffffff;
//unsigned long long k1 = 2*floor(m*(.5 - sqrt(3)/6.)) + 1;
//unsigned long long k2 = 8*ceil(m*atan(1.)/8.) + 5;
unsigned long long k1 = 2.*(m*(.5 - sqrt(3)/6.)) + 1;
unsigned long long k2 = 8.*(m*(atan(1.)/8.) + 5.);
public:
void init(unsigned long long seed);
int rand_exp(long long n,double *y,double u);
int rand_nor(long long n,double *y,double u,double sigma);
int uni_float(int n,double *y,int lo,int hi,int ndp);
int uni_int(int n,int *y,int lo,int hi);
double erf(double x);
double ierf(double x);
};//RAND_DIST
void RAND_DIST::init(unsigned long long seed){
randll = seed;
}//init()
double RAND_DIST::erf(double x){
double factorial[23] = {
1.,
1.,
2.,
6.,
24.,
120.,
720.,
5040.,
40320.,
362880.,
3628800.,
39916800.,
479001600.,
6227020800.,
87178291200.,
1307674368000.,
20922789888000.,
355687428096000.,
6402373705728000.,
121645100408832000.,
2432902008176640000.,
51090942171709440000.,
1124000727777607680000.};
int i,n = 20;
double s = 0,k;
for(i=0;i<23;i++){
k = (1./(2*i+1))*pow(-1,i)/factorial[i];
s += k * pow(x,2*i+1);
}
return s*2./sqrt(_PI);
}//erf()
double RAND_DIST::ierf(double x){
//inverse error function
if(x==0)return 0;
//if(x <= -1 || x >= 1.)return 999;
int k,m,nc = 1000;
double einv = 0.,q = sqrt(_PI)*x/2.;
#include "coef.txt" //c[]
//double c[nc];
//for(k=0;k<nc;k++)c[k] = 0;
//c[0] = 1;
//for(k=0;k<nc;k++){
//for(m=0;m<k;m++)
// c[k]+= (c[m]*c[k-1-m])/((m+1)*(2*m + 1));
// }
for(k=0;k<nc;k++)einv += ((c[k])/(2*k+1))*pow(q,(2*k+1));
return einv;
}//ierf()
int RAND_DIST::rand_exp(long long n,double *y,double u){
// array of random numbers from an exonentual distribution
// p represents random probabilities >0 <1
if(u <= 0 || n <= 0)return -1;
double p;
long long i;
for(i=0;i<n;i++){
randll = randll*k2 + k1;
p = (double)(randll%m+1)/(double)(m+2);
// quantile - probability to variable for a given mean
y[i] = -1.* log(1. - p) * u;
}
return 0;
}//rand_exp()
int RAND_DIST::rand_nor(long long n,double *y,double u,double sigma){
// array of random numbers from a normal distribution
// URAND random number algorithm
// p random probabilities >0 <1
if(n <= 0)return - 1;
double q,p,k1;
int i,k,nc = 0;
q = sqrt(2)*sigma;
for(i=0;i<n;i++){
randll = randll*k2 + k1;
// p >0 <1
//p = (double(rand())+.01)/double(RAND_MAX + 1);
p = (double)(randll%m+1e-6)/(double)(m+1e-3);
p = 2.*p-1.;
y[i] = u + ierf(p)*q; //normal qantile
}
return 0;
}//rand_nor()
int RAND_DIST::uni_float(int n,double *y,int lo,int hi,int ndp){
//uniform distribution floats > lo< hi + 1
if(lo >= hi || n <= 0)return -1;
int xint,i,nnums = (hi-lo),dp = pow(10,ndp);
double xdec;
for(i=0;i<n;i++){
randll = randll*k2 + k1;
xdec = (double)(randll%m)/(double)(m+1);
if(ndp>0){
xdec = xdec*dp;
xdec = int(xdec);
xdec = xdec/dp;
}
randll = randll*k2 + k1;
xint = (randll%m)%nnums + lo;
y[i] = xint + xdec;
}
return 0;
}//uni_float()
int RAND_DIST::uni_int(int n,int *y,int lo,int hi){
//uniform distribution ints lo -> hi-1
if(lo >= hi || n <= 0)return -1;
int i,nnums = hi-lo;
for(i=0;i<n;i++){
randll = randll*k2 + k1;
if(hi < 2)y[i] = (randll%m)%nnums;
y[i] = (randll%m)%nnums + lo;
}
return 0;
}//uni_int(
1) smooth frequency with boxes
You seem to admit that the function returns variates correctly, what difference do the internals matter? When p is uniform on (0, 1), it and (1-p) provide exactly equivalent distributions.
