Re: Bell shaped random number generation

From: Robert Dodier (robert_dodier_at_yahoo.com)
Date: 09/28/04 

Date: 27 Sep 2004 17:20:56 -0700

harish_recian@yahoo.com (Harish) wrote:

> I need to write a computer program in which I need to get random
> numbers in a range, that follow an approximation of a bell curve. I
> tell the program that the average is say 20 and also say that most of
> values lie between 15 & 25.
> The range of the values is between 0,100. I have been suggested that I
> would need to look at beta distribution to generate such numbers but I
> could not figure out how to do it.

Yes, using a beta distribution makes sense to me. A beta distribution
is defined on the interval (0,1) so just multiply by 100 to get
something in the range (0,100).

There are different ways to sample from a beta distribution. One
relatively efficient algorithm is the rejection method. Another
method, less efficient but simpler to program, is to sample two
numbers X1 and X2 from chi square distributions and then X1/(X1+X2)
has a beta distribution.

I've appended some code below. I also have an implementation of the
rejection method; email me for it. If I were going to need a lot
(say 1,000,000's) of samples, I'd forego Java entirely and use the
function rbeta in R (http://www.r-project.org/).

For what it's worth,
Robert Dodier

// by Robert Dodier, distributed under terms of GNU GPL
// usage: java BetaChiSquare A B n
// prints n samples from a beta distribution w/ parameters A and B
// mean value is A/(A+B), variance is AB/((A+B+1) (A+B)^2)
// disclaimer: this program is slow for large A or B

import java.util.*;

public class BetaChiSquare
{
    double A, B;
    Random rand;

    public BetaChiSquare( double A, double B )
    {
        this.A = A;
        this.B = B;
        System.err.println( " A: "+(int)A+" B: "+(int)B );
        rand = new Random();
    }

    public double next()
    {
        int i;
        double z;

        double chisq_a = 0;
        for ( i = 0; i < 2*A; i++ )
        {
            z = rand.nextGaussian();
            chisq_a += z*z;
        }

        double chisq_b = 0;
        for ( i = 0; i < 2*B; i++ )
        {
            z = rand.nextGaussian();
            chisq_b += z*z;
        }

        return chisq_a/(chisq_a+chisq_b);
    }

    public static void main( String[] args )
    {
        int i, n = Integer.parseInt( args[2] );
        BetaChiSquare bcs = new BetaChiSquare(
Double.parseDouble(args[0]), Double.parseDouble(args[1]) );

        for ( i = 0; i < n; i++ )
            System.out.println( bcs.next() );
    }
}

Loading