05-ex.pdf - Problem Set 5 MS&E 221 Due Friday:59PM Note on...

This preview shows page 1 out of 4 pages.

Unformatted text preview: Problem Set 5 MS&E 221 Due: Friday, February 22 11:59PM Note on Problem 5.1 and 5.2: You can do either Problem 5.1 or Problem 5.2. problems have equal credits. These two Remember to correctly tag the one you select when submitting on Gradescope (leaving the other one blank). Note on Problem 5.4: Part (b) amd (c) use material we will cover on Tuesday, Feb 19. Question 5.1 (Predicting Hurricanes via MCMC): NB(r, p) with parameter r > 0 (r is an integer) and The negative binomial random variable p ∈ (0, 1) H∼ has the following probability mass function:   k+r−1 P (H = k) = (1 − p)r pk , k = 0, 1, 2, . . . k It can be seen as the number of successes in a sequence of i.i.d. Bernoulli(p) trials before r failures occur. The negative binomial distribution also arises as a continuous mixture of Poisson distributions. More precisely (you do not have to show this), a NB(r, p) r.v. Poisson(Λ) r.v., where the parameter Λ follows a Gamma(r, p/(1 has the same distribution as a − p)) distribution. Due to this extra variability, the negative binomial model is sometimes more suitable than Poisson distribution in real world applications. More precisely, when the data has a larger variance than the mean (often referred to as overdispersion), the Poisson distribution cannot capture this overdispersion since Poisson(λ) (with a xed In all the problems below, we consider λ) has mean and H ∼ NB(r, 1/2). variance both equal to λ. We have observed the number of major hurricanes in the East Pacic Basin from the past 14 1 years : H1 = h1 , . . . , H14 = h14 . We incorporate our prior knowledge about hurricane behavior by assigning a prior distribution on r: P(r = rj ) = qj , (we assume qj > 0 for all for j = 1, . . . , m. j = 1, . . . , m) (a) Compute the posterior pmf P(r = rj |H1 = h1 , . . . , H14 = h14 ) in terms of rj , qj and h1 , . . . , h14 . We now wish to see how Markov chain Monte Carlo (MCMC) can approximate samples from the posterior distribution computed in part (a). Take 0.1. m=5 The r1 = 3, r2 = 4, r3 = 5, r4 = 6, r5 = 7 observed values h1 , . . . , h14 are with and q1 = .1, q2 = .3, q3 = .3, q4 = .2, q5 = 2, 6, 1, 2, 5, 2, 6, 5, 1, 9, 11, 6, 4, 10 1 Data resource: 1 (b) Implement a Metropolis-Hastings algorithm and run it up to time proposal distribution Q(x, y) = 1 m. n = 100, 000 Plot a histogram of the Markov chain's empirical state frequencies (i.e. the proportion of time spent in each state up to time and time n2 = 100, 000 for the uniform ni ) for time n1 = 25, 000 alongside the density computed in part (a). What do you see? (c) Use the method of batch means to produce a condence interval for the posterior probability of having 3 or fewer hurricanes in a year. Question 5.2 (Ising model): Consider the Ising model that we discussed in class, in which one has a two-dimensional Ising model consisting of nine sites in total (3 by 3). The corners each have two neighbors, the center site has 4 neighbors, and the sites in the middle of each boundary have 3 neighbors. Assume that if site are not neighbors, then Ji,j = 0. i and site j are neighbors, then We also assume that hj = 0 Ji,j = a, whereas if j . Assume we for all sites i and j start from a conguration where the three sites on the top edge have spin equal to 1, while all the other sites have spin equal to -1. (a) Run the Gibbs sampler for 1,000,000 time steps with a=1 to estimate the probability that 8 or more sites will have the same spin. Use the method of batch means to provide a condence interval. (b) Repeat part (a) with a = 4. Question 5.3 (Cayley Tree): The Cayley tree of degree k and depth d is a tree with undirected d−1 nodes, where every non-leaf node has edges, on 1 + k + k(k − 1) + · · · + k(k − 1) (Assume k>2 and k adjacent edges. 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 note that any non-leaf node has exactly of v in the graph; i.e., p(v) k neighbors. For a leaf is the unique node adjacent to 2 v. N (i) denote the neighbors of i; node v , let p(v) denote the parent Figure 1: A Cayley tree with k=4 d = 3; note that it has 1 + k + k(k − 1) + k(k − 1)2 = 53 k = 4 adjacent edges. and nodes, and that each non-leaf node has (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. Question 5.4 (Inventory Shortage): Let Xn Consider the (s, S)-inventory be the inventory position at the end of each day. model with (we assume Zn+1 's are independent from i.i.d. Zi ∼ and S = 7. We assume that the delivery arrives immediately at the end of each day. Suppose that the daily demands ∗ distribution with unknown parameter p so that s=3 Zn+1 follow the geometric P (Z1 = k) = (1 − p∗ )k p∗ ∗ Geom(p ) where Xn 's). Suppose we observe the inventory X0 = 7, X1 = 7, X2 = 7, X3 = 7, X4 = 7, X5 = 6, X6 = 4, X7 = 7, X8 = 5, X9 = 5, X10 = 5, X11 = 5, X12 = 7, X13 = 6, X14 = 6, X15 = 7. (a) Compute the MLE for the parameter (b) Varying the number of bootstraps p∗ . b = 10, 50, 100, 500, 1000, compute a 90% condence interval b? ∗ for p and plot it. How does the width of the condence interval change with We wish to compute the expected long-run inventory shortage cost n−1 1X [Zi+1 − Xi ]+ n→∞ n c = 10 · lim i=0 3 c. The quantity c is dened as This will allow us decide whether we should change our inventory strategy. (If the expected long-run cost c is too large, then to mitigate the eects of inventory shortage, we will increase more frequently and build a larger warehouse to have a larger capacity (c) Give a 90% condence interval for c histogram and provide your values of using the parametric bootstrap with z1 and s to replenish S ). b = 100. Plot the z2 . In MATLAB, you can use the functions fminunc, fmincon and fminsearch to optimize. In Julia, you may use the Optim.jl package to optimize. In Python, you may use minimize in package scipy.optimize to optimize. Hint: To implement the likelihood function, you may dene a function that takes both the observed data and the parameter as arguments. You should then do a partial evaluation with respect to the data, and pass the function of merely the unknown parameter to the optimizer. To do this, in MATLAB you can use the function handle @, and in Python you can use partial in the functools package. Julia handles the partial evaluation elegantly so there is nothing special you need to do. 4 ...
View Full Document

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture