ST3247程序语言代写、代写R课程编程

日期：2022-03-31 10:05

Simulation, Semester II 2021-2022

ST3247: Assignment 21

Due date: April 10th (Sunday of Week 12)

There are THREE questions in this question sheet.

Your answer sheet should contain answers to ALL questions. Do not leave answers in

your code file. Marking mainly depends on your answer sheet.

Please submit an answer sheet with an R file in the submission folder separately, and name

them with the same name. Please do NOT use compressed files.

The R file should contain your coding work as a supplementary file. Please make it ready

for reproductions.

For all the questions, you can only use standard uniform random variables unless

there are particular statements in the problem.

1. Suppose X1, X2, · · · , Xn ～ N(0, θ) where θ ～ Gamma(3, 0.5), that

p(θ) =

θ2e?θ/2

23Γ(3)

, θ > 0.

Note. For this question, you can use any random variable generation function in R, such

as rnorm(), rgamma() and so on. You can also use gamma() to calculate Γ(x) for any x.

(a) Find the posterior distribution that p(θ|x1, x2, · · · , xn), subject to a constant cp.

(b) We are interested in the mean of the posterior distribution, which is

E =

∫ ∞

0

θ ? p(θ|x1, x2, · · · , xn)dθ.

The posterior distribution is hard to sample even with the rejection algorithm. So

we want to use the importance sampling algorithm.

Suppose we have the random variables Y1, Y2, · · · , Ym ～ q, where q is a distribution.

Please write out the estimation in terms of m, p(θ) = p(θ|x1, x2, · · · , xn), q, and Yi’s.

1All rights reserved by NUS. Reproduction or distribution of lecture notes/tutorials/quiz/exam without the

written permission of the sponsor is prohibited.

1

(c) To get rid of the constants, let p = cpp0 and q = cqq0, where cp and cq are constants

and p0, q0 are kernel functions. Prove that∫

f(θ)p(θ)

q(θ)

q(θ)dθ =

∫ f(θ)p(θ)

q(θ) q(θ)dθ∫ p(θ)

q(θ)q(θ)dθ

=

∫ f(θ)p0(θ)

q0(θ)

q(θ)dθ∫ p0(θ)

q0(θ)

q(θ)dθ

.

(d) According to part (b) and (c), find the estimation in terms of m, p0, q0, and Yi’s.

Further, express them in terms of log p0, log q0 and Yi’s.

(e) We set the function q to be Gamma(cs2, c), i.e.

g(θ) =

ccs

2

θcs

2?1e?cθ

Γ(cs2)

, θ > 0.

Here, s2 is the sample variance and c > 0 is a parameter to decide. We use c to make

the function flat.

With the density functions of p and q, write R functions to calculate log p0 and log q0.

(f) With all the previous definitions of p, q and the discussion, write a pseudo-code to

estimate E.

(g) Generate 100 data points with distribution N(0, 5), denoted by x1, x2, · · · , x100. It

means the ground truth is θ = 5.

With your generated data, run your code to estimate E with c ranges from 0.04 to

20 with step size 0.04, m = 1000. Draw one plot of E? versus c and comment.

2

2. Messages arrive at a communications facility in accordance with a Poisson process having

a rate of λ = 2 per hour. The facility consists of three channels, and an arriving message

will either go to a free channel if any of them are free or else will be lost if all channels

are busy. The amount of time that a message ties up a channel is a random variable

that depends on the weather condition at the time the message arrives. Specifically, if the

message arrives when the weather is “good,” then its processing time is a random variable

having distribution function

F (x) = x, 0 < x < 1

whereas if the weather is “bad” when a message arrives, then its processing time has

distribution function

F (x) = x3, 0 < x < 1.

Initially, the weather is good, and it alternates between good and bad periods— with the

good periods having fixed lengths of 2 hours and the bad periods having fixed lengths of

1 hour. (Thus, for example, at time 5 the weather changes from good to bad.) Suppose

we are interested in the distribution of the number of lost messages by time T = 100.

(a) Set up a vector with length T , that decide the weather at time point 1 ≤ t ≤ T .

(b) Define the events and variables that enable us to use the DES approach.

(c) Write the pseudo code of how to update the variables for each iteration

(d) Implement the R code. Run your program to give an interval estimate with 95%

confidence of the expected number of lost messages in the first 100 hours of operation.

(No need to include code here.)

(e) Now I want to check the code under a simplified case. For the following, assume that

we are always in GOOD weather and there is only ONE channel.

i. What is the inter-arrival time between two messages?

ii. The status of this channel will be Busy, Idle,Busy, Idle, · · · . Busy means it is

processing a message and Idle means there is no message under processing now.

Let tB denote the length of one busy period. What is the distribution of tB?

What is E[tB]?

iii. Let tI be the length of one idle period. What is the distribution of tI?

Hint. tI is the time until next arriving message.

iv. The messages arrived during the busy period are lost. Let N denote the number

of lost messages. Let TB be the total busy time of the channel until T = 100.

3

Prove that E[N ] = λE[TB].

v. According to the property of stochastic processes, in the long run, the proportion

of all busy periods is approximately E[tB]/(E[tB] + E[tI ]). Let TB be the total

time that the channel is busy until T = 100. What is E[TB]?

vi. Find E[N ] according to parts (iv) and (v).

vii. Adjust your code for this case. Run 1000 simulations for T = 100. Compare

your results with the results in part (vi).

4

3. Let Q be a symmetric matrix, that is, qij = qji for all i, j. It satisfies the requirements for

a transition probability matrix, that every entry is non-negative and the row sum is 1.

Consider a Markov chain based onQ, where the transitions due to the following procedure,

Step 1. Generate Y ～ qi, where i is the current state. qi is the distribution defined by

i-th row of Q.

Step 2. Generate U ～ Unif(0, 1).

Step 3. If U ≤ bYbi+bY , set next step to be Y

Step 4. Otherwise, set next step to be the current state i.

Here, bj ’s are specified positive numbers.

(a) Say X0 = i, what is the distribution of X1?

(b) Find the transition probability matrix of the current Markov Chain.

(c) Actually this MC is time reversible. Please set up the local balanced equations and

find the corresponding limiting distribution pi.

(d) Consider the distribution on the permutations of (1, 2, 3, · · · , 10). The probability of

a permutation (x1, x2, · · · , x10) is that

P ((x1, x2, · · · , x10)) = C ?

10∑

i=1

xi/i,

where C is the normalising constant.

Please identify a procedure according to the results of this problem. Identify Q and

bi’s, and then write down the pseudo code.

(e) With the pseudo code in part (d), please run simulations to find the expectation of∑10

i=1 i ? xi. Write down the steps you get your samples (take that in part (d) as a

single function) and the result.