**Unformatted text preview: **Copyright c Peter Glynn All rights reserved. MS&E 221 Spring 2020 Problem Set 3 MS&E 221
Due: Friday, May 15 11:59PM (Perishable Product)
Consider a grocery store selling a perishable product. We denote the number of items of the
product at 8 am of the n-th day as X . The store's manager starts working at 8am, and she will
pick out and discard any spoiled product before the store opens at 9am. The number of spoiled
products follows a Binomial(X , p) distribution while the daily demand of the product follows a
P oisson(λ) distribution. If the product becomes out-of-stock, the store incurs lost sales. The store
follows an (s, S) policy and checks the inventory Y after the store closes at 9pm. Specically, if the
manager nds the inventory Y < s, a replenishment order will be placed immediately for S − Y
items, where s = 4 and S = 20. The order arrives before 8am the next day. Assume X = 19,
p = 0.2, and λ = 4.
(a) Compute the transition probabilities of X . Report P (X = 10|X = 15), P (X =
15|X = 10), and P (X
= 20|X = 6).
(b) To generate a better contract with the product supplier, we want to know the expected time
until the arrival of the rst replenishment order. Compute this quantity using FTA.
(c) To assess the product's protability, we are interested in the ratio between the total discarded
product D and the total sold product hS over
a selling season of 30 days. Simulate X for
i
100 trials and report your estimate of E
with a 90% condence interval.
Question 3.1 : n n n n n 0 n n n+1 n+1 n n+1 n 30 30
D30
S30 n Answer: (a) Let Z , Z , Z be the inventory level at 8am, 9am, 9pm on the n-th day respectively.
By denition, we know X = Z . We can compute the transition probabilities for Z ,
(
p (1 − p) if j ≤ i
P , P (Z
= j|Z = i) =
0
if j > i e
if 0 < j ≤ i P
P , P (Z
= j|Z
= i) = 1 −
e
if j > i 0
if j = 0
3n 3n+1 3n+2 n 1 2 3n n 3n+1 3n+2 3n 3n+1 i
j i−j j λi−j −λ
(i−j)!
i
λi−j −λ
j=1 (i−j)! 1 Copyright c Peter Glynn All rights reserved.
P3 , P (Z3n+3 = j|Z3n MS&E 221 Spring 2020 1
= i) = 1 0 The transition probability matrix for X is P P P .
n if j = i ≥ s
if j = S, i < s
otherwise 1 2 3 P (Xn+1 = 10|Xn = 15) = 0.1287, P (Xn+1 = 15|Xn = 10) = 0,
0.8971. (b) Dene T = min{n ≥ 0 : X
For x ∈ C = {s, . . . , S − 1}, n = S}, u(x) = Ex T, C c = {S}.
u(x) = 1 + X and P (X n+1 = 20|Xn = 6) = We know u(S) = 0 by denition. P (x, y)u(y) y∈C where P is the transition probability matrix for X dened in the rst part. Here the reason
that the set C starts at 4 is X will never take a value smaller than s. We solve it numerically
and get u(19) = 2.837.
(c) The estimate is 0.6891 and the condence interval is [0.67228, 0.70585].
n n (Cayley Tree) The Cayley tree of degree k and depth d is a tree with undirected
edges, on 1 + k + k(k − 1) + · · · + k(k − 1) nodes, where every non-leaf node has k adjacent edges.
(Assume k > 2 and d > 0.)
Assume that a particle does a random walk on a Cayley tree of degree k and depth d: from any
node, the particle moves to one of the neighbors of that node chosen uniformly at random. The
following notation may be helpful in this question. For a node i, let N (i) denote the neighbors of i;
note that any non-leaf node has exactly k neighbors. For a leaf node v, let p(v) denote the parent
of v in the graph; i.e., p(v) is the unique node adjacent to v.
Question 3.2 : d−1 2 Copyright c Peter Glynn All rights reserved. MS&E 221 Spring 2020 Figure 1: A Cayley tree with k = 4 and d = 3; note that it has 1 + k + k(k − 1) + k(k − 1)
nodes, and that each non-leaf node has k = 4 adjacent edges.
(a) Describe the movement of the particle as a Markov chain.
(b) Is this Markov chain irreducible?
(c) Is this Markov chain aperiodic?
(d) Compute the stationary distribution of the random walk. 2 = 53 Answer: (a) For any non-leaf node i, P (i, j) = for j ∈ N (i), and 0 otherwise. For any leaf node i,
P (i, j) = 1 for j = p(i), and 0 otherwise.
(b) Yes. There is a path with positive probability from any node to any other node.
(c) No. Any path from a node i back to itself must take an even number of steps.
(d) It suces to check the detailed balance conditions. By symmetry, all nodes at depth i have at
same equilibrium probability, which we denote by π . The detailed balance is equivalent to
1
k i 1
1
= π1 ·
k
k
1
1
π1 · = π2 ·
k
k
...
1
πd−1 · = πd · 1
k
π0 · 3 Copyright c Peter Glynn All rights reserved. MS&E 221PSpring 2020 From these we immediately get π = π = · · · = π = kπ . Let c = π , then by
1, where N is the number of nodes at depth i, we have:
0 1 d−1 0 d d
i=0 πi Ni = i
1
1 + k + k(k − 1) + · · · + k(k − 1)d−2 c + k(k − 1)d−1 c = 1,
k which indicates that h
i−1
c = 1 + k + k(k − 1) + · · · + k(k − 1)d−2 + (k − 1)d−1
. Thus the chain is reversible, since we have found a solution to the detailed balance conditions.
Remark: If you only solve it for k = 4, d = 3, you will get full credits.
(Expected failure via MCMC) Consider a video streaming service with a faulty
server. Let X denote the availability of the service (1 if available, 0 if server is down), which we
model as a Bernoulli distribution with probability θ. We incorporate our prior knowledge about
server failures by viewing θ as a random variable and assigning a prior distribution on θ:
P(θ = θ ) = p , for j = 1, . . . , m.
(we assume p > 0 for all j = 1, . . . , m)
Suppose that we have observed server availabilities X = x , . . . , X = x .
(a) Compute the posterior pmf P(θ = θ |X = x ) in terms of θ , p and X , . . . , X 's.
We now wish to see how Markov chain Monte Carlo (MCMC) can approximate samples from the
posterior distribution computed in part (a).
(b) Consider the Metropolis-Hastings algorithm with some proposal distribution Q = (Q(z , z ) :
z , z ∈ S) where S = {θ , . . . , θ }. Construct a Markov chain Z = (Z : n ≥ 0) with state
space S and stationary distribution π(θ ) = P(θ = θ |X = x ). Argue that
lim P(Z = θ |Z = z) = π(θ ) for any z ∈ S
(1)
if Q(x, y) > 0 for all x, y ∈ S.
Take m = 4 with θ = .2, θ = .4, θ = .6, θ = .8 and p = .4, p = .4, p = .1, p = .1. Suppose
that we have observed server availabilities
Question 3.3 : j j j 1 n
1 j 1 n
1 n j n j 1 n 1 1 2 1 m j n n→∞ 1 2 j 3 2 n n
1 j 0 4 n
1 j 1 2 3 4 X1 = 1, X2 = 1, X3 = 0, X4 = 1, X5 = 0,
X6 = 1, X7 = 0, X8 = 1, X9 = 0, X10 = 0. (c) Implement a Metropolis-Hastings algorithm and run it up to period n = 100, 000 for the uniform
proposal distribution Q(x, y) = . Plot a histogram of the Markov chain path from time
n = 50, 000 to n = 100, 000 alongside the posterior density computed in part (a). What do
you see?
(d) Generate the same plot for (n , n ) = (50, 100), (500, 1000), (5000, 10000). What is your suggested run length for the procedure in part (c) based on these plots?
4
1
m 1 2 1 2 Copyright c Peter Glynn All rights reserved. MS&E 221 Spring 2020 Answer: (a) From Bayes rule and since X , . . . , X are i.i.d. conditional on θ, we have
1 π(θk ) := P(θ = θk |X1n = xn1 ) n Q
P(θ = θk ) ni=1 P(Xi = xi |θ = θk )
P(X1n = xn1 |θ = θk )P(θ = θk )
Qn
P
= m
=
P(X1n = xn1 )
i=1 P(Xi = xi |θ = θj )
j=1 P(θ = θj )
Pn pk θk =P
m i=1 n Pn (1 − θk )n− Pn j=1 pj θj (b) Consider the Markov chain Z = (Z Xi i=1 : n ≥ 0) Xi i=1 (1 − θj )n− Xi Pn i=1 Xi . (2) with transition matrix given by
π(z2 )Q(z2 , z1 )
P (z1 , z2 ) = Q(z1 , z2 ) min 1,
π(z1 )Q(z1 , z2 ) where π(θ ) is as dened in (2). The Markov chain Z satises detailed balance with the
distribution π so that its stationary distribution is indeed given by π as claimed. In particular,
Z is positive recurrent.
Note that π(θ ) > 0 for all k = 1, . . . , m. Since Q is irreducible with Q(z, z) > 0 for some
z ∈ S , Z is also irreducible with P (z, z) > 0. Hence, Z is an irreducible, positive recurrent, and
aperiodic, we have the limit (1) holds.
(c) See the attached code and the gure. The histogram of the sample path of the chain from time
n = 50, 000 to n = 100, 000 more or less matches the posterior.
(d) Judging from the gure, we recommend using run size at least 1000.
k k 5 Copyright c Peter Glynn All rights reserved. MS&E 221 Spring 2020 50000 to 100000 50 to 100
0.8 histogram
posterior 0.7
0.6 0.6
probability probability 0.5
0.4
0.3 0.4
0.3
0.2 0.1 0.1
0.0
0.2 0.3 0.4 0.5
theta 0.6 0.7 0.8 0.2 0.3 0.4 500 to 1000 0.8 0.5
theta 0.6 0.7 0.8 5000 to 10000
histogram
posterior 0.7 histogram
posterior 0.7
0.6 0.6 0.5 0.5 probability probability 0.5 0.2 0.0 histogram
posterior 0.7 0.4
0.3 0.4
0.3 0.2 0.2 0.1 0.1 0.0
0.2 0.3 0.4 0.5
theta 0.6 0.7 0.8 0.0 0.2 0.3 0.4 0.5
theta 0.6 0.7 0.8 Figure 2: Posterior vs Histogram Julia Code
using PyPlot;
function target(z, z_prop, p, p_prop, num_succ, num_fail)
return (p_prop / p) * (z_prop / z)^num_succ
* ((1.-z_prop) / (1.-z))^num_fail;
end
X = [1., 1., 0., 1., 0., 1., 0., 1., 0., 0.];
thetas = [.2, .4, .6, .8];
# ps = [.3, .1, .4, .2];
ps = [.4, .4, .1, .1];
# thetas = [.1, .2, .3, .4, .5, .6, .7, .8, .9];
# ps = [.1, .2, .05, .1, .05, .2, .15, .1, .05];
mm = length(thetas);
num_succ = sum(X);
num_fail = length(X) - num_succ; 6 Copyright c Peter Glynn All rights reserved. MS&E 221 Spring 2020 n_iter = 100000;
z = thetas[4];
p = ps[1];
mcpath = zeros(n_iter+1);
mcpath[1] = z;
accept = 0;
for ii = 1:n_iter
# Sample from uniform proposal
# uniformly sample a number between 1, .., 4
a = rand(1:mm, 1)[1];
z_prop = thetas[a];
p_prop = ps[a];
t = min(1.0, target(z, z_prop, p, p_prop, num_succ, num_fail));
# flip coin
u = rand(1)[1];
if u < t
z = z_prop
accept += 1;
end
mcpath[ii+1] = z;
end
# Compute posterior
normalizer = 0;
for ii = 1:mm
normalizer += ps[ii] * (thetas[ii])^num_succ * (1.-thetas[ii])^num_fail;
end
post = ps .* thetas.^num_succ .* (1 .- thetas).^num_fail / normalizer; # Plot posterior and histogram together
plot_start = 50000;
hist = zeros(mm);
for jj = 1:mm
hist[jj] = sum(mcpath[plot_start+1:end] .== thetas[jj])
/ (n_iter - plot_start + 1);
end
figure(1); 7 Copyright c Peter Glynn All rights reserved. MS&E 221 Spring 2020 ax = gca();
title(string(plot_start, " to ", n_iter));
plot(thetas, hist, color = "r", label = "histogram");
plot(thetas, post, color = "k", label = "posterior");
xlabel("theta");
ylabel("probability");
legend(["histogram", "posterior"]);
savefig("./figures/long-run-mcmc.pdf");
close("all");
plot_start = 50;
plot_end = 100;
hist = zeros(mm);
for jj = 1:mm
hist[jj] = sum(mcpath[plot_start+1:plot_end] .== thetas[jj])
/ (plot_end - plot_start);
end
figure(1);
ax = gca();
title(string(plot_start, " to ", plot_end));
plot(thetas, hist, color = "r", label = "histogram");
plot(thetas, post, color = "k", label = "posterior");
xlabel("theta");
ylabel("probability");
legend(["histogram", "posterior"]);
savefig(string("./figures/short-", plot_start, "-mcmc.pdf"));
close("all"); plot_start = 500;
plot_end = 1000;
hist = zeros(mm);
for jj = 1:mm
hist[jj] = sum(mcpath[plot_start+1:plot_end] .== thetas[jj])
/ (plot_end - plot_start);
end
figure(1);
ax = gca();
title(string(plot_start, " to ", plot_end));
plot(thetas, hist, color = "r", label = "histogram");
plot(thetas, post, color = "k", label = "posterior");
xlabel("theta");
ylabel("probability");
legend(["histogram", "posterior"]); 8 Copyright c Peter Glynn All rights reserved. MS&E 221 Spring 2020 savefig(string("./figures/short-", plot_start, "-mcmc.pdf"));
close("all"); plot_start = 5000;
plot_end = 10000;
hist = zeros(mm);
for jj = 1:mm
hist[jj] = sum(mcpath[plot_start+1:plot_end] .== thetas[jj])
/ (plot_end - plot_start);
end
figure(1);
ax = gca();
title(string(plot_start, " to ", plot_end));
plot(thetas, hist, color = "r", label = "histogram");
plot(thetas, post, color = "k", label = "posterior");
xlabel("theta");
ylabel("probability");
legend(["histogram", "posterior"]);
savefig(string("./figures/short-", plot_start, "-mcmc.pdf"));
close("all"); (Epidemic Modeling) Consider an epidemic spreading through the Stanford campus. The disease is infectious for only one period of time. Let X denote the number of newly
infected individuals at time n and let Y denote the number of people infected in the next period
by the i-th infected individual at time n. Assume that Y 's are i.i.d. random variables taking
values in {0, 1, 2, . . .} for all i ≥ 0 and n ≥ 0. Suppose there is one infected individual at time 0.
(a) Write down an expression for X in terms of Y 's for n ≥ 1. Is X = (X : n ≥ 0) a Markov
chain?
(b) If E[Y ] = µ ∈ [0, ∞) for all i ≥ 0 and n ≥ 0, derive an expression for the expected number of
infected individuals at time n in terms of µ for n ≥ 1. Using this answer, derive an expression
for the expected total number of infected individuals at time n in terms of µ for n ≥ 1.
(c) (Bonus) Now, assume that for all i ≥ 0, n ≥ 0, Y ∼ Pois(λ) for λ < 1. Derive an expression
for the expected time until the epidemic dies outmeaning there are no remaining newly
infected individuals to infect other peoplein terms of the probability generating function
G(z) := E[z ] = e
for the Poisson distribution Y ∼ Pois(λ).
Hint: Show that
G (z) := E [z ] = G
· · ◦ G}(z)
| ◦ ·{z
times
and use P (X = 0) = G (z).
Question 3.4 : n n, i n,i n n−1,i n n, i i.i.d. n,i Y λ(z−1) n 1 Xn n 1 n n 9 Copyright c Peter Glynn All rights reserved. MS&E 221 Spring 2020 Answer: (a) We have Xn−1 Xn = X Yn−1, i i=1 for all n ≥ 1. The process X = (X : n ≥ 0) is a branching process; as covered in class, it is a
Markov chain.
(b) We begin by noting that X is independent of Y for all n ≥ 0 and i ≥ 1. Hence
n n n, i Xn−1 X Yn−1, i E1 [Xn ] = E i=1 = ∞
X Xn−1 X E1 x=0 i=1 ∞
X " Yn−1, i Xn−1 = x P1 (Xn−1 = x) x
X # =
E1
Yn−1, i Xn−1 = x P1 (Xn−1 = x)
x=0
i=1
" x
#
∞
X
X
=
E1
Yn−1, i P1 (Xn−1 = x)
x=0 = ∞
X i=1 xµP1 (Xn−1 = x) x=0 = µE[Xn−1 ]. Here, the fourth equality follows from independence and the fth equality follows from the fact
that the Y 's are i.i.d. with mean µ. Since X = 1, we conclude that E [X ] = µ .
By linearity of expectation, we then have
#
"
(
X
X
X
if µ 6= 0
E
X =
E [X ] =
µ =
n+1
if µ = 1
(c) Dene T := min(n ≥ 0 : X = 0). We wish to compute E [T ]. First, we note that since T = n
happens if and only if X = 0 and X > 0, we have P(T = n) = P(X = 0) − P(X = 0).
Denoting by G (z) = E [z ] the probability generating function (PGF) for X , we have
P(T = n) = G (0) − G
(0). Now, recall the recursive formula for G (z)
j, i 0 n n 1 n 1 k k=0 1 k=0 k k 1 n n µn+1 −1
µ−1 k=0 n 1 n n n n−1 n Xn n−1 n n−1 n Gn (z) = G
| ◦G
{z· · · G}(z) where G(z) = E[z
We conclude that Y1,1 ] = expλ(z−1) E1 [T ] = ∞
X
n=1 . n nP(T = n) = times ∞
X
n=1 n(Gn (0) − Gn−1 (0)). (3)
10 Copyright c Peter Glynn All rights reserved. Recall that for Y ∼ Geom(p) with
, the probability generating function is given by G(s) = E[z ] =
. Solution for geometric distributions P(Y = k) = Let MS&E 221 Spring 2020 Yn,i ∼ Geom(p) (1 − p)k p Y
0
p
A=
,
−(1 − p) 1 Using induction, one can show that Gn (s) = p
1−s(1−p)
an bn
A =:
.
cn dn
n an s + bn
.
cn s + dn (4) For example, if p = , then we can show inductively that
1
2
1 −(n − 1)
n
.
A = n
−n
n+1
2
n Plugging in the preceeding display in (4) and in turn in the expression (3), we arrive at an
expression for E[T ]. 11 ...

View
Full Document