This is the normal desintegration decay distribution.
The problem you face is not how to determine that this is the right distribution but to produce a random generator that has this property. Normally you work taking as a base another random number generator, and try to apply a function that from that generator of know properties you get numbers distributed in the new fashion.
For your approach I've used several times the standar random generator (a flat number generator between [0..1), as source) if I apply a function to it, let's say (just a guess) the log(x)
function you will get a graph approx like this:
|-> x
|-> x
|-> x
|-> x
|-> x
|-> x
+====================================================================
and it seems to match approximately the accumulation of points desired (indeed, it is the right approach)
The integral to this function (that gives you the probability of having X <= a)
is given by the function represented below:
|-> X
|-> X
|-> x
|-> X
|-> X
|-> X
+====================================================================
0
and it's inverse (a logarithm function is given below, mirrored on the XY bisectrix):
| |
| X |
| |
| X |
|X |
+=============
0 1
which is -log(1-x)
(but -log(x)
will suffice, as x is random between 0 and 1).
double y = -A * log(unit_random());
in which A is some constant specific to produce a larger or smaller average. y
will have the expected response, always the unit_random()
is a flat distribution with equally probability distribution in the unit interval.
Below is a complete sample of a little program to generate points with the requested distribution. The function
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#include <stdio.h>
#include "gr.h"
double
geometric_random(void)
{
double n = random();
return -log(fabs(n) / INT_MAX);
} /* geometric_random */
produces the desired random values with an average value of 1.0
. To get an average of M
, just multiply the values returned by it by M
.
A complete example is published in github, with a montecarlo test program to test for distribution properties. To execute it, run:
$ make
$ gr -n 10000000 | mc -n 100 -b 20 >mc.out
and you'll see how the ratios of the counters for the different subintervals are constant over the whole range of values (well, not at low frequecies, as expected)
Output is:
n: 10000000
sum_x: 10007209.0488715283572674
sum_x2: 20023758.8835431151092052
avg_x: 1.0007209048871528
sdev_x: 1.0004667205707121
min_x: 0.0000002696179218
max_x: 17.2248827198513297
below_A: 0
[0.0000000000000000, 0.2000000000000000]1811635
[0.2000000000000000, 0.4000000000000000]1484598: 0.8194796413184775
[0.4000000000000000, 0.6000000000000001]1213219: 0.8172037144061894
[0.6000000000000001, 0.8000000000000000]994937: 0.8200802987754066
[0.8000000000000000, 1.0000000000000000]813669: 0.8178095698521615
[1.0000000000000000, 1.2000000000000000]666035: 0.8185576690275775
[1.2000000000000000, 1.3999999999999999]545997: 0.8197722341918968
[1.3999999999999999, 1.5999999999999999]447841: 0.8202261184585263
[1.5999999999999999, 1.7999999999999998]365854: 0.8169283294740768
[1.7999999999999998, 1.9999999999999998]300525: 0.8214342333280489
[1.9999999999999998, 2.1999999999999997]246141: 0.8190366857998502
[2.1999999999999997, 2.3999999999999999]201303: 0.8178361183224250
[2.3999999999999999, 2.6000000000000001]164134: 0.8153579430013462
[2.6000000000000001, 2.8000000000000003]134768: 0.8210852108642938
[2.8000000000000003, 3.0000000000000004]110564: 0.8204024694289446
[3.0000000000000004, 3.2000000000000006] 90139: 0.8152653666654607
[3.2000000000000006, 3.4000000000000008] 74252: 0.8237499861325287
[3.4000000000000008, 3.6000000000000010] 60746: 0.8181059096051285
[3.6000000000000010, 3.8000000000000012] 49701: 0.8181773285483818
[3.8000000000000012, 4.0000000000000009] 40559: 0.8160600390334198
[4.0000000000000009, 4.2000000000000011] 33344: 0.8221109987918834
[4.2000000000000011, 4.4000000000000012] 27077: 0.8120501439539347
[4.4000000000000012, 4.6000000000000014] 22338: 0.8249806108505373
[4.6000000000000014, 4.8000000000000016] 18370: 0.8223654758707136
[4.8000000000000016, 5.0000000000000018] 15016: 0.8174197060424605
...
std::geometric_distribution
or if you need real numbersstd::exponential_distribution
– Mgetz Sep 19 '19 at 18:00