AMS 241, Fall 2010, Homework 2

Vadim von Brzeski vvonbrze@ucsc.edu November 8, 2010

1 Gibbs Sampling Equations

The model we are studying is:

y_{i}|θ_{i}, ϕ ∼ k_{N}(y_{i}; θ_{i}, ϕ), i = 1, ..., n
θ_{i}|G ∼ G, i = 1, ..., n
G|α, μ, τ^{2} ∼ DP (α, G_{0} = N(μ, τ^{2}))
α ∼ Gamma(a_{α}, b_{α})
μ ∼ N(a_{μ}, b_{μ})
τ^{2} ∼ InvGamma(a_{τ^{2}}, b_{τ^{2}})
ϕ ∼ InvGamma(a_{ϕ}, b_{ϕ})

1.1 Posterior θ_{i}

p(θ_{i}|y_{i}, \{θ_{k}, k ≠ i\}, α, ϕ, μ, τ^{2}, data) = \frac{q_{0}h(θ_{i}|ϕ, μ τ^{2}, y_{i}) + ∑_{j=1}^{n^{∗−}} n_{j}^{−}q_{j}δ_{θ_{j}^{∗−}}(θ_{i})}{q_{0} + ∑_{j=1}^{n^{∗−}} n_{j}^{−}q_{j}}  (1)

where

q_{j} = N(y_{i}; θ_{j}^{∗−}, ϕ)
q_{0} = α ∫_{−∞}^{∞} k(y_{i}; θ_{i}, ϕ)g_{0}(θ_{i}|μ, τ^{2})dθ_{i}
= α ∫_{−∞}^{∞} N(y_{i}; θ_{i}, ϕ)N(θ_{i}|μ, τ^{2})dθ_{i}
= αN(y_{i}; μ, ϕ + τ^{2})

and

h(θ_{i}|ϕ, μ τ^{2}, y_{i}) = C · k_{N}(y_{i}|θ_{i}, ϕ)g_{0}(θ_{i}|μ, τ^{2})
= C · N(y_{i}; θ_{i}, ϕ)N(θ_{i}|μ, τ^{2})
= N(θ_{i}; m, v)

where

m = \frac{μ/τ^{2} + y_{i}/ϕ}{1/τ^{2} + 1/ϕ} v = \frac{1}{1/τ^{2} + 1/ϕ}

So θ_{i} will either be a new distinct θ value drawn from N(θ_{i}; m, v) with probability proportional to q_{0}, or equal to an existing θ_{j} value with probability proportional to n_{j}^{−}q_{j}.

1.2 Posterior θ_{j}^{∗}

p(θ_{j}^{∗}|w, n^{∗}, μ, τ^{2}, ϕ, data) = C · [ ∏_{i:w_{i}=j} k(y_{i}|θ_{j}^{∗}, ϕ)]g_{0}(θ_{j}^{∗}|μ, τ^{2})

= C · [ ∏_{i:w_{i}=j} N(y_{i}|θ_{j}^{∗}, ϕ)]N(θ_{j}^{∗}|μ, τ^{2})

= N(θ_{j}^{∗}; m^{∗}, v^{∗})  (2)

where

m^{∗} = \frac{μ/τ^{2} + n_{j}ȳ_{j}/ϕ}{1/τ^{2} + n_{j}/ϕ} v = \frac{1}{1/τ^{2} + n_{j}/ϕ}

So for each individual cluster (component) j, we indepedently draw a value for θ_{j} from the distribution defined in equation [2].

1.3 Posterior ϕ

p(ϕ|θ, data) = C · p(ϕ) ∏_{i=1}^{n} k(y_{i}|θ_{i}, ϕ)

= C · InvGamma(a_{ϕ}, b_{ϕ}) ∏_{i=1}^{n} N(y_{i}|θ_{i}, ϕ)

= C · ϕ^{−(a_{ϕ}+1)}exp(−b_{ϕ}/ϕ) ∏_{i=1}^{n} ϕ^{−1/2}exp[\frac{−1}{2ϕ}(y_{i} − θi)^{2}]

= C · ϕ^{−(a_{ϕ}+n/2+1)}exp(−b_{ϕ}/ϕ)exp[ ∑_{i=1}^{n} \frac{−1}{2ϕ}(y_{i} − θi)^{2}]

