Re: Bell shaped random number generation

From: Alan Miller (amiller_at_bigpond.net.au)
Date: 09/28/04 

Date: Tue, 28 Sep 2004 01:45:00 GMT

"Harish" <harish_recian@yahoo.com> wrote in message
news:a3ef9093.0409252303.74881b55@posting.google.com...
> Hello,
> 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.
>
> I would be thankful to you if anyonce can provide any further input in
> this regard
>
> Thanks
> Harish

And here is another generator of random deviates from the beta distribution.

FUNCTION random_beta(aa, bb, first) RESULT(fn_val)

! Adapted from Fortran 77 code from the book:
! Dagpunar, J. 'Principles of random variate generation'
! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9

! FUNCTION GENERATES A RANDOM VARIATE IN [0,1]
! FROM A BETA DISTRIBUTION WITH DENSITY
! PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).
! USING CHENG'S LOG LOGISTIC METHOD.

! AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
! BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)

REAL, INTENT(IN) :: aa, bb
LOGICAL, INTENT(IN) :: first
REAL :: fn_val

! Local variables
REAL, PARAMETER :: aln4 = 1.3862944
REAL :: a, b, g, r, s, x, y, z
REAL, SAVE :: d, f, h, t, c
LOGICAL, SAVE :: swap

IF (aa <= zero .OR. bb <= zero) THEN
  WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)'
  STOP
END IF

IF (first) THEN ! Initialization, if necessary
  a = aa
  b = bb
  swap = b > a
  IF (swap) THEN
    g = b
    b = a
    a = g
  END IF
  d = a/b
  f = a+b
  IF (b > one) THEN
    h = SQRT((two*a*b - f)/(f - two))
    t = one
  ELSE
    h = b
    t = one/(one + (a/(vlarge*b))**b)
  END IF
  c = a+h
END IF

DO
  CALL RANDOM_NUMBER(r)
  CALL RANDOM_NUMBER(x)
  s = r*r*x
  IF (r < vsmall .OR. s <= zero) CYCLE
  IF (r < t) THEN
    x = LOG(r/(one - r))/h
    y = d*EXP(x)
    z = c*x + f*LOG((one + d)/(one + y)) - aln4
    IF (s - one > z) THEN
      IF (s - s*z > one) CYCLE
      IF (LOG(s) > z) CYCLE
    END IF
    fn_val = y/(one + y)
  ELSE
    IF (4.0*s > (one + one/d)**f) CYCLE
    fn_val = one
  END IF
  EXIT
END DO

IF (swap) fn_val = one - fn_val
RETURN
END FUNCTION random_beta 

--
Alan Miller
Retired
Formerly with CSIRO,
Division of Mathematics & Statistics
http://www.users.bigpond.net.au/amiller
Loading