Generates random samples using rejection sampling (W. van den Berg).
||What to plot (
||Size of each random sample; no default, must be set|
||Calculation defining the probability density function f(x) to sample; no default, must be set|
||Data structure used inside
||Lower bound of the region in which f(x) is non-negligible; default -10|
||Upper bound of the region in which f(x) is non-negligible; default 10|
||Calculation defining the probability density function g(x) used to generate the sample; default
||Calculation to sample from the probability density g(x) used to generate the sample (note,
||Multiplier M used in the definition of the envelope M × g(x) that must always be greater than f(x); default 10|
||Number of random samples to take in each sampling step; default
||Saves each random sample|
||Seed to use for the random numbers used to generate each random sample; default 0|
GREJECTIONSAMPLE generates random samples from a probability density function, using rejection sampling. The density function f(x), which need not integrate to one (e.g. a posterior in a Bayesian analysis), is specified as a Genstat expression using the
PRDENSITY option. The
X option specifies the identifier of the data structure used to represent the x-coordinate of the density in the calculation.
The method operates by generating a random sample from a probability density g(x), selected so that
f(x) ≤ M × g(x)
where M is a suitably chosen multiplier. The density g(x) can be defined using the
PRENVELOPE option, as a Genstat expression which uses the same identifier
X as the
PRDENSITY expression. If
PRENVELOPE is not set,
GREJECTIONSAMPLE uses the probability density function of a t-distribution with 60 degrees of freedom i.e.
PRENVELOPE = !e( PRT(X; 60) )
The multiplier M is specified by the
MULTIPLIER option; default 10.
If you set
PRENVELOPE, you must also set the
GRENVELOPE option to define an expression to take a random value z from the distribution g(x).
GREJECTIONSAMPLE also takes a random value from the uniform distribution U[0,1], and accepts z if the uniform random number u is less than or equal to
f(z) / (M × g(z))
The process works more efficiently if several values are sampled at once from g(x) and U[0,1].
GREJECTIONSAMPLE can decide the number automatically; or you can define it yourself using the
NTRIES option. If you are defining
GRENVELOPE you should set
NTRIES to the identifier of a Genstat scalar, and use this in the expression defined by
GRENVELOPE. You set the scalar to a missing value if you still want
GREJECTIONSAMPLE to decide the number. For example
X=xcoord; PRENVELOPE=!e( PRT(xcoord; 5) );\
GRENVELOPE=!e( GRT(ntry; 5) ); NTRIES=ntry]
GRENVELOPE = !e( GRT(NTRIES; 60) )
The random samples can be saved, in variates, using the
NUMBERS parameter. You can use the
SEED parameter to supply a seed to use in the
CALCULATE directive for the sequences of the random numbers used to generate the random values from g(x) and U[0,1]. The default,
SEED=0, continues an existing sequence of random numbers, if any of the random-number functions has already been used in the current Genstat run. If, however, this is the first time that these functions have been used, Genstat picks a random seed.
PLOT option controls the graphs produced by
GREJECTIONSAMPLE, with settings:
||to plot the density f(x) and the envelope M×g(x), and|
||to plot a histogram of the selected sample from f(x).|
By default these are both plotted.
For further details of the method see e.g. Carlin & Louis (2000) or Robert & Casella (2004).
Carlin, B.P. & T.A. Louis (2000). Bayes and Empirical Bayes methods for Data Analysis. Chapman & Hall, London.
Robert, C.R. & Casella, G. (2004). Monte Carlo Statistical Methods, Second Edition. Springer, New York.
Commands for: Calculations and manipulation.
CAPTION 'GREJECTIONSAMPLE examples'; STYLE=meta " 100 samples from the density p(x) = 0.5 for -1 <= x <= 1 " GREJECTIONSAMPLE [NVALUES=100; X=x; XLOWER=-1; XUPPER=1; MULTIPLIER=2.5;\ PRDENSITY=!e(0.5*(ABS(x).LE.1))] s1; SEED=171246 " 100 samples from the density p(x) = 1/3 for -1 <= x <= 0 = 2/3 for 0 < x <= 1 " GREJECTIONSAMPLE [NVALUES=100; X=x; XLOWER=-1; XUPPER=1; MULTIPLIER=3;\ PRDENS=!e((x.LE.0.AND.x.GE.-1)/3+(x.GT.0.AND.x.LE.1)*2/3)] s2 " 100 samples from the density EXP(-1*ABS(x))*ABS(SIN(x)) " GREJECTIONSAMPLE [NVALUES=100; X=x; PRDENSITY=!e(EXP(-1*ABS(x))*ABS(SIN(x)));\ PRENVELOPE=!e(PRT(x; 1.5)); GRENVELOPE=!e(GRT(200; 1.5));\ MULTIPLIER=3] s3 LPHISTOGRAM s1,s2,s3