Professor of Stochastic Modelling
	School of Mathematics & Statistics
	Newcastle University

Rejection sampling

• An example of rejection sampling

• Generating beta random quantities

Exercises

1.

Write a vectorised R function for beta variate generation.

2.

Write a wrapper so that it generates approximately the number you ask for.

3. (for postgrads)

See if you can modify the wrapper so that it always gives exactly the number you ask for.

4. (for keen undergrads and all postgrads)

For the Normal model with unknown mean and variance, use rejection sampling to simulate from the posterior based on candidates from the prior (Bayes theorem via the rejection method). Note that you will need to know the maximum likelihood estimators for the normal unknown mean and variance model in order to do this!

5. (for postgrads)

See if you can write some C code for Bayes theorem via the rejection method.

Project 1

Use the information in your lecture notes on the envelope method to generate t random variates with given degrees of freedom. Verify that the scaled Cauchy density with density proportional to 1/(4+x^2) forms an envelope for the t density if both are scaled to have a maximum of 1 at zero. You should investigate this with some overlay plots in Maple. The following commands carry out the verification for 2 degrees of freedom. You should also check the general case analytically.

with(stats);
tpdf:=(x,v)->statevalf[pdf,studentst[v]](x);
scpdf:=x->1/(4+x^2);
plot({scpdf(x)/scpdf(0),tpdf(x,2)/tpdf(0,2)},x=-6..6);

The t density is proportional to (1+(x^2)/v)^(-(v+1)/2), where v is the degrees of freedom. Include overlay plots demonstrating that the envelope works for a range of degrees of freedom. You may wish to demonstrate theoretically that the envelope is valid. Assuming that the envelope is valid, describe how the rejection algorithm is implemented. Implement the algorithm in R (for general v), and compare a histogram of roughly 20000 simulated values for v=2 with the actual density - the R commands lines and dt may be useful here. Include your actual R functions in your report. What is your empirical acceptance probability for the v=2 case? Also work out the theoretical acceptance probability for this scheme (you may wish to use Maple to help with the integrals), and make sure your empirical value is very close to this.

Your report should be coherent - not just a list of answers. It should be 5-10 sides of A4, including all plots, tables and code, and produced in Word or LaTeX. The submission deadline is the end of week 11, but bear in mind that you will get another project in two weeks time.

		darren.wilkinson@ncl.ac.uk
		http://www.staff.ncl.ac.uk/d.j.wilkinson/