ReferencesCS/CNS/EE 155: Probabilistic Graphical ModelsProblem Set 3Handed out: 9 Nov 2009Due: 23 Nov 2009Gibbs samplingPreamble. Suppose that we are given a full joint distribution p(Y1, . . . , Yn) over randomvariables, but are only interested in drawing samples Yi∼ p(Yi). The natural solution to thisproblem would be to actually calculate the marginal by integrating over {Y1, . . . Yn} − Yi:p(Yi) =X{Y1,...Yn}−Yip(Y1, . . . , Yn),from which samples Yimight be drawn.Unfortunately, in many cases (e.g., where the joint distribution has high treewidth), the re-quired marginalization is intractable. In Homework 1, you showed that if p(Y1, . . . , Yn) isgiven as a graphical model, and Y1, . . . , Ynare all discrete, and you have observations of allvariables except Yk, computing the conditional distributionp(Yk| Y1= y1, . . . , Yk−1= yk−1, Yk+1= yk+1, . . . , Yn= yn)is very efficient. In this case, sampling from this conditional distribution only requires sam-pling from a Bernoulli random variable. Gibbs sampling is an approximate inference methodpeculiarly well suited to such situation, where one can sample from the conditional distribu-tionsp(Yk|{Y1, . . . Yn} − Yk) ∀k.As outlined in both [KF09, Bis06], the algorithm for Gibbs sampling is quite straightforward:1. Generate an initial guess (t = 0) about the state (Y1= yt=01, . . . , Yn= yt=0n).2. Repeat, over t, until convergence (a very loaded statement—see mixing times, stationarydistributions, etc):Draw yt+11∼ p(Y1|Y2= yt2, . . . Yn= ytn)Draw yt+12∼ p(Y2|Y1= yt+11, Y3= yt3. . . Yn= ytn)...Draw yt+1n∼ p(Yn|Y1= yt+11, Y2= yt+12. . . Yn−1= yt+1n−1)As it turns out, these samples {y1i, y2i, . . . , yTi} not only converge (in the sense that given a longenough burn in time, the samples are all coming from the same distribution), but converge tothe “correct” distribution, which in this case is the marginal p(Yi).1Problem.1. Write a Gibbs sampler to estimate the marginal distribution generated from the condi-tional distributionsp(x|y) ∝ ye−yx, 0 < x < B < ∞p(y|x) ∝ xe−xy, 0 < y < B < ∞where B is a known positive constant (as we will see later, the restriction to the intervalx, y ∈ (0, B) ensures that the marginal p(x) exists).For B = 5, and for sample sizes T = 500, 5000, 50000, plot the histogram of values for x.Note that although the simulation of p(x) is straightforward using a Gibbs sampler, themarginal is not obvious from inspection of the above marginals.2. Importantly, simulation based approximations provide straightforward methods of calcu-lating the moments of the marginal distributions. Provide an estimate of the expectationof X, Ep(X)[X] by using the 500, 5000 and 50000 samples from 1).3. [Extra Credit:] Consider the bivariate case, where we are given pX|Y(X|Y ) and pY |X(Y |X).As it turns out, we can derive a fixed point integral equation to determine the marginaldensity pX(X), and thus the joint density, using repeated applications of the chain ruleand marginalization: first, notice thatpX(X) =ZpXY(X , Y )dY =ZpX|Y(X|Y )pY(Y ) dY.However, we can similarly substitute for pY(Y ):pX(X) =ZpX|Y(X|Y )ZpY |X(Y |Z)pX(Z) dZ dY (Z is a dummy variable)=ZZpX|Y(X|Y )pY |X(Y |Z) dYpX(Z) dZ=Zh(X , Z)pX(Z) dZ.Using the normalized conditionals from problem 1p(x|y) =ye−yx1 − e−By, 0 < x < B < ∞p(y|x) =xe−xy1 − e−Bx, 0 < y < B < ∞calculate px(x). Hint: px(x) ∝ (1 − e−Bx)/x – find the normalizer and verify that pxsolves the fixed point equations.4. [Extra Credit:] Now that we have an analytic form for the marginal px(x) from problem3, we can compare with our earlier estimates.(a) Plot px(x) and the histogram of collected samples of px(x).(b) Compute Epx(x)[x] using the analytic form found in problem 3.2References[Bis06] C. Bishop. Pattern Recognition and Machine Learning. Springer Science+BusinessMedia, LLC, New York, NY, 2006.[KF09] D. Koller and N. Friedman. Probabilistic Graphical Models: Principles and Techniques.The MIT Press, Cambridge, MA,