In simulations, you frequently need to generate random number following a distributions other than uniform distributions. In previous exercises, you have already created a rather simple way to draw a random number from the binominal distribution and geometric distribution (remember the coin flipping functions in the homeworks?).

Here I briefly discuss about approaches to create non-uniform random numbers from uniform PRNG.

In reality, I generally use pre-made functions from a library called GNU Scientific Library (GSL).

e.g., To create a random number *y* from exponential
distribution, draw a unifrom PRN *z* [0,1), and use the following
transformation.

*g(x)*: envelope function, whose shape is similar to the actual distribution (*f(x)*). This function should lie*above**f(x)*, and the inverse of can be calculated easily.- Draw a first uniform PRN, and transform it with the inverse function of
to create a candidate
*j*. - Draw a second uniform PRN between 0 and
*g(j)*. - If the 2nd PRN is less than f(j), accept it as a PRN drawn from
a distribution
*f(x)*. If not, reject it, and draw another candidate until you find an acceptable one.

- Example
#include <stdio.h> #include <gsl/gsl_rng.h> #include <gsl/gsl_randist.h> const gsl_rng *gBaseRand; /* global rand number generator */ int main (void) { unsigned long randSeed; int i, k, n = 10; double lambda = 3.0; /* specifying to use Mersenne twister MT-19937 as the uniform PRNG */ gBaseRand = gsl_rng_alloc(gsl_rng_mt19937); srand(time(NULL)); /* initialization for rand() */ randSeed = rand(); /* returns a non-negative integer */ gsl_rng_set (gBaseRand, randSeed); /* seed the PRNG */ /* print n random variates chosen from the poisson distribution with mean parameter lambda */ for (i = 0; i < n; i++) { k = gsl_ran_poisson (gBaseRand, lambda); printf (" %u", k); } printf ("\n"); gsl_rng_free(gBaseRand); return 0; }

- Explanation of key elements:
- You need to use the two header files
`const gsl_rng *gBaseRand;`

Making a

*special*global variable (I'm calling it`gBaseRand`), which stores what kind of uniform random number generator you want to use.`gBaseRand = gsl_rng_alloc(gsl_rng_mt19937);`

We are setting this variable to use

`gsl_rng_mt19937`

.

There are many other PRNGs which you can specify (e.g., you can use the familiar`drand48()`if you specify`gsl_rnd_rand48`

).`gsl_rng_set (gBaseRand, randSeed);`

Seeding the specified PRNG (

`gBaseRand`) with a seed (`randSeed`).`k = gsl_ran_poisson (gBaseRand, lambda);`

Use

`gBaseRand`as the base uniform PRNG, and draw a random number from a poisson distribution (whose parameter is specified as`lambda`).`gsl_rng_free(gBaseRand);`

Free the memory allocated for the PRNG

- Compilation (You can copy
`/home/progClass/src/prng/prng.c`.)gcc prng.c -lgsl -lgslcblas -lm

`-l`means*libraryName***link**with the library.- Linking means that extract the functions used in your source code from the library and attach them to your object code.
- To use the random number generators, you need to link with 3 libraries: gsl, gslcblas (basic linear algebra subprograms), and m (math).

- Other useful random number functions
The first argument for all the functions are the variable containing the type of uniform PRNG (here I'm using

`gBaseRand`).- Discrete distribution

Retuns`unsigned int`.

`gsl_ran_binomial (gBaseRand, double p, unsigned int n)`

`gsl_ran_geometric (gBaseRand, double p)`

`gsl_ran_negative_binomial (gBaseRand, double p, double n)`

`gsl_ran_poisson (gBaseRand, double lambda)`

Returns

`unsigned long int``gsl_rng_uniform_int(gBaseRand, unisigned long int n)`

*n-1*inclusive. - Continuous distributions

All returns`double`.`gsl_ran_gaussian (gBaseRand, double sigma)`

Normal distribution with mean of 0 `gsl_ran_exponential (gBaseRand, double theta)`

Exponential distribution `gsl_ran_flat (gBaseRand, double a, double b)`

Uniform (flat) distribution, *a <= x < b*, [a, b)`gsl_rng_uniform(gBaseRand)`

Uniform, 0 <= *x*< 1, [0, 1)`gsl_rng_uniform_pos(gBaseRand)`

Uniform, 0 < *x*< 1, (0, 1)

- Discrete distribution

There are many more distributions in GSL, and you can use them in the same way as the example above.

- Create a program to print out the random numbers drawn from a
whatever distribution implemented in GSL function. Take a
look at the
manual
if you want to try other kinds of dstributions not listed above.
Print out many random numbers (say 10000, each line contains one
random number). Capture the output to a file (e.g
`out1.txt`). Then use R to plot the distribution.`$ R`

> dat <- read.table("out1.txt")

> hist(dat[,1])

> mean(dat[,1])

> var(dat[,1])

> q()Compare the observed mean and variance with the expectation if you are using one of the distribution we covered in the class.

Try to use several different distributions. Also change the parameters of distributions, and observe how the shape changes.

- Model of minority cytotype exclusion in polyploidy.

When there is a population with mixed ploidy, inter-ploidy cross can create infertile offspring. This create a competition between ploidy levels. Implement a stochastic model simulating this process, and analyze the outcome. For examples, you can analyze how the initial frequencies influence the final outcome, how time till exclusion by one ploidy is influenced by different parameter values, whether coexistence of two ploidy levels can be possible under certain parameter combinations. Or you can start from a diploid population and introduce a tetraploid individual, and see when the tetraploid individual can invade into the population.*t*tetraploids and*d*diploids in a population (*N = t + d*), and N stays constant.- Each plant makes only 1 flower (like tulips). Number of
ovules per flower is normally distributed with mean
*o*and standard deviation*w*. - Proportion
*s*of the ovules are self-fertilized (binominal distribution). The rest are fertilized by other plants (outcrossed). - Assume that the fitness of outcrossed and selfed seeds are the same (no inbreeding depression).
- Number of pollinator visits per flower is the Poisson distribution
with parameter
*v*. - Each pollinator visit results in deposition of 100 pollen
grains from the flower visitation immediately before this flower.
Assuming that the polinator visit plants randomly, the probability
that the previous visit was diploid is
*d/N*. - If
*x*ovules are available for outcrossing and if the flower received*a*diploid pollen and*b*tetraploid pollen, the number of ovule fertilized by diploid pollen is binominally distributed (number of trials =*x*and probability of success is*a/(a+b)*. - For a simplicity, let's assume that triploid seeds get aborted after fertilization (not biologically realistic). In other words, the ovules used for interploidy crosses are wasted.
- After one generation, there will be
*l*diploid seeds, and*m*tetraploid seeds. To create the next generation, you can use the binominal distribution.