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