03-ex-solution 2020.pdf - Copyright c Peter Glynn All...

This preview shows page 1 out of 11 pages.

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

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture