# How to Gibbs sample proportional to a probability

I am reading this tutorial on Hierarchical Chinese Restaurant Process. On pdf page 141 (slide title: MCMC Problem Specification for N-grams) it says:

$$F(s_{1,k})=\frac{\alpha^{S'_1+s_{1,k}}}{(\alpha)_{S'_1+S_2+s_{1,k}}}S^{s_{1,k}+s_{2,k}}_{r_k,0}S^{t_{1,k}+t_{2,k}+t_{3,k}}_{s_{1,k},0}.$$ Gibbs by sampling $s_{1,k}$ proportional to this for all

$$1\leq s_{1,k}\leq t_{1,k}+t_{2,k}+t_{3,k}.$$

• ignore the terms definition as it is not critical for understanding the question. The key is that we can compute all these quantities and $s_{1,k}$ is a count, so we can easily use all values in that inequality range.

My question basically is about how this sampling is done? This is my guess: do we take all possible values in the range specified by the inequality and plug them into the equation and compute it to get a corresponding probability (it's not a true probability but the normalizer seems to be the same for all quantities, hence we can turn the computed F to probability once we had all the computed values) for each sampled value? Then use the cumulative probability to sample from these?

For example, for $$1\leq s_{1,k}\leq 3$$ do we consider $F(s_{1,k}=1),F(s_{1,k}=2),F(s_{1,k}=3)$ and then define a cumulative probability $(F(s_{1,k}=1), F(s_{1,k}=1)+F(s_{1,k}=2), 1)$ and then generate a random number and based on the interval it falls into, we decide which sampled value for $s_{1,k}$ to take?

Yes, that works. To see this look at the following R code. I set $F(s = 1) = .2$, $F(s = 2) = .3$ and $F(s = 3) = .5$, and then I use your propose methods to draw values of $s$. I do 1000 repetitions, and calculated the proportion of times I see each of 1, 2, 3. These proportions should be close to .2, .3 and .5 respectively.

set.seed(10)

## P(s_1 = 1) = .2
## P(s_2 = 2) = .3
## P(s_3 = 3) = .5

reps <- 1000
count <- numeric(length = reps)

for(i in 1:reps)
{
## draw a uniform random number between 0 and 1
random <- runif(1)
if (random <= .2)
{
count[i] <- 1
} else {
if(random <= .5)
{
count[i] <- 2
} else {
count[i] <- 3
}
}
}

sum(count == 1)/reps
# [1] 0.196
sum(count == 2)/reps
# [1] 0.303
sum(count == 3)/reps
# [1] 0.501
• Nice! - I wonder if my understanding is correct though. Not sure if there is any other way of interpreting the statement in that tutorial. – user3639557 Jul 20 '16 at 15:03
• Your understanding is correct I think. – Greenparker Jul 20 '16 at 15:05