1. Home
2. GREJECTIONSAMPLE procedure

# GREJECTIONSAMPLE procedure

Generates random samples using rejection sampling (W. van den Berg).

### Options

`PLOT` = string tokens What to plot (`density`, `sample`); default `dens`, `samp` 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 `PRDENSITY` for the x-coefficient of the density function f(x) no default, must be set 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 `!e(PRT(X; 60))` 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 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 `*` i.e. determined automatically

### Parameters

`NUMBERS` = variates Saves each random sample 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 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.

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
```
Updated on March 7, 2019