V152: Random Numbers According to Any Function

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

FUNLUX 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: FUNLUX, FUNLXP
Internal Entry Names: FUNPCT, FUNLZ
Files Referenced: Printer
External References: RADAPT, RANLUX, user-supplied FUNCTION subprogram
COMMON Block Names and Lengths: /FUNINT/ 1

Usage:

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

CALL FUNLUX(FSPACE,XRAN,LEN) (for each vector of random numbers)

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 200.

XLOW

(REAL) Lower limit of the requested interval.

XHIGH

(REAL) Upper limit of the requested interval.

XRAN

(REAL) A vector of random numbers returned by FUNRAN.

LEN

(INTEGER) Length of the vector XRAN.

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

Method:

In FUNLXP, the 100 percentiles of the integral of are calculated using a combination of trapezoidal and Gaussian integration to a rather high accuracy, which is printed out by FUNLXP. Then both the left-hand and right-hand 2 percentiles are expanded to 50 percentiles each in order to cater for functions with long tails. If the desired accuracy is not obtained, a warning is printed in addition.

Subroutine FUNLUX finds the desired random number by calling RANLUX (V115) 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 FUNLXP, FINT contains the integral of F from XLOW to XHIGH.

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