代写lectures留学生、代做R编程语言、代写histogram、R程序设计调试

日期：2019-10-12 11:40

Lab 1: Simulation:

Method of inversion and rejection sampling

Hand-in date: Wednesday, 9th of October at 12pm.

Aim of this practical

This practical is focused on simulating from exponential and gamma distributions using the

techniques we’ve learned in lectures, namely the inversion method and rejection sampling.

Remember that the basic ’building block’ for each of these methods is the ability to sample

from the uniform distribution, U(0, 1). In R the function runif(n) generates n pseudo-random

numbers from U(0, 1). Note, that simulating from many distributions is already built into R, but

for the purposes of this lab assignment, let’s assume that this isn’t the case!

Exponential and gamma distributions

The density of the exponential distribution is

fλ(x) = λ exp(−λx),

for x > 0 and λ > 0. The CDF is written as

Fλ(x) = 1 − exp(−λx).

The gamma density can be written as

f(α,λ)(x) = 1

for x > 0 and α, λ > 0. The exponential distribution is actually a special case of the gamma

distribution with α = 1. Moreover, for X1, . . . , Xk idd ∼ Ex(λ), and

Y = X1 + · · · + Xk

Y ∼ Ga(k, λ).

Method of inversion

Recall from lectures, that if U ∼ U(0, 1), then

− log(1 − U)/λ ∼ Ex(λ).

We can therefore generate an exponential random deviate as

1 lambda <− 1 # D ef i n e t h e r a t e

2 u <− r u n i f ( 1 ) # G e n e r a t e a U( 0 , 1 ) number

3 x <− −log(1−u )/lambda # T ra n sfo rm to an e x p o n e n t i a l

4 x # Output to s c r e e n

✝ ✆

We can extend this code to sample more than one realisation:

✞ ☎

1 n <− 20 # D ef i n e t h e sam ple s i z e

2 lambda <− 1 # D ef i n e t h e r a t e

3 u <− r u n i f ( n ) # V e c to r of n U( 0 , 1 ) numbe rs

4 x <− −log(1−u )/lambda # T ra n sfo rm to v e c t o r of e x p s

5 x # Output to s c r e e n

✝ ✆

Note, that it is also possible to generate a single exponential random variate and use a loop to

repeat this each time. This is very inefficient, however. The code above makes use of vectors –

R is much more efficient if everything is vectorised.

We can also wrap our code inside a function as follows:

✞ ☎

1 i n v . exp <− func tion ( n , lambda ) { #D ef i n e t h e f u n c t i o n i n v . exp

2 u <− r u n i f ( n ) # G e n e r a t e a U( 0 , 1 ) number

3 x <− −log(1−u )/lambda # T ra n sfo rm to an exp

4 x # R e t u r n x

5 }

6 i n v . exp ( n=20, lambda=1)

✝ ✆

We can now, informally at least, check that our code samples correctly, by comparing a histogram

of sampled values to a plot of the density:

✞ ☎

1 x <− i n v . exp ( n=1000, lambda =0.5) # G e n e r a t e a l a r g e sam ple

2 h i s t ( x , f r e q=FALSE ) # P l o t a h i s t o g r am

3 t <− 0:150 /10 # C r e a t e a s e q u e n c e

0 , 0 . 1 , . . . , 1 5

4 l i n e s ( t , dexp ( t , r a t e =0.5) , lwd=2) # Add t h e d e n i s t y to t h e

h i s t o g r am

✝ ✆

We could also use a Kolmogorov-Smirnov test to test the null hypothesis that our sample is from

an Ex(0.5) distribution:

✞ ☎

1 k s . t e s t ( x , pexp , r a t e =0.5) # K−S t e s t

✝ ✆

As before, we have that Y ∼ Ga(k, λ) with k ∈ R

Y = X1 + · · · + Xk

2

where Xi

iid ∼ Ex(λ). We can use this fact to simulate from a Ga(k, λ) distributions with

integer-valued k:

✞ ☎

1 k <− 4

2 lambda <− 1

3 x <− i n v . exp ( n=k , lambda=lambda ) # Draw an exp number

4 y <− sum( x )

5 y

✝ ✆

We can also draw n random Ga(k, λ) random numbers by first creating a matrix with n rows and

k columns, consisting of the nk random exponential numbers, and then add up each row:

✞ ☎

1 n <− 20 # d e f i n e n

2 k <− 4 # d e f i n e k

3 lambda <− 1 # d e f i n e lambda

4 x <− matrix ( i n v . exp ( n=n∗k , lambda=lambda ) , ncol=k )

5 y <− apply ( x , 1 , sum) # sum up eac h row

6 y # o u t p u t y

✝ ✆

We can also wrap this code inside a function

✞ ☎

1 i n v .gamma. i n t <− func ti on ( n , k , lambda ) {

2 x <− matrix ( i n v . exp ( n=n∗k , lambda=lambda ) , ncol=k )

3 apply ( x , 1 , sum)

4 }

5 y <− i n v .gamma. i n t ( 2 0 , 4 , 1 ) # t e s t o u t t h e f u n c t i o n

✝ ✆

3

Exercises

1. A factory contains 8 machines. Any one of the machines can breakdown independently of

the other machines following an exponential distribution with a mean value of 10 days.

(a) What is the distribution of the time until the first machines breakdowns? [10]

(b) Implement an algorithm to simulate from this distribution, the time until the first

machine breaks down, given a random variate u from U[0, 1]. [10]

(c) Estimate the expected time until the first machine breakdowns using Monte Carlo

sampling. What is the analytic (true) value of this expected value? [10]

(d) Estimate, using Monte Carlo sampling, the probability that the time until the first

machine breakdowns is at least 5 days. [10]

2. Consider the following: If X ∼ Po(λ) and Xi ∼ Ex(λ), i ∈ N, then

P(X = 0) = P(X1 > 1)

and

P(X = k) = P(X1 + · · · + Xk ≤ 1 < X1 + · · · + Xk+1), for k = 1, 2, . . . ,

(a) Based on this identity, implement an algorithm to simulate from a Po(λ) distribution

using simulations from an Ex(λ) distribution. [15]

(b) Generate a sample of size 100 drawn from this algorithm for each of the following

values of λ = 2, 5, 10. Examine estimated means and variances in each case and

comment. [5]

(c) What are limitations of this algorithm? [5]

3. Consider the standard normal density f(x) = (2π)

−1/2

exp(−x

2/2) and a standard doubleexponential

density, h(x) = 1

2

exp(−|x|), for x ∈ R.

(a) Describe how to perform rejection sampling algorithm to sample X ∼ f using draws

Y ∼ h, ie, where h(x) is used as an envelope function. Describe in particular how

this algorithm may be implemented to minimise the probability of rejection. [15]

(b) Is it possible to use rejection sampling to simulate Y ∼ h using f(x) as an envelope

function? [5]

(c) Implement in R the algorithm developed in (a). [15]

4