GR5206-Homework 6 Solved

Goals: Simulating probability distributions using the Inverse Transform Method and AcceptReject Method. Built in R distribution functions. Estimate the mathematical constant π using Monte Carlo techniques. More practice writing functions and plotting.

First set your seed as 0, i.e., set.seed(0) Part 1: Inverse Transform Method

Consider the Cauchy random variable X with probability density function

 .

1.    Let U be a uniform random variable over [0,1]. Find a transformation of U that allows you to simulate X from U. Note: I will solve this in class.

2.    Write a R function called cauchy.sim that generates n simulated Cauchy random variables. The function should have the single input n and should use the inversetransformation from Part 1. Test your function using 10 draws.

3.    Using your function cauchy.sim, simulate 1000 random draws from a Cauchy distribution. Store the 1000 draws in the vector cauchy.draws. Construct a histogram of the simulated Cauchy random variable with fX(x) overlaid on the graph. Note: when plotting the density curve over the histogram, include the argument prob = T. Part 2: Reject-Accept Method

Problem 2
Let random variable X denote the temperature at which a certain chemical reaction takes place. Suppose that X has probability density function

(1)                                              

                                                                                    0     otherwise

Perform the following tasks:

4.    Write a function f that takes as input a vector x and returns a vector of f(x) values.

Plot the function between −3 and 3. Make sure your plot is labeled appropriately.

5.    Determine the maximum of f(x) and find an envelope function e(x) by using a uniform density for g(x). Write a function e which takes as input a vector x and returns a vector of e(x) values.

6.    Using the Accept-Reject Algorithm, write a program that simulates 10,000 draws from the probability density function f(x) from Equation 1. Store your draws in the vector f.draws.

7.    Plot a histogram of your simulated data with the density function f overlaid in the graph. Label your plot appropriately.

Problem 3: Reject-Accept Method Continued
Consider the standard normal distribution X with probability density function

 .

In this exercise, we will write a function named normal.sim that simulates a standard normal random variable using the Accept-Reject Algorithm.

Perform the following tasks:

8. Write a function f that takes as input a vector x and returns a vector of f(x) values.

Plot the function between −5 and 5. Make sure your plot is labeled appropriately. 9. Let the known density g be the Cauchy density defined by pdf

 .

Write a function e that takes as input a vector x and constant alpha (0 < α < 1) and returns a vector of e(x) values. The envelope function should be defined as e(x) = g(x)/α.

10.    Determine a “good” value of α. To show your solution, plot both f(x) and e(x) on the interval [−10,10].

11.    Write a function named normal.sim that simulates n standard normal random variables using the Accept-Reject Algorithm. The function should also use the InverseTransformation from Part 1. Test your function using n=10 draws.

12.    Using your function normal.sim, simulate 10,000 random draws from a standard normal distribution. Store the 10,000 draws in the vector normal.draws. Construct

a histogram of the simulated standard normal random variable with f(x) overlaid on the graph. Note: when plotting the density curve over the histogram, include the argument prob = T.

Part 3: Simulation with Built-in R Functions

Consider the following “random walk” procedure:

•     Start with x = 5

•     Draw a random number r uniformly between −2 and 1.

•     Replace x with x + r

•     Stop if x ≤ 0

•     Else repeat

Perform the following tasks:

13.    Write a while() loop to implement this procedure. Importantly, save all the positive values of x that were visited in this procedure in a vector called x.vals, and display its entries.

14.    Produce a plot of the random walk values x.vals from above versus the iteration number. Make sure the plot has an appropriately labeled x-axis and and y-axis. Also use type="o" so that we can see both points and lines.

15.    Write a function random.walk() to perform the random walk procedure that you implemented in question (9). Its inputs should be: x.start, a numeric value at which we will start the random walk, which takes a default value of 5; and plot.walk, a boolean value, indicating whether or not we want to produce a plot of the random walk values x.vals versus the iteration number as a side effect, which takes a default value of TRUE. The output of your function should be a list with elements: x.vals, a vector of the random walk values as computed above; and num.steps, the number of steps taken by the random walk before terminating. Run your function twice with the default inputs, and then twice times with x.start equal to 10 and plot.walk = FALSE.

16.    We’d like to answer the following question using simulation: if we start our random walk process, as defined above, at x = 5, what is the expected number of iterations we need until it terminates? To estimate the solution produce 10,000 such random walks and calculate the average number of iterations in the 10,000 random walks you produce. You’ll want to turn the plot off here.

17.    Modify your function random.walk() defined previously so that it takes an additional argument seed: this is an integer that should be used to set the seed of the random number generator, before the random walk begins, with set.seed(). But, if seed is NULL, the default, then no seed should be set. Run your modified function random.walk() function twice with the default inputs, then run it twice with the input seed equal to (say) 33 and plot.walk = FALSE. Part 4: Monte Carlo Integration

The goal of this exercise is to estimate the mathematical constant π using Monte Carlo Integration. Consider the function

                                                                              (√1  − x2      0 ≤ x ≤ 1

(2)                                             g(x) =

                                                                                       0     otherwise

The above function traces out a quartile circle over the interval [0,1].

Perform the following tasks:

18.    Run the following code:

g <- function(x) {

return(sqrt(1-x^2))

}

plot(seq(0,1,.01),g(seq(0,1,.01)),type="l",col="purple")

The above code should produce the plot of a quarter circle.

19.    Identify the true area under the curve g(x) by using simple geometric formulas.

20.    Using Monte Carlo Integration, approximate the mathematical constant π within a 1/1000 of the true value. When performing this simulation, make sure to choose a probability density function that has support over the unit interval, i.e. uniform(0,1) or beta(α,β).