R Distribution Sampling and Visualization

Tags: R, Data Sampling, Visualization, Histogram, Monte Carlo, ggplot, rnorm, rejection sampler

Goal: Sampling from different distributions via the inverse-cdf method and the rejection sampler.

Fsor this homework, submit your clearly commented R code for this assignment electronically via Canvas. Note the deadline on Canvas. In addition to the R code, you must also turn in your typed-in answers. Please submit your R code and the typed up answers as separate documents. Points will be deducted if R code is very poorly formatted.

  1. Draw 10,000 independent samples from a Pareto( = 3; = 5) density by using the inverse cdf method. The Pareto pdf is

​ f(x) =x+1 ; a < x < 1; > 0; > 0 (1)

​ (a) Clearly describe your algorithm in "pseudo-code", i.e., write out your algorithm for this problem briey, systematically, in words, filling in mathematical details where necessary. Pseudo-code example (Note, this is not pseudo-code for this problem):

i. Initialize the sequence by setting X0 = 1:

ii. For i = 1 to N (where N = 1000 here), repeat: ...

iii. Set T =PNi=1 Xi

iv. Return ^ = T=N.

​ (b) Plot a histogram of your samples, along with a smoothed density plot. You will need commands: plot, lines, and density. Make sure you label the plot well so it is easy to read.

​ (c) Suppose X Pareto( = 3; = 5). Estimate E(X) and E(X2) via Monte Carlo. Report your estimate along with the Monte Carlo standard error of your estimate.

  1. Draw 10,000 independent samples from a Beta( = 2; = 8) density by using a rejection sampler. The Beta pdf is

    ​ f(x) =􀀀( + )􀀀()􀀀() x􀀀1(1 􀀀 x)􀀀1; x 2 (0; 1); > 0; > 0 (2)

    You are only allowed to use uniform random variables (runif). (There is a function rbeta in R that would solve this problem for you, but defeat the purpose of the assignment which is to learn how to construct a rejection sampler...)

    (a) Before you write your algorithm, make sure you state the dierent conditions that you have satised in order for your algorithm to work. For example, the proposal pdf q needs to satisfy certain conditions. Make sure you show that the conditions are satised.

    (b) Clearly describe your algorithm in pseudo-code.

    (c) Plot a histogram of your samples, along with a smoothed density plot. Make sure you label the plot well so it is easy to read.

    (d) Suppose Y Beta( = 2; = 8). Use Monte Carlo to estimate E(Y 2) and E(Y 3). Report your estimates along with associated Monte Carlo standard errors

  2. Compare dierent ways to simulate a normal random variable. You will simulate a standard normal random variable via transformations, rejection sampling, and standard R function. You can time events by using this following commands:

time.start = Sys.time()

... CODE HERE ...

time.end = Sys.time()

time.end - time.start

For n = 100, n = 1000, n = 10; 000, and n = 100; 000, record the time it takes to simulate the normal random variables using the following 3 methods. Plot the times for the three methods all on the same graph. Each method should be a dierent color and connected by a dierent line (three separate lines, each a new color). Plotting can be down using the plot, points and lines commands, or via ggplot in the ggplot2 package, which is sometimes dicult to learn, but if you wish, a help page is found here: https://ggplot2.tidyverse.org/.

(a) Simulate n standard normal random variables using rejection sampling. Use a cauchy random variable (rcauchy in R) as your proposal density, since the cauchy random variable has heavier tails than a normal random variable.

(b) Simulate n standard normal random variables using the fact that if U1 U(0; p 1) and U2 U(0; 1) are independent, then Z1 = 􀀀2 log(U1) cos(2U2) and Z2 =p􀀀2 log(U1) sin(2U2) are standard normal random variables and are independent of each other.

(c) Simulate n standard normal random variables using rnorm in R.