STAT - S611代做、C/C++编程调试、data留学生代写、c++语言代写

日期：2020-04-22 10:41

STAT - S611 ASSIGNMENT 6

Due Sunday April 26. PLEASE TURN YOUR WORK ELECTRONICALLY VIA CANVAS

• 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 files containing your computer code to

Canvas. Include instructions detailing how to use them in your report.

1. The double exponential function has density.

(a) Derive an algorithm to generate variates from a double-exponential distribution using the inverse

CDF 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 acceptreject

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” file,

uploaded to the ”files/code” section in Canvas. To use it, download the file and save it in the same

directory of your code. Include the following lines after your ”include” statements in your source code

file.

#include"MersenneTwister.h"

#include<math.h>

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();

If your object rng is global, a call to rng.rand() will generate a different pseudo-random number each

time.

Test your code with the following main function (you need to include the header file stdio.h for this

to work)

1

int main(){

FILE* f = fopen("numbers.dat", "w");

for (int i = 0; i < 1000; i++){

fprintf(f, "%1.15e\n", rnormal());

}

fclose(f);

}

This program will generate 1000 normal variates using your rnorm function (assuming your function

is correct!) and will write them into a text file called ”numbers.dat”.

Load these numbers into an object in R by using the scan function (read about this function in the R

help files).

>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 (specified by the user) of random variates using your rnormal function, and return

them to R. Don’t forget to include the “extern"C"{}” 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 file to learn about this function). Comment your results.

2