# How can I sample multivariate binary variables such that sum of them follows a gamma distribution?

Edit:

Since the original question was confusing as whuber pointed out, let me rephrase the question with a Poisson distribution instead of a gamma distribution.

The energy term of a Poisson distribution ($\mathrm{Pois}(\lambda)$ defined as http://en.wikipedia.org/wiki/Poisson_distribution) will be written as:

$$E(k) = -k \log{\lambda} + \sum_{k^\prime=1}^{k} \log{k^\prime}$$ Here I assume $k = \sum_d x_d$.

Following the discussion of the original question, I wrote the following code to sample multivariate binary variables $\mathbf x$ such that $\sum_d x_d \sim \mathrm{Pois}(\lambda)$ and $\lambda = 25$.

function energy(x::Vector{Float64}, lambda::Float64)
k = sum(x)
en = -k * log(lambda)
for k′ in 1:k
en += log(k′)
end
en
end

srand(1234)
# dimension
D = 100
# shape and rate parameters
lambda = 25.
# generate 10 samples
for n in 1:10
# initialize variables
x = convert(Vector{Float64}, rand(D) .> 0.5)
# sample 1,000 times
for i in 1:1000
for d in 1:D
x[d] = 1
e1 = energy(x, lambda)
x[d] = 0
e0 = energy(x, lambda)
p1 = 1 / (1 + exp(e1 - e0))
if p1 > rand()
x[d] = 1
else
x[d] = 0
end
end
i += 1
end
@show sum(x)
end


The output was:

sum(x) => 33.0
sum(x) => 41.0
sum(x) => 43.0
sum(x) => 39.0
sum(x) => 38.0
sum(x) => 51.0
sum(x) => 45.0
sum(x) => 35.0
sum(x) => 40.0
sum(x) => 39.0


This is far from the expected mean ($\mathrm{E}[k] = \lambda = 25$).

I want to know why this happened and how to fix the problem. I really appreciate your advice.

Original Question:

I want to sample multivariate binary variables $\mathbf x$ ($x_d \in \{0, 1\}$) such that sum of them follows a given gamma distribution (i.e. $\sum_d x_d \sim \mathrm{Gamma}(\alpha, \beta)$ where $\alpha$ is a shape parameter and $\beta$ is a rate parameter defined as http://en.wikipedia.org/wiki/Gamma_distribution).

I thought a Gibbs sampler would be suitable for this purpose, and derived the following formula.

\begin{align} p(\mathbf x) &= \frac{1}{Z} \exp\{(\alpha - 1) \log \sum_d x_d - \beta \sum_d x_d\}\\ &= \frac{1}{Z} \exp{\{-E(\mathbf x)\}} \end{align} where $Z$ is a normalizing constant and $E(\mathbf x) = -(\alpha - 1) \log \sum_d x_d + \beta \sum_d x_d$.

The ratio of probabilities of $x_d$ given $\mathbf x_{-d}$ is \begin{align} \frac{p(x_d = 1 | \mathbf x_{-d})}{p(x_d = 0 | \mathbf x_{-d})} &= \frac{\exp{\{-E(x_d=1, \mathbf x_{-d})\}}}{\exp{\{-E(x_d=0, \mathbf x_{-d})\}}} \\ &= \exp{[-\{E(x_d=1, \mathbf x_{-d}) - E(x_d=0, \mathbf x_{-d})\}]} \\ &= \exp{(-\Delta E)} \end{align} where $\Delta E = E(x_d=1, \mathbf x_{-d}) - E(x_d=0, \mathbf x_{-d})$.

Since $p(x_d = 1|\mathbf x_{-d}) + p(x_d = 0|\mathbf x_{-d}) = 1$, we can see that

$$\frac{p(x_d = 1 | \mathbf x_{-d})}{1 - p(x_d = 1 | \mathbf x_{-d})} = \exp{(-\Delta E)}.$$

Hence $$p(x_d = 1|\mathbf x_{-d}) = \frac{1}{1 + \exp{(\Delta E)}}.$$

Then I wrote the following code to check the behavior with $\alpha = 1.0$ and $\beta = 0.1$ (written in the Julia language):

function energy(x::Vector{Float64}, alpha::Float64, beta::Float64)
-(alpha - 1) * log(sum(x)) + beta * sum(x)
end

srand(1234)
# dimension
D = 50
# shape and rate parameters
alpha = 1.0
beta = 0.1
# generate 10 samples
for n in 1:10
# initialize variables
x = convert(Vector{Float64}, rand(D) .> 0.5)
# sample 1,000 times
for i in 1:1000
for d in 1:D
x[d] = 1
e1 = energy(x, alpha, beta)
x[d] = 0
e0 = energy(x, alpha, beta)
p1 = 1 / (1 + exp(e1 - e0))
if p1 > rand()
x[d] = 1
else
x[d] = 0
end
end
end
@show sum(x)
end


sum(x) => 15.0
sum(x) => 23.0
sum(x) => 25.0
sum(x) => 26.0
sum(x) => 26.0
sum(x) => 25.0
sum(x) => 26.0
sum(x) => 23.0
sum(x) => 27.0
sum(x) => 27.0


Apparently this does not follow the expected mean value of $\mathrm E[X] = \frac{\alpha}{\beta} = 10$ and when D was increased, sum(x) was also increased.

I guess I misunderstand something fundamental, but I couldn't figure out what it is for several hours. I will appreciate any advice.

Thank you for reading this long question.

• It is impossible for any function of a finite number $0,1$ variables (which would be discrete) to follow a Gamma distribution (which is continuous). This makes your question impossible to understand. – whuber Jan 30 '15 at 19:48
• Ah, that may may be a point I missed. I thought that I can use a continuous distribution can be used to model a discrete data and the fitted model "looks like" same shape of the continuous one. – bicycle1885 Jan 30 '15 at 20:28