It's been awhile so I figured: it is time for a random post. Seriously, this is a post about random numbers.
If you code in Javascript, then Math.random() is probably already familiar to you. It does a nice job of generating (pretty) uniform variables from 0 to 1. A little multiplying and rounding and you can get integers from 0 to x, or -x to x, or whatever your needs.
But what if you want NON-uniform randoms? For statistics freaks: what if you need normal distribution? or Gamma? or Beta? Exponential? Chi-squared?
Well you are in luck! I found myself in need of some common distributions for a little side game I'm making, and went the extra mile to read some thesis papers about how to generate these useful distributions.
Rather than bore you with how it all works, here is the code for a RANDOM object you can insert into your code. I'm still working on Betas, they are bit trickier.
RANDOM.norm(n)
This returns an (approximately) normal bell curve centered on mean 0.5. The parameter "n" can be between 2 and 5 (default 2). Higher values of "n" makes the bell curve more narrow. The value returned is between 0 and 1.
RANDOM.inorm(n)
This
returns an inverse bell curve centered on mean 0.5. The
parameter "n" can be between 2 and 5 (default 2). Higher n makes output more likely to be near 0 or 1.
RANDOM.exp(n)
A nice easy exponential variable from 0 to 1. Increasing "n" makes the graph skewed towards 0. Default value for "n" is 2.
RANDOM.invexp(n)
This is just RANDOM.exp(n) times -1.
RANDOM.gamma(alpha,beta)
This is a conglomeration of a few generators. For alpha > 1 and < 0, it uses an algorithm developed by Kundu and Gupta. No beta is used for this algorithm.
For alpha >= 1 it uses the algorithms proposed by Martino and Luengo. It generates very quick results due to a low rejection rate. Default values are 1 and 1.
I could not have completed this without those great algorithms! If you use this code, please leave the credits they deserve.
TIP: RANDOM.gamma(2,2) returns a very nice gamma distribution. Be aware: the values returned may be higher than 1 (especially for high alpha and beta values). Alpha values below 1 resemble exponential distributions.
RANDOM.igamma(alpha,beta)
Simply RANDOM.gamma() times -1.
RANDOM.chi2(k)
A commonly used distribution resembling gamma distributions (I use this one a lot). Again, returned values may be higher than 1, especially for high "k".
RANDOM.coinFlip(weight)
A handy weighted coin flip that returns 0 or 1. Default weight is 2 (50/50) and MUST be > 1. Weights over 2 favor a positive outcome (i.e. weight of 3 will have 2:3 odds of returning 1).
The code:
As usual, the code is free to use under the GNU. Please give credit though, and I'd love to hear how this has been useful for anyone! Also, feel free to play around with the fiddle! (The code is adjusted slightly for display).
I'll post a follow-up if I ever get around to adding Beta randoms too!
var RANDOM = { //returns random between 0 and 1, normal distribution centered on 0.5
"norm": function(n) {
if (!(n > 2 && n <= 5)) n = 2;
var nrand = 0;
n = Math.floor(n);
for (var i = 1;i<=n*2;i++){
nrand += Math.random();
}
return nrand/(2*n);
},
"inorm": function(n) { //returns random between 0 and 1
if (!(n > 2 && n <= 5)) n = 2;
var nrand = 0;
n = Math.floor(n);
for (var i = 1;i<=n*2;i++){
nrand += Math.random();
}
return ((1 - Math.abs((nrand-n) / n))*(Math.abs(nrand-n)/(nrand-n)) + 1)/2;
},
"exp": function(n) { //returns greater than 0
if (!(n > 2 && n <= 5)) n = 2;
var nrand = Math.random();
for (var i = 2;i<=n;i++){
nrand *= Math.random();
}
return 2*nrand;
},
"invexp": function(n) { //returns less than 0
return -RANDOM.exp(n);
},
"gamma3": function(alpha) { //Kundu and Gupta algorithm 3 http://home.iitk.ac.in/~kundu/paper120.pdf
if (!alpha || Math.abs(alpha) > 1) alpha = 1; //alpha between 0 and 1
var d = 1.0334 - (0.0766*Math.exp(2.2942*alpha));
var a = Math.pow(2,alpha)*Math.pow(1-Math.exp(-d/2),alpha);
var b = alpha*Math.pow(d,alpha-1)*Math.exp(-d);
var c = a + b;
var U = Math.random();
var X = (U <= a/(a+b)) ? -2*Math.log(1-(Math.pow(c*U,1/alpha)/2)) : -Math.log(c*(1-U)/(alpha*Math.pow(d,alpha-1)));
var V = Math.random();
if (X <= d) {
var mess = (Math.pow(X,alpha-1)*Math.exp(-X/2))/(Math.pow(2,alpha-1)*Math.pow(1-Math.exp(-X/2),alpha-1));
if (V <= mess) return X;
else return this.gamma3(alpha);
} else { //X > d
if (V <= Math.pow(d/X,1-alpha)) return X;
else return this.gamma3(alpha);
}
//output is > 0 and possibly > 1
},
"gamma": function(alpha,beta) { //Martino and Luengo http://arxiv.org/pdf/1304.3800.pdf luca@tsc.uc3m.es luengod@ieee.org
if (!alpha || alpha <= 0) alpha = 1; //alpha >= 1 if negative or 0
if (alpha > 0 && alpha < 1) return this.gamma3(alpha); // use different algorithm
if (!beta || beta <= 0) beta = 1; //beta > 0
var alphap = Math.floor(alpha);
var X = Math.random();
for (var i=2;i<=alphap;i++){
X *= Math.random();
}
var betap = (alpha < 2) ? beta/alpha : beta*(alphap-1)/(alpha-1);
X = -Math.log(X)/betap;
var Kp = (alpha < 2) ? Math.exp(1-alpha)*Math.pow(alpha/beta,alpha-1) : Math.exp(alphap-alpha)*Math.pow((alpha-1)/beta,alpha-alphap);
//then accept with prob p(X)/pi(X)
if (alphap >= 2) {
if (Kp*Math.pow(X,alphap-1)*Math.exp(-betap*X) >= Math.pow(X,alpha-1)*Math.exp(-beta*X)) return X/alpha;
else return this.gamma(alpha,beta);
}
else if (alphap < 2) {
if (Kp*Math.exp(-betap*X) >= Math.pow(X,alpha-1)*Math.exp(-beta*X)) return X/alpha;
else return this.gamma(alpha,beta);
}
},
"igamma": function(alpha,beta) { // returns less than 0
return -RANDOM.gamma(alpha,beta);
},
"chi2": function(k) { // returns greater than 0
var nrand = RANDOM.norm(2);
nrand = nrand*nrand;
if (!k || k <= 1) return nrand;
for (var i=2;i<=k;i++){
var krand = RANDOM.norm(2);
krand = krand*krand;
nrand += krand;
}
return nrand;
},
"coinFlip": function(weight){
if (!weight || weight < 1) weight = 2;
if (Math.random() > 1/weight) return 1;
else return 0;
}
};
//Copyright 2016 Trevor Gast: codeandcompose.com
//RANDOM, non-uniform random generators (5 May 2016)
//GNU General Public License
Dice gif courtesy of
http://bestanimations.com/Games/Dice/Dice.html.png)
0 comments: