Generates random samples using rejection sampling (W. van den Berg).
Options
PLOT = string tokens |
What to plot (density, sample); default dens, samp |
|---|---|
NVALUES = scalar |
Size of each random sample; no default, must be set |
PRDENSITY = expression structure |
Calculation defining the probability density function f(x) to sample; no default, must be set |
X = identifier |
Data structure used inside PRDENSITY for the x-coefficient of the density function f(x) no default, must be set |
XLOWER = scalar |
Lower bound of the region in which f(x) is non-negligible; default -10 |
XUPPER = scalar |
Upper bound of the region in which f(x) is non-negligible; default 10 |
PRENVELOPE = expression structure |
Calculation defining the probability density function g(x) used to generate the sample; default !e(PRT(X; 60)) |
GRENVELOPE = expression structure |
Calculation to sample from the probability density g(x) used to generate the sample (note, PRENVELOPE and GRENVELOPE must either be both set, or both unset); default !e(GRT(NTRIES; 60)) |
MULTIPLIER = scalar |
Multiplier M used in the definition of the envelope M × g(x) that must always be greater than f(x); default 10 |
NTRIES = scalar |
Number of random samples to take in each sampling step; default * i.e. determined automatically |
Parameters
NUMBERS = variates |
Saves each random sample |
|---|---|
SEED = scalars |
Seed to use for the random numbers used to generate each random sample; default 0 |
Description
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
SCALAR ntry
GREJECTIONSAMPLE [...\
X=xcoord; PRENVELOPE=!e( PRT(xcoord; 5) );\
GRENVELOPE=!e( GRT(ntry; 5) ); NTRIES=ntry]
By default,
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.
The PLOT option controls the graphs produced by GREJECTIONSAMPLE, with settings:
density |
to plot the density f(x) and the envelope M×g(x), and |
|---|---|
sample |
to plot a histogram of the selected sample from f(x). |
By default these are both plotted.
Options: PLOT, NVALUES, PRDENSITY, X, XLOWER, XUPPER, PRENVELOPE, GRENVELOPE, MULTIPLIER, NTRIES.
Parameters: NUMBERS, SEED.
Method
For further details of the method see e.g. Carlin & Louis (2000) or Robert & Casella (2004).
References
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.
See also
Directive: CALCULATE.
Procedures: GRCSR, GRLABEL, GRMULTINORMAL, GRTHIN, GRTORSHIFT, SAMPLE, SVSAMPLE.
Functions: GRBETA, GRBINOMIAL, GRCHISQUARE, GRF, GRGAMMA, GRHYPERGEOMETRIC, GRLOGNORMAL, GRNORMAL, GRPOISSON, GRSAMPLE, GRSELECT, GRT, GRUNIFORM.
Commands for: Calculations and manipulation.
Example
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