1
$\begingroup$

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

But the result was sad:

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.

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – whuber Jan 30 '15 at 19:48
  • $\begingroup$ 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. $\endgroup$ – bicycle1885 Jan 30 '15 at 20:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.