= C · ϕ^{−(a_{ϕ}+n/2+1)}exp[\frac{−1}{ϕ}(b_{ϕ} + \frac{nv}{2} )]

= InvGamma(a_{ϕ} + n/2, b_{ϕ} + nv/2)  (3)

where

nv = ∑_{i=1}^{n}(y_{i} − θ_{i})^{2}

1.4 Posterior μ, τ^{2}

p(μ, τ^{2}|θ^{∗}) = C · p(μ)p(τ^{2}) ∏_{i=1}^{n^{∗}} g_{0}(θ_{j}^{∗}|μ, τ^{2})

= C · N(μ|a_{μ}, b_{μ})InvGamma(τ^{2}|a_{τ^{2}}, b_{τ^{2}}) ∏_{i=1}^{n^{∗}} N(θ_{j}^{∗}|μ, τ^{2})

We are not in a conjugate prior setting in this case for p(μ, τ^{2}|θ^{∗}), but we can easily compute the conditional distributions needed for Gibbs sampling:

p(μ|τ^{2}, θ^{∗}) = C · N(μ|a_{μ}, b_{μ}) ∏_{i=1}^{n^{∗}} N(θ_{j}^{∗}|μ, τ^{2})
= N(m_{μ}, s_{μ}^{2})  (4)

where

m_{μ} = \frac{a_{μ}/b_{μ} + n^{∗}θ̄^{∗}/τ^{2}}{1/b_{μ} + n^{∗}/τ^{2}} s_{μ}^{2} = \frac{1}{1/b_{μ} + n^{∗}/τ^{2}}

where n^{∗} is the number of distinct θ^{∗} values and θ̄^{∗} is the mean of the distinct θ^{∗} values.

p(τ^{2}|μ, θ^{∗}) = C · InvGamma(τ^{2}|a_{τ^{2}}, b_{τ^{2}}) ∏_{i=1}^{n^{∗}} N(θ_{j}^{∗}|μ, τ^{2})

Using the same algebraic manipulation as was done for the posterior of ϕ above, we have:

p(τ^{2}|μ, θ^{∗}) = InvGamma(a_{τ^{2}} + n^{∗}/2, b_{τ^{2}} + n^{∗}ν/2)  (5)

where

n^{∗}ν = ∑_{i=1}^{n^{∗}}(θ_{i}^{∗} − μ)^{2}

1.5 Posterior α

Here we just re-iterate the sampling scheme for α using the auxiliary variable scheme outlined by Escobar and West, 1995. First we draw η from :

p(η|α, data) = Beta(α + 1, n)  (6)

Define probability p = (a_{α} + n^{∗} − 1)/(n(b_{α} − log(η)) + a_{α} + n^{∗} − 1), and draw:

with probability p :  α|η, n^{∗}, data ∼ Gamma(a_{α} + n^{∗}, b_{α} − log(η))
with probability 1 − p :  α|η, n^{∗}, data ∼ Gamma(a_{α} + n^{∗} − 1, b_{α} − log(η))  (7)

1.6 Overview of Gibbs Sampling Scheme

A brief overview of the Gibbs sampling scheme is shown below. In our sampler, we maintain a state structure which consists of the θ_{i} values, the cluster (component) membership index for each y_{i}, and the values of all the parameters α, μ, τ^{2}, ϕ.

1.6.1 Initialization

At initialization, we assign each of the N y_{i} observations to belong to a different cluster, thus starting out with N separate clusters (components). Each cluster is initialized with a θ_{j}^{∗} drawn from g_{0} = N(μ, τ^{2}) distribution, given the starting values of μ, τ^{2}.

1.6.2 MCMC Simulation

We typically run for 1000 burn-in iterations, followed by 4000 monitoring iterations. We arrived at the 1000/4000 values after preliminary experiments that showed that increasing iterations (e.g. 2000/10000) did not appreciably change the results, indicating that the chain had reached convergence.

The sampling scheme consists of :

• For i = 1..n, sample θ_{i} from the mixed distribution defined by (1). Adjust cluster (component) membership indeces of the data if an existing component was dropped or if a new component was added.

• For j = 1..n^{∗}, sample θ_{j}^{∗} from (2).

