V151: Random Numbers According to Any Function

Author(s): F. James	Library: MATHLIB
Submitter:	Submitted: 27.11.1984
Language: Fortran	Revised:

FUNRAN generates random numbers distributed according to any (one-dimensional) distribution f(x). The distribution is supplied by the user in the form of a FUNCTION subprogram. If the distribution is known as a histogram only, HISRAN (V150) should be used instead.

Structure:

SUBROUTINE subprograms
User Entry Names: FUNRAN, FUNPRE
Internal Entry Names: FUNZER
Files Referenced: Printer
External References: GAUSS, RNDM, user-supplied FUNCTION subprogram
COMMON Block Names and Lengths: /FUNINT/ 1

Usage:

CALL FUNPRE(F,FSPACE,XLOW,XHIGH) (once for each function)

CALL FUNRAN(FSPACE,XRAN) (for each random number)

F

(REAL) A name of a FUNCTION subprogram declared EXTERNAL in the calling program. This subprogram must calculate the (non-negative) density function , for all X in the interval .

FSPACE

(REAL) One-dimensional array of length 100.

XLOW

(REAL) Lower limit of the requested interval.

XHIGH

(REAL) Upper limit of the requested interval.

XRAN

(REAL) A random number returned by FUNRAN.

A call to FUNPRE calculates the percentiles of F between XLOW and stores them into the array FSPACE.

Method:

In FUNPRE, the percentiles are calculated using a combination of trapezoidal and Gaussian integration to a rather high accuracy, which is printed out by FUNPRE. If the desired accuracy is not obtained, an warning is printed in addition.

Subroutine FUNRAN finds the desired random number by calling RNDM (V104) and doing a 4-point interpolation on FSPACE to transform the uniform random number to the distribution specified. This method produces quite accurately distributed numbers even when the function F is badly skew or spiked as long as the width of a spike is not less than 1/1000 of the total range.

Error handling:

An error message is printed

- if the integral of the user-supplied function F is zero or negative,

- if ,

- if somewhere between XLOW and XHIGH.

Notes:

Some additional information which may be of use is contained in

    COMMON / FUNINT/ FINT

After a call to FUNPRE, FINT contains the integral of F from XLOW to XHIGH.

After a call to FUNRAN, FINT contains the integral of F from XLOW to XRAN, divided by the total integral to XHIGH (i.e., it will be a number uniformly distributed between zero and one).