How do I get random numbers that follow standard distribution?

[sinan]

I need to generate random numbers; but they have to follow the standard distribution. The usual random number generation routines do not do this.

In other words, if I have a range of 1 to 100, I want the numbers around 50 to occur frequently, and the numbers about 1 or about 100 to occur with a very low frequency. With the usual random number generation routines, all numbers have the same frequency.

There must be a simple algorithm that create standard distribution numbers using the usual generator. (i.e. An algorithm that converts such numbers)
Does anybody know any function that does this? Can anybody write a simple algorithm for this purpose?

Please do not tell me any non-standard modules because I have trouble installing them.

TIA,
Sinan

Originally posted as a Categorized Question.

[tilly]

If you don't need a perfect normal distribution then you can get a very good approximation by choosing several random numbers and adding them together. You can then multiply then add by appropriate constants to get the desired average and variance.

For most purposes, adding a dozen is good enough.

If you want better than that, try installing Math::Random.

DISCLAIMER: I am one dissertation short of a PhD in math. (It is hand-written, needed an income so I never cleaned it up and typed it.) Therefore for all intents and purposes, I am a mathematician. :-)

[lhoward]

How about something like this:

sub f{
  my $t=0;
  $t+=rand() for(1..5);
  return int($t*20+1);
}

You can vary the number of iterations (1..5) and scaling factor (20) based on the exact distribution curve you want.

There are several other ways to approach this problem. If efficiency were of paramount importance, I would try using statistical method to derive a function that would aproximate the curve of how in frequency you would map numbers from 0 to 0.9999 (the range of perl's rand function) into the range you want (1 to 100).

[grackle]

The question posted asks for a non-module solution, and an excellent one was provided by lhoward and tilly. Here is a solution using a CPAN module.

I looked for Math::Random on CPAN and couldn't find it. There are modules for random numbers, such as Crypt::TrulyRandom, but I believe they will only give you better uniformly distributed random variables (correct me if I'm wrong). For a normal distribution, you can use Math::CDF (Cumulative Distribution Functions). The relevant function is Math::CDF::qnorm(), which will give you what you want when fed a random number x, 0<x<1, with a uniform distribution (i.e., what you get from rand()). (It supposedly returns a value for the inputs x=1 and x=0, but I don't know what those values could be, you'll have to try it).

So the function you want is:

$myrand = qnorm(rand());

(or replace rand() with your function of choice).

One caveat: I don't know how the Math::CDF module calculates these values. If it provides more precision than you need, it could be a big waste of processing time. In that case lhoward and tilly's solution, summing several random variables, works fine. It would be nice to have bounds for the errors, but I can't help you there. Calculate, experiment, or find a reference.... Perhaps this would be a good place to post the results.

[tilly]

You misunderstood one thing. The point of Math::TrulyRandom is not to have a more accurate distribution, it is to have a random number that is not guessable. Normal random numbers are guessable based on the time your program runs.

If the random numbers are used for any cryptographic purpose this is an important consideration.

[grackle]

The following sentence was truncated due to author error:

The relevant function is Math::CDF::qnorm(), which will give you what you want when fed a random number x, 0 < x < 1.

My apologies.

All exact science is dominated by the idea of approximation.
-- Bertrand Russell

[Roy Johnson]

The Perl Cookbook provides this code:

sub gaussian_rand {
    my ($u1, $u2);  # uniformly distributed random numbers
    my $w;          # variance, then a weight
    my ($g1, $g2);  # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w = $u1*$u1 + $u2*$u2;
    } while ( $w >= 1 );

    $w = sqrt( (-2 * log($w))  / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;
    # return both if wanted, else just one
    return wantarray ? ($g1, $g2) : $g1;
}

and notes:

The gaussian_rand function [generates] two numbers with a mean of 0 and a standard deviation of 1 (i.e., a Gaussian distribution). To generate numbers with a different mean and standard deviation, multiply the output of gaussian_rand by the new standard deviation, and then add the new mean

[Roy Johnson]

This page claims that the while condition should be } while ($w>=1 or $w==0);

[hv]