• Sample ϕ from (3).

• Sample μ|τ^{2} from (4) and τ^{2}|μ from (5).

• Sample α using (6) and (7).

1.7 Predictive Distribution Calculations

The posterior predictive distribution for a new observation θ_{0} is given by the Polya-Urn scheme:

p(θ_{0}|n^{∗}, w, θ^{∗}, α, ϕ) = \frac{α}{α + n}G_{0}(θ_{0}|ϕ) + \frac{1}{α + n} ∑_{j=1}^{n^{∗}} n_{j}δ_{θ_{j}^{∗}}(θ_{0})  (8)

The posterior predictive distribution for a new y_{0} is given by :

p(y_{0}|data) = ∫ ∫ k(y_{0}; θ_{0}, ϕ)p(θ_{0}|n∗, w, θ^{∗}, α, ϕ)p(n^{∗}, w, θ^{∗}, α, ϕ|data)  (9)

Thus given B sets of samples from our MCMC output, we can obtain samples from p(y_{0}|data) as follows: for each set b of posterior parameter values in the MCMC output, first draw θ_{0,b} from p(θ_{0}|n∗, w, θ^{∗}, α, ϕ), and then draw y_{0,b} from p(y_{0}|θ_{0,b}, ϕ_{b}) = N(·; θ_{0,b}, ϕ_{b})

Samples from the prior predictive distribution p(y_{0}) are obtained as follows:

• Draw (μ, τ^{2}) from p(μ)p(τ^{2}). Since p(μ) and p(τ^{2}) are independent of each other, p(μ, τ^{2}) does not form a conjugate prior for G_{0} = N(μ, τ^{2}). So we independently sample b times:

• τ_{b}^{2} ∼ InvGamma(a_{τ^{2}}, b_{τ^{2}})

• μ_{b} ∼ N(a_{μ}, b_{μ})

• θ_{0,b} ∼ G_{0} = N(μ, τ^{2})

• y_{0,b} ∼ ∫ N(y_{0,b}; θ_{0,b}|ϕ)p(ϕ)dϕ = t_{ν}(θ_{0,b}, s^{2}), a t distribution with ν = 2a_{ϕ} degrees of freedom, mean θ_{0,b} and scale s^{2} = b_{ϕ}/a_{ϕ} The derivation of y_{0,b} ∼ t_{ν}(θ_{0,b}, s^{2}) is shown below:

∫ N(y; θ|ϕ)IG(ϕ; a, b)dϕ = ∫ \frac{1}{\sqrt{2πϕ}}exp[\frac{−1}{2ϕ}(y − θ)^{2}]\frac{b^{a}}{Γ(a)}ϕ^{−(a+1)}exp[\frac{−b}{ϕ}]dϕ

= C · ∫ ϕ^{−(a+1+1/2)}exp[\frac{−1}{2ϕ}((y − θ)^{2} + 2b)]dϕ

Letting z = \frac{A}{2ϕ}
ϕ = \frac{A}{2z}
|dϕ| = \frac{A}{2z^{2}} dz
where A = ((y − θ)^{2} + 2b) we get

= C · ∫ (\frac{A}{2z} )^{−(a+1+1/2)}exp[−z]\frac{A}{2z^{2}} dz

= C · (\frac{A}{2})^{−(a+1/2)} ∫ z^{(a+1/2−1)}exp(−z)dz

= C · (\frac{A}{2})^{−(a+1/2)}Γ(a + 1/2)

= C · Γ(a + 1/2)[\frac{1}{2}(y − θ)^{2} + b]^{−(a+1/2)}

= C · Γ(a + 1/2)b^{−(a+1/2)}[1 + \frac{(y − θ)^{2}}{2b} ]^{−(a+1/2)}

which we recognize as the kernel of a t-distribution. Letting a = a_{ϕ} = ν/2 and b = b_{ϕ} = νs^{2}/2 and filling in the constants, we get the familiar form of the t-distribution, namely:

y_{0,b} ∼ \frac{Γ(\frac{ν+1}{2})}{Γ(\frac{ν}{2})} \frac{1}{\sqrt{πν}s}[1 + \frac{(y − θ)^{2}}{νs^{2}} ]^{− \frac{ν+1}{2}}