Generating random numbers of exponential distribution

Question

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double ran_expo(lambda){
    double u;

    u = rand() / (RAND_MAX + 1.0);

    return -log(1- u) / lambda;
}

I am using this (I do not put other parts of the codes here) to generate random numbers of an exponential distribution lambda = 0.05. (Probability density function being lambda * exp(-lambda * x)).

However, I always get very small numbers, such as 0.000041, or something like 1.#INF00 (what is this?).

In fact, for exponential distribution with lambda = 0.05, numbers generated should be quite large, namely most are not far from 30. My results are very weird. In addition, the degree of precision is not satisfactory, only to 10^(-6). I tried long double, but it is still like this.

I am using DEV C++ under windows.

The mean of your distribution is 30, but that doesn't mean most of the values you generate should be near that. The mode is actually 0. — Adam Acosta, Jan 01 '16 at 17:33

score 5 · Accepted Answer · answered Jan 01 '16 at 17:36

5

This is because you did not declare the type for lambda, having corrected that the results are in the range you seek. If undeclared, old compilers will assume it to be int.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

double ran_expo(double lambda){
    double u;
    u = rand() / (RAND_MAX + 1.0);
    return -log(1- u) / lambda;
}

int main(void)
{
    int i;
    srand((unsigned)time(NULL));
    for (i=0; i<20; i++)
        printf("%f\n", ran_expo(0.05));
    return 0;
}

Program output:

answered Jan 01 '16 at 17:36

Weather Vane

33,872
7
36
56

Thanks, I was doubtful at first about your range of random numbers (concerned about `log(0)`) but I see you subtracted from `1`. – Weather Vane Jan 01 '16 at 17:42
There is a glaring fault in my answer. There is a difference between assuming that `int` is the default type, and declaring `int lambda`, but I don't know what it is. In the first case, every result is `-0.000000`. But with the explicit `int` type every result is `1.#INF00`, as you would expect, the result of div 0. (MSVC). – Weather Vane Jan 01 '16 at 18:06

score 0 · Answer 2 · answered Jan 01 '16 at 17:33

I threw your code verbatim into code chef, it looks like it works fine. Compiler is gcc c++ 4.9.2 . This code prints out 36.6751:

#include <iostream>
#include <string>
#include <cstdlib>
#include <cmath>

using namespace std;

double ran_expo(double lambda){
    double u;

    u = rand() / (RAND_MAX + 1.0);

    return -log(1- u) / lambda;
}


int main() {
    cout << ran_expo(0.05);
    return 0;
}

Works the same way in C, too (again, gcc 4.9.2):

#include <stdlib.h>
#include <math.h>

double ran_expo(double lambda){
    double u;

    u = rand() / (RAND_MAX + 1.0);

    return -log(1- u) / lambda;
}


int main() {
    printf("%f\n", ran_expo(0.05));
    return 0;
}

Must be your compiler/libraries?

Generating random numbers of exponential distribution

2 Answers2

Linked