举例：FUNLUX, FUNLXP---------Random Numbers According to Any Function
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).
本帖转自52MC论坛