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.
 A: 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[1], sd=sig[1]^2))
  logPYgminus = log(pnorm(q=Y[i], mean=mu[2], sd=sig[2]^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.)
A: Answer to the follow-up question
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.
