# Gibbs Sampling and Probability Notation

Problem 1

I am trying to implement Gibbs Sampling for the following problem:

There is a grid measuring 3 x 3 sites, each "site" can be designated in a state, $X$, of 1 or -1. The sites are numbered 1--9 and have corresponding observed values, $Y$. Similar to the Ising model, we can define the probability of the states via $$\mathbb{P}(X) \propto \prod_{(i,j) \in G} \exp{(J \cdot X_i X_j)}$$. Additionally, "the data $Y$ are conditionally independent given $X$", $$Y_i|X_i = s \sim N(\mu_s, \sigma^2_s)$$ where $i$ is the site index, and $s$ is the state. We are given $\mu_1$, $\sigma_1$, $\mu_{-1}$, $\sigma_{-1}$.

I need to implement a sampler that has $\mathbb{P}(X|Y)$ as target. I know for the Gibbs sampling algorithm we need the complete conditional probabilities, given by $$\mathbb{P}(X_i|X_{[-i]},Y) \propto \exp{\left\lbrace X_i \sum_{j \in N_i} X_j + \log \mathbb{P}(Y_i|X_i) \right\rbrace}$$ where $N_i$ are the neighbors of the $i$th site.

Solution 1

I have begun implementing the algorithm but am confused regarding some of the probability notation, specifically $\log \mathbb{P}(Y_i|X_i)$. To calculate this value, do I simply evaluate as $$\log \frac{1}{\sigma \sqrt{2 \pi}} \exp{\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)}$$ where $x$ = $Y_i$, $\sigma = \sigma_{X_i}$, $\mu = \mu_{X_i}$, so for example, if $X_1$ = 1, $Y_1$ = 2, $\mu_1 = 2$, $\sigma_1 = 1$, $$\log \frac{1}{\sqrt{2 \pi}} \exp{\left(-\frac{(2 - 2)^2}{2}\right)}$$

My question is, when I calculate $\mathbb{P}(X_i|X_{[-i]},Y)$, I get non-integer values. However, $X$ must be an integer of 1 or -1 (to represent the states), so is there another step after calculating $\mathbb{P}(X_i|X_{[-i]},Y)$ to get an integer representation, such as a Bernoulli distribution?

If it helps, here is the R code I'm working, which includes the actual values for all parameters.

site = c(1, 2, 3, 4, 5, 6, 7, 8, 9)  # site position
Y = c(2, 2, 2, 2, 0, 0, 1, 2, 1)  #  vegetation indices

N = list(c(2, 4), c(1, 3, 5), c(2, 6),
c(1, 5, 7), c(2, 4, 6, 8), c(3, 5, 9),
c(4, 8), c(5, 7, 9), c(6, 8))  # neighbor list

mu = c(2, 0.5)
sig = c(1, 0.5)

gibbsJ = function(J,n) {
X = matrix(1, n, 9)
for (t in seq(2,n)) {
X[t,] = X[t-1,]
for (i in seq(9)) {
logPYgX = log(pnorm(q=Y[i], mean=mu[X[t,i]], sd=sig[X[t,i]]^2))
PXi = exp(J*X[t,i]*sum(X[t,N[[i]]]) + logPYgX)
if (PXi < runif(1)){
X[t,i] = -X[t,i]
}
}
}
return(list(X=X, p=colMeans(X)))
}

soln = gibbsJ(J=0.2, n=1000)
soln$X; soln$p


Followup Question

I need to calculate the log conditional density (up to a normalizing contant) $$f(X^t|Y) = \sum_{(i,j) \in G} JX_i^tX_j^t + \log\mathbb{P}(Y|X^t)$$ using results from a simulation of the system. To do the above calculation, my instinct is to use the a formula similar to that used for the complete conditional probability: $$f(X^t|Y) = \prod_{i=1}^9 JX_i\sum_{j \in N_i} X_j^t + \log\mathbb{P}(Y_j|X^t_j)$$ or in R

X = gibbsJ(J=0.2, n=1000)
J = 0.2
fXtgY = numeric(1000)
for (t in seq(1000)) {
logPYgX = log(pnorm(q=Y[i], mean=mu[X[i]], sd=sig[X[i]]^2))
fXtgY[t] = fXtgY[t] + J*X[t,i]*sum(X[N[[i]]]) + logPYgX
for (i in seq(2,9)) {
logPYgX = log(pnorm(q=Y[i], mean=mu[X[i]], sd=sig[X[i]]^2))
fXtgY[t] = fXtgY[t]*(J*X[t,i]*sum(X[N[[i]]]) + logPYgX)
}
}


I have the product there because it is supposed to be the log conditional likelihood for a given time step, and I figure the likelihood of a given configuration is the product of the individual likelihoods. Is the approach/method correct?

I am also looking for an MCMC estimate for $\mathbb{P}(X_i|Y)$ and the entropy (modulo a constant) of the system, $$\sum_X f(X|Y)\mathbb{P}(X|Y)$$ For the former, I believe the MCMC estimates of $\mathbb{P}(X_i|Y)$ is simply the probability of the posterior distribution of the simulation and can be calculated by averaging the chains over the sites, (i.e., colMeans(X)). For the latter, I am less confident. Based on the definition, this should return a scalar; however, from the calculations about $f(X|Y)$ is a vector over the time series, and $\mathbb{P}(X|Y)$ is a vector over the sites. What am I missing?

It seems my main problem is my misunderstanding of the probability notation and confusion about the indices (or lack of) which would help describe what I'm looking for.

Very clean and clear presentation of the issue! This is exactly a hidden Ising model.

Your probabilities $\mathbb{P}(X_i=1|X_{[-i]},Y_i)$ and $\mathbb{P}(X_i=-1|X_{[-i]},Y_i)$ give you the way to simulate $X_i$ conditional on the others and on the observations. That those probabilities are non-integer is not an issue.

Maybe what is confusing you is the proportional $\propto$ symbol in $$\mathbb{P}(X_i|X_{[-i]},Y) \propto \exp{\left\lbrace X_i \sum_{j \in N_i} X_j + \log \mathbb{P}(Y_i|X_i) \right\rbrace}$$ Indeed, without the $\propto$ symbol , we have \begin{align*} \mathbb{P}(X_i=1|X_{[-i]},Y) = \exp{\left\lbrace 1 \sum_{j \in N_i} X_j + \log \mathbb{P}(Y_i|1) \right\rbrace} \Big/ \\ \left[\exp{\left\lbrace 1 \sum_{j \in N_i} X_j + \log \mathbb{P}(Y_i|1) \right\rbrace}+ \exp{\left\lbrace -1 \sum_{j \in N_i} X_j + \log \mathbb{P}(Y_i|-1) \right\rbrace}\right] \end{align*} In your R code, you therefore have to compute both values:

for (i in seq(9)) {
logPYgplus = log(pnorm(q=Y[i], mean=mu, sd=sig^2))
logPYgminus = log(pnorm(q=Y[i], mean=mu, sd=sig^2))
PXiplus = exp(J*sum(X[t,N[[i]]]) + logPYgplus)
PXiminus = exp(-J*sum(X[t,N[[i]]]) + logPYgminus)
X[t,i]=1
if (PXiplus < runif(1)*(PXiplus+PXiminus)){
X[t,i] = -1}
}


where I assumed the index 1 corresponds to $X_i=1$ and the index 2 to $X_i=-1$ in your definition of the vectors $\mu$ and $\sigma$. (It is a matter of convention.)

• You are spot on. I did not understand how to deal with the proportional sign. It also didn't help that I was returning the solution prematurely. After implementing your (incredibly useful help) with the change to return(), everything is running great. I won't be in France anytime in the near future, but I believe I owe you a bottle of wine. Thanks again, Xi'an! Nov 21, 2014 at 7:21
• I updated my question with followup question(s), any additional help you can provide would be appreciated. Nov 21, 2014 at 22:28

Unfortunately, your intuition about the smoothing distribution $\mathbb{P}(x|y)$ is not correct as $\mathbb{P}(x|y)$ is not the product of the full conditionals: $$\prod_{i=1}^9 \mathbb{P}(x_i|x_{[-i]},y) \ne \mathbb{P}(x|y)$$ even up to a proportionality constant. And even less on a log scale. Note also that the time step $t$ should not appear in the formula as time is only connected with Gibbs iterations in your problem, not with the initial statistical problem.
Now, you know $f(x|y)$ up to a constant. This constant is actually $$\sum_{x\in\{-1,1\}^9} \mathbb{P}(x|y)$$ which remains a manageable sum, hence can be calculated.
My first solution for the marginal posterior density of $X_i$ (given $Y$) would be to use the decomposition \begin{align*} \mathbb{P}(x_i|y) &= \sum_{x_{[-i]}} \mathbb{P}((x_i,x_{[-i]})|y)\\ &= \sum_{x_{[-i]}} \mathbb{P}(x_i|x_{[-i]},y) \mathbb{P}(x_{[-i]}|y)\\ &= \mathbb{E}\left[ \mathbb{P}(x_i|X_{[-i]},y) | y \right] \end{align*} where the expectation is understood as the one of $X_{[-i]}$ given $y$. Since the Gibbs sampler returns simulations from all subsets of $X$ given $y$, a converging approximation to the above is $$\hat{\mathbb{P}}(x_i|y) = \frac{1}{T} \sum_{t=1}^T \mathbb{P}(x_i|x_{[-i]}^t,y)$$ where $x^t$ denotes the value of the Markov chain at the $t$-th iteration.
However, since your state space $\{,-1,1\}^9$ is small, you can also derive this marginal exactly: $$\mathbb{P}(x_i|y) = \sum_{x_{[-i]}\in\{-1,1\}^8} \mathbb{P}((x_i,x_{[-i]})|y) .$$ Hence an exact derivation of the entropy as well.