STAT - 611 ASSIGNMENT 4
Due Tuesday Nov 06, 2018 in class
Remember (from syllabus): Solutions must be neatly word-processed (using Latex or a word processor)
and stapled. Please annotate your work with brief, clear sentences explaining your approach and interpreting
your results (you are not expected to write full-blown data analysis reports however). For assignments
involving mathematical manipulations, students can write answers by hand provided penmanship is neat;
illegible answers will be marked as incorrect.
*All these questions require you write computer programs. Please write a report with your code, and an
explanation of what you’ve done. Additionally upload the les containing your computer code to Canvas.
Include instructions detailing how to use them in your report.
1. The double exponential function has density
fL(x) = 2 exp( jxj)
for > 0.
(a) Derive an algorithm to generate variates from a double-exponential distribution using the inversion
method.
(b) Derive an accept-reject algorithm for obtaining samples from the N(0;1) distribution using the
double exponential with = 1 as a proposal. Make sure that you use the optimal enveloping
constant.
(c) Implement a function in R that generates variates from the N(0;1) distribution using your accept-
reject algorithm. Use it for generating 1000 variates. Plot some interesting plots (histogram,
density, etc.). Test your samples for normality using a Kolmogorov-Smirnov test ks.test (see R
help). Comment your results.
2. Implement your accept-reject standard normal sampler in C++. You will have to implement a function
rnormal with prototype,
double rnormal();
that generates a single draw from a standard normal distribution using your rejection algorithm.
For this you will need to generate uniform. variates. Don’t use the default prng from the C library|
rand(); it is a low-quality algorithm, unsuitable for statistical applications. Instead you will use the
Mersenne twister algorithm. You can found a C++ implementation in the "MersenneTwister.h" le,
uploaded to the " les/code" section in Canvas. To use it, download the le and save it in the same
directory of your code. Include the following lines after your "include" statements in your source code
le.
#include"MersenneTwister.h"
#include
MTRand rng;
These lines declare and initialize a global object called rng, of class MTRand, that can be used to
generate pseudo-random numbers. To get a pseudo-random number on the interval (0;1) you can use
the class method "rand()" within any of your functions. For example:
double u = rng.rand();
1
If your object rng is global, a call to rng.rand() will generate a di erent pseudo-random number each
time.
Test your code with the following main function (you need to include the header le stdio.h for this
to work)
int main(){
FILE* f = fopen("numbers.dat", "w");
for (int i = 0; i x <- scan("numbers.dat")
Again generate some plots and perform. a test of normality. Comment your results.
3. Now you will write an R interface to your C++ code. Write a function my norm, callable from R
using the .C mechanism (return type void, and pointer arguments). This function should generate an
arbitrary number (speci ed by the user) of random variates using your rnormal function, and return
them to R. Don’t forget to include the \extern"C"fg" declaration when compiling with g++.
Write and R wrapper function to call your C function from R. Use it to generate 1,000,000 samples
and compare speed with your R implementation from the previous question using the system.time
function (read the corresponding help le to learn about this function). Comment your results.