$w == 0 will occur only when $u1 == $u2 == 1/2, which will happen rarely enough that I can't see it skewing the results noticeably - even if you only have 32-bit random numbers, you'll hit this case only 1.27 in 2^64 times.

(That 1.27 is based on my rusty attempts to calculate how many times we loop until we get $w < 1. I think the probability of a success is arcsin(1)/2, which is about 0.78.)

Hugo

[polettix]

I do support the previous answer: the code should include the test for $w being zero, otherwise you risk a division by zero in the very next line.

The proposed solution roots in Numerical Recipes in C, chapter 7, section 2; probably it' better to update the Perl Cookbook and the answer in Q&A as well.

[Anonymous Monk]

sub norm{
    my $v1,$v2,$r;
    if( defined $v2 ){
        $v1 = $v2;
        undef $v2;
    }else{
        do{
            $v1=rand(2)-1;
            $v2=rand(2)-1;
            $r = $v1*$v1+$v2*$v2;
        }while( $r >= 1 || $r == 0 );
        $r = sqrt(-2*(log $r)/$r);
        $v1 *= $r;
        $v2 *= $r;
   }
   return $v1;
}

[jimbojones]

Just to let others know: this method appears to be from "Numerical Recipes in C", which is now available, full-form, online. See c7-2.pdf.
this also features recipes for other distributions, plus lots of other good stuff ....
- j

[Roy Johnson]

This is a buggy translation from Knuth. $v2 will never be defined; it should have been declared outside the sub.

[sinan]

I finally did something that solves the problem partially.

The two following functions do that. The first function, rands will return a "standard" randomized integer. It has two input values: The first input variables denote the range of the output. (i.e. 0 to input) The second one will denote the precision: Higher it is, closer the result will be to the normal distribution. 5 is more than enough for many purposes.

The second function (random) is a little bit more tricky. The first input is the range, just like in the first function. The second input is a string which can either be "low" or "high". The third input is a limit. If the second input is "low", then the third input denotes the lowest possible number. If it is high, then it denotes the highest possible number. All the results that are outside either the lowest or the highest frequency will be ignored and "re-rolled".

Here is the snippet:


sub rands{
    my $t=0;
    $t = 0;
    $t+= ( rand()* ($_[0]/$_[1]) ) for(1..$_[1]);
    return int($t+1);
}

sub random{
    my ($no,$type,$limit) = @_;
    
    $t = rands($no,5);

    if($type eq "low"){
        while ($t<$limit){
            $t = rands($no,5)
        }
    }
    elsif($type eq "high"){
        while ($t>$limit){
            $t = rands($no,5)
        }
    }
}

I was not able to "skew" my results, so I used this _harsh_ method of cutting off instead.

One final note: You can set the second input variable of rands to a low number (2,3) to get a more homogenous distribution, i.e. a distribution with "high variance".

Sinan

[jlistf]

i think you're going to have to use statistics. i don't believe Perl has anything built-in to do it, though there might be a module on CPAN. Perl's rand($n) will give you numbers (more or less) uniformly distributed between 0 and 1. its up to you to convert this to a normalized distribution.

[chas]

By "standard distribution", I think you mean "standard Normal distribution". One way to get a random sample from such a distribution is as follows: Get a random sample x_1,...x_n from a uniform distribution on the unit interval. (this is easily done using the usual random number generators you refer to.) Then for each x_i, compute z_i satisfying
x_i=(1/sqrt(2 pi))\int_{-\infty}^z_i e^{-t^2/2}dt. (I.e. x_i=F(z_i) where F is the distribution function for a N(0,1).) Then z_1,...z_n will be a random sample from a N(0,1). You have to write a bit of code to solve for z_i, of course. I've done similar things in c, but not perl so I can't display perl code to achieve this.
The previous suggestion to add the results of several uniformly distributed random numbers is interesting; one needs independence to make use of the Central Limit Theorem, so I'm not sure exactly how to set that up, but it should be possible.
The method I described first is pretty standard.
chas

[Anonymous Monk]

3

Originally posted as a Categorized Answer.