#include<iostream> #include<ctime> #include<cmath> using namespace std; class rand_gen{ public: rand_gen(); double uniform_dist(double a, double b); double exponential_dist(double lamda); double powerlaw_dist(double alpha, double lamda); int binomial_dist(int n, double p); int poisson_dist(double lambda); double normal_dist(double sigma, double mu); }; rand_gen::rand_gen(){ srand(time(NULL)); } double rand_gen::uniform_dist(double a, double b){ double r=rand()/(RAND_MAX+1.0); return a+(b-a)*r; } /* To simulate exponential distribution exp(-lambda*x), the inverse method is used. The cumulative distribution function for the exponential distribution is:1-exp(-lambda*x). The inversion function is -log(1-U)/lambda. The simplified form is -log(U)/lambda, where U is a uniform [0,1] random variable. http://cg.scs.carleton.ca/~luc/rnbookindex.html */ double rand_gen::exponential_dist(double lambda){ return(-1*log(uniform_dist(0.,1.))/lambda); } /* Simulate power law with cut-off x^(-alpha)*exp(-lambda*x) To simulate power-law with cutoff , one can generate an exponentially distributed random number using the formula above (as k>0 and integer, so k start at 1) and then accept or reject it with probability p or 1 - p respectively (i.e accept U1 <p or reject U1>p, and U1 is a uniform [0,1] random variable), where p = (x/x_min)^(-alpha) and x_min=1. http://www.santafe.edu/~aaronc/powerlaws/ */ double rand_gen::powerlaw_dist(double alpha, double lamda) { double x; do{ x = exponential_dist(lamda); } while (pow(x,-1*alpha) < uniform_dist(0.,1.)); return (x); } int rand_gen::binomial_dist(int n, double p){ int m=0; for( int i=0;i<n;i++){ if(uniform_dist(0.,1.)<p) m++; } return m; } /* Poisson distribution f(n, lambda) = lambda^n exp(-lambda)/n! Algorithm Poisson random number (Knuth): init: Let L =exp(-lambda), k=0 and p =1 do: k = k+ Generate uniform random number u in [0,1] and let p =p*u while p > L. return k - 1. */ int rand_gen::poisson_dist(double lambda){ double L=exp(-1*lambda); int k; double p=1.; k = 0; do{ k = k + 1; p = p*uniform_dist(0.,1.); } while (p > L); return k - 1; } double rand_gen::normal_dist(double sigma, double mu) { double V1, V2, S; double X; do { double U1 = uniform_dist(0.,1.); double U2 = uniform_dist(0.,1.); V1 = 2 * U1 - 1; V2 = 2 * U2 - 1; S = V1 * V1 + V2 * V2; } while(S >= 1 || S == 0); X = V1 * sqrt(-2 * log(S) / S); return sigma*X+mu; } int main(){ rand_gen myrand_gen; cout<<"uniform distribution between 0 and 4"<<endl; for(int i=1; i<10; i++) { cout <<myrand_gen.uniform_dist(0.,4.)<<endl; } cout<<"exponential distribution exp(-2*x)"<<endl; for(int i=1; i<10; i++) { cout <<myrand_gen.exponential_dist(2.)<<endl; } cout<<"power law distribution x^-2*exp(-2*x)"<<endl; for(int i=1; i<10; i++) { cout <<myrand_gen.powerlaw_dist(2.,2.)<<endl; } cout<<"binomial distribution"<<endl; for(int i=1; i<20; i++) { // flip the coin 4 times, count number of heads cout <<myrand_gen.binomial_dist(4,0.5)<<endl; } cout<<"Poisson distribution"<<endl; for(int i=1; i<20; i++) { cout <<myrand_gen.poisson_dist(4)<<endl; } cout<<"Normal distribution"<<endl; for(int i=1; i<20; i++) { // flip the coin 4 times, count number of heads cout <<myrand_gen.normal_dist(1.,0.)<<endl; } return(0); }
Online computer courses, code, programming tutorial and sidebar information for monitoring Canadian S&P/TSX index. Build friendship and networking. Welcome to visit my blogs often!!! I also have two other sites: YouTube Channel and Google site.
Adsense
Popular Posts
- PHP: add a download as pdf file button in report page
- How to blend adsense inside your post?
- Formatting my post
- PHPWind-- A PHP forum script applcaition in China
- Notepad++ - Add C++ compiler
- PHP connect IBM db2 database in XAMPP
- Datatable export excel wraptext and newline
- phpexcel toggle expand and hide column in EXCEL and summary
- Sweet Alert JS library - beautiful replacement of JavaScript Alert
- ActionScript 3.0 demo: create a falling snow in flash CS6
Saturday, June 26, 2010
random number generator -C++
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment