C programming:Gamma pseudo random distribution

Nagara · Oct 7, 2015

Hello All.🙂
I need help with some help with generation of Gamma pseudo random distribution using on pure C. I found this. But it's too rough. And they consider not all cases in C-code snippet I guess. Besides, I can not understand what kind of Pseudo random generator they use.
I found how to generate gauss variable here. But I guess tha there is some thing more. Some kind of transform should be applied to that

Thank a lot.()🙂

sm625 · Oct 7, 2015

Seems like a bunch of nonsense to me. I just call rand() and seed it with

Code:

srand( (unsigned)time( NULL) );

If I want something more random than that, I seed the random number generator using time, then call rand() repeatedly to generate an array of random seed values. I'll use rand() to select one of the seed values from the array as follows

Code:

srand( seed_array[rand()&0xF] );

I'll also throw random bitwise operations in there

Code:

value = rand();
if(value & (rand() & 0xF))
value <<= rand()&0xF;
else
value >>= rand()&0xF;

At some point this has to produce something that is about as close to true random as you can possibly get.

K7SN · Oct 7, 2015

Gamma Distributions are easy if you think of shape and scale as Mu(mean) and Sigma(Deviation). Generating them correctly is not so easy. Below is a C# way I've done it.

Code:

[FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]public[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]static[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] [] get_Gamma_Distribution([/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] n, [/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] shape, [/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] scale)[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]{ [/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]  int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] i, k, m1 = 0;[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]  double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] a, x, x1;[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]  double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] s1,s2,y;[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]  double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2][] data = [/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]new[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2][n];[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]  k = ([/FONT][/SIZE][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2]) shape;[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]  a = shape - ([/FONT][/SIZE][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2]) k;[/SIZE][/FONT]
 
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#008000][FONT=Consolas][SIZE=2][COLOR=#008000][FONT=Consolas][SIZE=2][COLOR=#008000]  // routine for gamma deviates [/COLOR][/SIZE][/FONT]
[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT]
[FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]  do[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT]
[SIZE=2][FONT=Consolas][SIZE=2][FONT=Consolas]  {[/FONT][/SIZE]
[SIZE=2][FONT=Consolas]       s1 = [/FONT][/SIZE][/FONT][/SIZE][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Math[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Pow([/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Utility[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Rand().NextDouble(), 1.0 / a);[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]       s2 = [/FONT][/SIZE][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Math[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Pow([/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Utility[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Rand().NextDouble(),(1.0 / (1.0-a)));[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]       if[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] ((s1 + s2) < 1.0) [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]       {[/FONT][/SIZE]
[SIZE=2][FONT=Consolas]            x1 = 1.0;[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]            if[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (k > 0) [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]            {[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]                 for[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (i = 0; i < k; i++)[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]                 x1 *= [/FONT][/SIZE][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Utility[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Rand().NextDouble();[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]            }[/FONT][/SIZE]
[SIZE=2][FONT=Consolas]            y = - scale * [/FONT][/SIZE][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Math[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Log(x1);[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]            x = - scale * (s1 * [/FONT][/SIZE][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Math[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Pow(s1 + s2, -1)) * [/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Math[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Log([/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Utility[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Rand().NextDouble());[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]            if[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (x + y > 0.000000001)  // Adjust to your precision needs[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]                if[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (x + y < 1e+31) // Adjust to your precision needs[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]                    data[m1++] = x + y;[/FONT][/SIZE]
[SIZE=2][FONT=Consolas]       }[/FONT][/SIZE]
[SIZE=2][FONT=Consolas]   } [/FONT][/SIZE][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]while[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2](m1 < n);[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]   return[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] data;[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]}[/FONT][/SIZE]
[/SIZE][/FONT]

In C# I have utility routines (shown below) to explain the called random functions in the code above: Below is my code, it is a lot easier in C

Code:

[FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]public[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]class[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Utility[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2]
[SIZE=2][FONT=Consolas]{[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]   static[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Random[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] rand = [/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]null[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2];[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]   public[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]static[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Random[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] Rand() [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]   {   [/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]        if[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2](rand == [/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]null[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2]) [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]           rand = [/FONT][/SIZE][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]new[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Random[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2]([/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]unchecked[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2](([/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2])[/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]DateTime[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2].Now.Ticks)); [/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]         return[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] rand;[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]   }[/FONT][/SIZE]
 
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]   private[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]static[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Random[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] random = [/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]new[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af][FONT=Consolas][SIZE=2][COLOR=#2b91af]Random[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2](); [/SIZE][/FONT]
[/SIZE][/FONT]
[FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]   public[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]static[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]double[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] NextDouble() [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]   { [/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]        lock[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (random) [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]        { [/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]             return[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] random.NextDouble(); [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]        } [/FONT][/SIZE]
[SIZE=2][FONT=Consolas]    } [/FONT][/SIZE]
 
[/SIZE][/FONT]
[FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]    public[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]static[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] Next([/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] val) [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]    { [/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]        lock[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (random) [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]        { [/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]              return[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] random.Next(val); [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]        } [/FONT][/SIZE]
[SIZE=2][FONT=Consolas]     } [/FONT][/SIZE]
[/SIZE][/FONT]
[FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]     public[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]static[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] Next([/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] min, [/SIZE][/FONT][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] max) [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]     { [/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]         lock[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (random) [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]         { [/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]              return[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] random.Next(min, max); [/SIZE][/FONT]
[SIZE=2][FONT=Consolas]      } [/FONT][/SIZE]
[SIZE=2][FONT=Consolas]  } [/FONT][/SIZE]
[SIZE=2][FONT=Consolas]} [/FONT][/SIZE]
[/SIZE][/FONT]

If you want a truly random gamma distribution, mot one for a given shape and scale, call your function with random values for shape and scale.

Now tell me why you are messing around with generation gamma distributions. I may have done that myself.

Essence_of_War · Oct 7, 2015

If you grab the latest gnu science library from here:

ftp://ftp.gnu.org/gnu/gsl/gsl-latest.tar.gz

You can check out their source code for the gamma distribution inside folder randist -> gamma.c

A relevant snippet w/ reference:

Code:

/* New version based on Marsaglia and Tsang, "A Simple Method for
 * generating gamma variables", ACM Transactions on Mathematical
 * Software, Vol 26, No 3 (2000), p363-372.
 *
 * Implemented by J.D.Lamb@btinternet.com, minor modifications for GSL
 * by Brian Gough
 */

double
gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b)
{
  return gsl_ran_gamma (r, a, b);
}

double
gsl_ran_gamma (const gsl_rng * r, const double a, const double b)
{
  /* assume a > 0 */

  if (a < 1)
    {
      double u = gsl_rng_uniform_pos (r);
      return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
    }

  {
    double x, v, u;
    double d = a - 1.0 / 3.0;
    double c = (1.0 / 3.0) / sqrt (d);

    while (1)
      {
        do
          {
            x = gsl_ran_gaussian_ziggurat (r, 1.0);
            v = 1.0 + c * x;
          }
        while (v <= 0);

        v = v * v * v;
        u = gsl_rng_uniform_pos (r);

        if (u < 1 - 0.0331 * x * x * x * x) 
          break;

        if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
          break;
      }
    
    return b * d * v;
  }
}

K7SN · Oct 7, 2015

Essence_of_War said:
If you grab the latest gnu science library from here:

ftp://ftp.gnu.org/gnu/gsl/gsl-latest.tar.gz

You can check out their source code for the gamma distribution inside folder randist -> gamma....

Been working in C# too long, the gnu library is probably where I got my C# code,, 15 years ago.

repoman0 · Oct 7, 2015

sm625 said:
Seems like a bunch of nonsense to me. I just call rand() and seed it with

Code:

srand( (unsigned)time( NULL) );

If I want something more random than that, I seed the random number generator using time, then call rand() repeatedly to generate an array of random seed values. I'll use rand() to select one of the seed values from the array as follows

Code:

srand( seed_array[rand()&0xF] );

I'll also throw random bitwise operations in there

Code:

value = rand(); if(value & (rand() & 0xF)) value <<= rand()&0xF; else value >>= rand()&0xF;

At some point this has to produce something that is about as close to true random as you can possibly get.

You're talking about trying to generate uniform random numbers, not a Gamma distribution, and in any case this way would only have marginally more entropy than just seeding with time and calling rand. Besides, rand as implemented in C uses an LCG as far as I know, far from being cryptographically secure ... i.e., very very far from about as close to true random as one can get.

McFox · Oct 7, 2015

If you need exactly Gamma probability distribution exactly in C you can try this.
It looks like that it's right implementation of marsaglia example for especially [FONT=MathJax_Main]1[/FONT][FONT=MathJax_Main]Γ[/FONT][FONT=MathJax_Main]([/FONT][FONT=MathJax_Math]a[/FONT][FONT=MathJax_Main])[/FONT][FONT=MathJax_Math]b[/FONT][FONT=MathJax_Math]a[/FONT][FONT=MathJax_Math]x[/FONT][FONT=MathJax_Math]a[/FONT][FONT=MathJax_Main]−[/FONT][FONT=MathJax_Main]1[/FONT][FONT=MathJax_Math]e[/FONT][FONT=MathJax_Main]−[/FONT][FONT=MathJax_Math]x[/FONT][FONT=MathJax_Main]/[/FONT][FONT=MathJax_Math]b[/FONT].

I guess that's the main problem it's to test it😉

http://www.hongliangjie.com/2012/12/19/how-to-generate-gamma-random-variables/

Nagara · Oct 8, 2015

Thank you all, guys!🙂
McFox, yeah pure C, good example, thanks. That's realy good question how will I know that this is realy gamma distribution and not something else:hmm:

Merad · Oct 8, 2015

sm625 said:
Seems like a bunch of nonsense to me. I just call rand() and seed it with

Code:

srand( (unsigned)time( NULL) );

If I want something more random than that, I seed the random number generator using time, then call rand() repeatedly to generate an array of random seed values. I'll use rand() to select one of the seed values from the array as follows

Code:

srand( seed_array[rand()&0xF] );

I'll also throw random bitwise operations in there

Code:

value = rand(); if(value & (rand() & 0xF)) value <<= rand()&0xF; else value >>= rand()&0xF;

At some point this has to produce something that is about as close to true random as you can possibly get.

Ugh, no. Your last sentence isn't even remotely true.

First of all, rand() is often a very poor quality PRNG. IIRC, the C standard places basically no requirements on it, so you really have no idea what you're getting without inspecting the implementation. Also, RAND_MAX is often quite small, limiting the range of numbers that will be produced.

Second, changing the seed (or even how you produce the seed) has no effect whatsoever on the quality of the PRNG. Each seed value produces a distinct sequence of numbers, and will always product that sequence when the same seed is set. Changing the seed changes the sequence of numbers, not the distribution.

Finally, doing bitwise operations isn't magic that makes things better. In fact, in the code snippet you posted, you're removing bits of randomness, and actually making the distribution worse.

Generating random numbers isn't a simple topic, and it's a very mathy topic. If you just need some numbers, rand() or whatever your language provides should be fine. If you're doing something where the quality of the numbers matters, or for whatever reason you thing the basic language PRNG isn't good enough, you should switch to a better PRNG like Mersenne Twister. If you're doing something where the quality really matters, you need to be using a cryptographically secure generator like /dev/random or a hardware RNG.

Crusty · Oct 8, 2015

Merad said:
Ugh, no. Your last sentence isn't even remotely true.

First of all, rand() is often a very poor quality PRNG. IIRC, the C standard places basically no requirements on it, so you really have no idea what you're getting without inspecting the implementation. Also, RAND_MAX is often quite small, limiting the range of numbers that will be produced.

Second, changing the seed (or even how you produce the seed) has no effect whatsoever on the quality of the PRNG. Each seed value produces a distinct sequence of numbers, and will always product that sequence when the same seed is set. Changing the seed changes the sequence of numbers, not the distribution.

Finally, doing bitwise operations isn't magic that makes things better. In fact, in the code snippet you posted, you're removing bits of randomness, and actually making the distribution worse.

Generating random numbers isn't a simple topic, and it's a very mathy topic. If you just need some numbers, rand() or whatever your language provides should be fine. If you're doing something where the quality of the numbers matters, or for whatever reason you thing the basic language PRNG isn't good enough, you should switch to a better PRNG like Mersenne Twister. If you're doing something where the quality really matters, you need to be using a cryptographically secure generator like /dev/random or a hardware RNG.

Not quite true re: /dev/random. When using Linux there is almost never any reason to choose /dev/random over /dev/urandom.

http://www.2uo.de/myths-about-urandom/

Essence_of_War · Oct 8, 2015

Not quite true re: /dev/random. When using Linux there is almost never any reason to choose /dev/random over /dev/urandom.

Yeah, this.

/dev/urandom uses the CSPRNG as well.

K7SN · Oct 8, 2015

Nagara said:
Thank you all, guys!🙂
McFox, yeah pure C, good example, thanks. That's realy good question how will I know that this is realy gamma distribution and not something else:hmm:

You could download a program I wrote if you can run a windows system?
http://www2.epa.gov/land-research/proucl-software

C programming:Gamma pseudo random distribution

Nagara

Junior Member

sm625

Diamond Member

K7SN

Senior member

Essence_of_War

Platinum Member

K7SN

Senior member

repoman0

Diamond Member

McFox

Junior Member

Nagara

Junior Member

Merad

Platinum Member

Crusty

Lifer

Essence_of_War

Platinum Member

K7SN

Senior member

TRENDING THREADS