# Question on Gibbs sampler, estimate the density function

I am reading the paper Explaining the Gibbs Sampler(Casella, George, and Edward I. George. "Explaining the Gibbs sampler." The American Statistician 46.3 (1992): 167-174). I am stuck with example 2, Figure 2, as how the solid line was created from the equation $$(2.9)$$.

http://www2.stat.duke.edu/~scs/Courses/Stat376/Papers/Basic/CasellaGeorge1992.pdf

The figure is here :

The paper says' Gibbs sampling can be used to estimate the density itself by averaging the final conditional densities from each Gibbs Sequence.

I know how to create the final Gibbs sequence for example 2. Here is an excellent example of how to create these sequences (Gibbs sampler for conditionals that are exponential: Example from Casella & George paper).

I know I can use the Kernel density method to estimate the density function using the 500 final sample points. However, this was not the method the paper used.

Here is the code to produce 500 final sample points, which I copied from the above URL.

set.seed(1010)
X = rep(0, 500)
Y = rep(0, 500)
k = 15

for (i in 1:500) {
x = rep(1, k)
y = rep(1, k)

for (j in 2:k) {
temp_x = 6
while(temp_x>5) {
x[j] = rexp(1,y[j-1])
temp_x = x[j]
}

temp_y = 6
while(temp_y>5) {
y[j] = rexp(1,x[j])
temp_y = y[j]
}
}

X[i] = x[k]
Y[i] = y[k]
}


From these values, we can calculate the density by equation $$(2.8a)$$.

In the equation $$(2.9)$$ we have a fixed $$x$$ value, such as if we let $$x=\xi$$ then equation $$(2.9)$$ becomes

$$\hat{f}(\xi)=\frac{1}{m}\left[f(\xi|y_1)+f(\xi|y_2)+f(\xi|y_3)+...+f(\xi|y_m)\right]$$

I think the Gibbs Sampler did not produce a fixed last value (15th) $$\xi$$ for the 500 sequences. I am stuck here. I hope you can help me.

Anyway, overall the question is how the solid line in figure 2 was produced.

• This question is difficult to follow without Figure 2. Do you think you could reproduce it or create some kind of facsimile?
– whuber
Commented Mar 25, 2023 at 13:26
• Thanks @ whuber, now I add a link for the paper and the figure 2. Commented Mar 26, 2023 at 3:05
• I suggest you erase the 500 values from the post as it is not helping. Commented Mar 27, 2023 at 10:26

The kernel density estimator $$\hat f(\cdot)$$ is a generic density estimation method that produces a density function out of a sequence of generations from the true distribution $$f$$. It is biased and converging (for the univariate case) at a speed of $$n^{-2/5}$$, in general.

The Rao-Blackwell(ised) estimator (2.9) is unbiased, meaning that for every entry $$x$$ $$f(x) = \mathbb E[\hat f(x)] = \mathbb E^Y[f(x|Y)],$$ and it converges at the parametric speed of $$n^{-1/2}$$. This is a gain exploiting the latent variable simulated by the Gibbs sampler.

Note that the above R code is not perfect, as the 500 chains are all starting from the value $$X_1=1$$, rather than exploiting the higher proximity of the previous chain with the stationary distribution at the end of the 15 iterations. Here is my alternate version with no parallel chain and a thinning of one out of 15:

X = Y= runif(500,0,5)
y=Y[1]
for (t in 1:500) {
for (j in 2:15) {
temp_x = temp_y = 6
while(temp_x>5)temp_x = x = rexp(1,y)
while(temp_y>5)temp_y = y = rexp(1,x)
}
X[t]=x;Y[t]=y}

plot(density(X,from=0,to=5),main="")
rb=seq(0,5,le=123)
for(i in 1:123)
rb[i]=mean(Y*exp(-rb[i]*Y))
lines(seq(0,5,le=123),rb,col="red2",lwd=2)


Here is a comparison of the density estimates based on the Gibbs chain $$(X_t)$$ using a standard non-parametric approach (through R density) [in black] and of the Rao-Blackwell estimate based on the Gibbs chain $$(Y_t)$$ [in red]:

The estimators of the densities are close enough to conclude that the Rao-Blackwell version is not completely off.

• Thank you very much @Xi'an. I will read the Rao-Blackwellised estimator and see if I can understand, thanks again. Commented Mar 26, 2023 at 3:08
• I am reading your book Introducing Monte Carlo Methods with R, chapter 7, especially page 229, I still cannot figure out how the marginal density can be estimated from the Gibbs sampler output (the 500 sample points) by the method, especially in terms of the fixed value $x$. Thank you for the help. Commented Mar 26, 2023 at 7:23
• This is a consequence of the Law of Total Expectation. It applies for all values of $x$, so I do not understand the difficulty you have with this formula. Commented Mar 26, 2023 at 9:28
• Thank you very much. I mean suppose, x=1 for the equation $(2.9)$ then I need to calculate $f(1)=\frac{1}{500}[f(1|y_1)+f(1|y_2)+....+f(1|y_m)]$, how can I make sure that the last value of the 15 long sequences is 1 for x? then what are $y_1,y_2,...y_m$?. It seems I can understand the formula, but I don't understand, how a Gibbs sampling can implement this. Commented Mar 26, 2023 at 10:00
• The 15th realisation of the Markov chain $(X_t)$ does not matter. The argument is that the chain $(Y_t)$ has reached stationarity and is thus distributed from the right marginal, meaning that the unbiasedness property holds for all $x$'s, not only for the corresponding realisation of $X_{15}$. Commented Mar 26, 2023 at 14:07