discrete Gibbs; Markov blanket with DB-like data; Dawid & Skene (1979) model
124 views
Skip to first unread message
Bob Carpenter
Oct 26, 2012, 3:17:53 PM10/26/12
to stan...@googlegroups.com
Time to do discrete parameters (and another Bob mail-novel).
First an example that's more indicative of coding
Markov blankets than LDA, which is relatively easy.
Dawid and Skene (1979) introduced a model of
data coding by multiple coders where the truth
wasn't known. Here's what it would look like in
Stan if Stan did discrete sampling:
---------------------------------
data {
int<lower=2> K; // number of categories
int<lower=1> I; // number of items
int<lower=1> J; // number of coders
int<lower=1> N; // number of observations
int<lower=1,upper=J> jj[N]; // coder for observation n
int<lower=1,upper=I> ii[N]; // item for observation n
int<lower=1,upper=K> y[N]; // label for observation n
}
parameters {
simplex[K] pi; // prevalence of categories
simplex[K] theta[J,K]; // response of anotator j to category k
int<lower=1,upper=K> z[I]; // true category of item i
}
model {
for (i in 1:I)
z[i] ~ categorical(pi);
for (n in 1:N)
y[n] ~ categorical(theta[jj[n],z[ii[n]]]);
}
-------------------------------------
It uses the standard DB-like encoding of the data with y[N],
ii[N] and jj[N]. In general, every item is only coded by
a subset of the coders.
The relevant lines of the model for sampling a given discrete
z[j] for j in 1:I are
* i s.t. (i == j) in the first loop, and
* n s.t. (ii[n] == j) in the second loop.
It's the second of these that's particularly problematic
because it's data dependent --- there's no way to annotate the
relevant loop indexes at compile time.
What I did to run the model(* see P.S.) in RStan 1.0.2 (I'm a very happy
RStan user) was sum out the z[i]:
---------------------------------
data {
int<lower=2> K; // number of categories
int<lower=1> I; // number of items
int<lower=1> J; // number of coders
int<lower=1> N; // number of observations
int<lower=1,upper=J> jj[N]; // coder for observation n
int<lower=1,upper=I> ii[N]; // item for observation n
int<lower=1,upper=K> y[N]; // label for observation n
}
parameters {
simplex[K] pi; // prevalence of categories
simplex[K] theta[J,K]; // response of anotator j to category k
}
model {
real cat_log[I,K]; // vector of log probs (up to const) for item i
for (k in 1:K)
for (i in 1:I)
cat_log[i,k] <- log(pi[k]);
for (n in 1:N)
for (k in 1:K)
cat_log[ii[n],k] <- cat_log[ii[n],k] + log(theta[jj[n],k,y[n]]);
for (i in 1:I)
lp__ <- lp__ + log_sum_exp(cat_log[i]);
}
------------------------------------
The log_sum_exp just helps with arithmetic stability.
This is the same kind of code as we'd need for Gibbs.
(This isn't actually how I coded the first loop -- it can
be optimized by computing log(pi[k]) once for i=1 and
reusing it, but the above is simpler to understand.)
We could instead of this summing out as in the current Stan,
declare the z[I] parameter and have a Gibbs block like this:
------------------------------
gibbs {
real cat_log[I,K]; // vector of log probs (up to const) for item i
for (k in 1:K)
for (i in 1:I)
cat_log[i,k] <- log(pi[k]);
for (n in 1:N)
for (k in 1:K)
cat_log[ii[n],k] <- cat_log[ii[n],k] + log(theta[jj[n],k,y[n]]);
for (i in 1:I)
z[i] ~ categorical_log(cat_log[i]);
}
--------------------------------
where the "_log" in "categorical_log" here just indicates the params
are on the (unnormalized) log scale rather than a simplex:
categorical_log(z|alpha) = categorical(z|softmax(exp(alpha)))
Implementing discrete sampling this way would be:
* very very easy to implement, and
* optimally efficient to sample (given we're doing Gibbs/slice).
Basically, we'd just do it the same way we'd do generated
quantities -- the ~ just indicates a randomly generated variate.
And with the Markov blanket in hand, we can't do much better
computation-wise. Users can even control shared subcomputations.
The big downside is that it's an enormous burden on the users compared
to the "natural" specification of the model in the first example.
One upside is that by the time they've figured out Gibbs,
they can probably marginalize out the discrete params!
You could add the same calcs to generated quantities and either
save the category distributions for z[i] or actually sample.
A serious problem with the suggested model (not the approach,
just the particular coding above) is that it computes all the probs
first, then samples them all. Usually Gibbs is implemented in a
rotated style, where you sample, e.g.,
x[1], y[1], z[1] init
x[2] | y[1], z[1]
y[2] | x[2], z[1]
z[2] | x[2], y[2]
x[3] | y[2], z[2]
...
What the above model codes is this:
x[2] | y[1], z[1]
y[2] | x[1], z[1]
z[2] | x[1], y[1]
Nice for parallelism, but not so nice for mixing.
We can solve that problem with the same kind of gibbs block.
To make the models efficient in these non-full data designs, we'd
need to compute DB-like indexes. All we need is to be able to find,
given j in 1:I, the n such that (ii[n] == j). It would even be
possible to compute these indexes in the transformed data block.
More on this later.
BUGS and JAGS compute the Markov blanket for each variable
by expanding out the plates of the graphical model
explicitly. If we restrict the expressiveness of Stan a bit,
we could do the same kind of thing.
- Bob
P.S. This model rocks in Stan, by the way -- it's what I
showed you the fits for. We've got about 30 datasets
for different words, with Mechanical Turk coders identifying
the dictionary meaning for a word in context. We have I=1000, K=7,
N=25,000 and it converges in fewer than 10 steps and mixes
very very well (n_eff is between 1/3 and 2/3 the number of draws!).
It takes about 2 minutes for 500 iterations of a chain, and
even runs on Becky's ancient 2GB 32-bit Windows machine. The discrete
form of the model is a bit of a strain on JAGS -- Becky couldn't
get it in memory and it took hours to run without good
mixing on my machine -- I haven't done the fair comparison
of using the marginalized out model in JAGS. The obvious
next step is to do something hierarchical over the 30
data sets -- I already fit hierarchical models for theta within
datasets. We also want to fit item effects, so I'm going to
move it all over to a more multilevel-like logistic structure
as Andrew and Jennifer first suggested about four years ago!
(I've done this in JAGS already, but didn't have enough data to
pin down the difficulty per item -- that's why we collected so
much data, 25 labels per word instance.)
First an example that's more indicative of coding
Markov blankets than LDA, which is relatively easy.
Dawid and Skene (1979) introduced a model of
data coding by multiple coders where the truth
wasn't known. Here's what it would look like in
Stan if Stan did discrete sampling:
---------------------------------
data {
int<lower=2> K; // number of categories
int<lower=1> I; // number of items
int<lower=1> J; // number of coders
int<lower=1> N; // number of observations
int<lower=1,upper=J> jj[N]; // coder for observation n
int<lower=1,upper=I> ii[N]; // item for observation n
int<lower=1,upper=K> y[N]; // label for observation n
}
parameters {
simplex[K] pi; // prevalence of categories
simplex[K] theta[J,K]; // response of anotator j to category k
int<lower=1,upper=K> z[I]; // true category of item i
}
model {
for (i in 1:I)
z[i] ~ categorical(pi);
for (n in 1:N)
y[n] ~ categorical(theta[jj[n],z[ii[n]]]);
}
-------------------------------------
It uses the standard DB-like encoding of the data with y[N],
ii[N] and jj[N]. In general, every item is only coded by
a subset of the coders.
The relevant lines of the model for sampling a given discrete
z[j] for j in 1:I are
* i s.t. (i == j) in the first loop, and
* n s.t. (ii[n] == j) in the second loop.
It's the second of these that's particularly problematic
because it's data dependent --- there's no way to annotate the
relevant loop indexes at compile time.
What I did to run the model(* see P.S.) in RStan 1.0.2 (I'm a very happy
RStan user) was sum out the z[i]:
---------------------------------
data {
int<lower=2> K; // number of categories
int<lower=1> I; // number of items
int<lower=1> J; // number of coders
int<lower=1> N; // number of observations
int<lower=1,upper=J> jj[N]; // coder for observation n
int<lower=1,upper=I> ii[N]; // item for observation n
int<lower=1,upper=K> y[N]; // label for observation n
}
parameters {
simplex[K] pi; // prevalence of categories
simplex[K] theta[J,K]; // response of anotator j to category k
}
model {
real cat_log[I,K]; // vector of log probs (up to const) for item i
for (k in 1:K)
for (i in 1:I)
cat_log[i,k] <- log(pi[k]);
for (n in 1:N)
for (k in 1:K)
cat_log[ii[n],k] <- cat_log[ii[n],k] + log(theta[jj[n],k,y[n]]);
for (i in 1:I)
lp__ <- lp__ + log_sum_exp(cat_log[i]);
}
------------------------------------
The log_sum_exp just helps with arithmetic stability.
This is the same kind of code as we'd need for Gibbs.
(This isn't actually how I coded the first loop -- it can
be optimized by computing log(pi[k]) once for i=1 and
reusing it, but the above is simpler to understand.)
We could instead of this summing out as in the current Stan,
declare the z[I] parameter and have a Gibbs block like this:
------------------------------
gibbs {
real cat_log[I,K]; // vector of log probs (up to const) for item i
for (k in 1:K)
for (i in 1:I)
cat_log[i,k] <- log(pi[k]);
for (n in 1:N)
for (k in 1:K)
cat_log[ii[n],k] <- cat_log[ii[n],k] + log(theta[jj[n],k,y[n]]);
for (i in 1:I)
z[i] ~ categorical_log(cat_log[i]);
}
--------------------------------
where the "_log" in "categorical_log" here just indicates the params
are on the (unnormalized) log scale rather than a simplex:
categorical_log(z|alpha) = categorical(z|softmax(exp(alpha)))
Implementing discrete sampling this way would be:
* very very easy to implement, and
* optimally efficient to sample (given we're doing Gibbs/slice).
Basically, we'd just do it the same way we'd do generated
quantities -- the ~ just indicates a randomly generated variate.
And with the Markov blanket in hand, we can't do much better
computation-wise. Users can even control shared subcomputations.
The big downside is that it's an enormous burden on the users compared
to the "natural" specification of the model in the first example.
One upside is that by the time they've figured out Gibbs,
they can probably marginalize out the discrete params!
You could add the same calcs to generated quantities and either
save the category distributions for z[i] or actually sample.
A serious problem with the suggested model (not the approach,
just the particular coding above) is that it computes all the probs
first, then samples them all. Usually Gibbs is implemented in a
rotated style, where you sample, e.g.,
x[1], y[1], z[1] init
x[2] | y[1], z[1]
y[2] | x[2], z[1]
z[2] | x[2], y[2]
x[3] | y[2], z[2]
...
What the above model codes is this:
x[2] | y[1], z[1]
y[2] | x[1], z[1]
z[2] | x[1], y[1]
Nice for parallelism, but not so nice for mixing.
We can solve that problem with the same kind of gibbs block.
To make the models efficient in these non-full data designs, we'd
need to compute DB-like indexes. All we need is to be able to find,
given j in 1:I, the n such that (ii[n] == j). It would even be
possible to compute these indexes in the transformed data block.
BUGS and JAGS compute the Markov blanket for each variable
by expanding out the plates of the graphical model
explicitly. If we restrict the expressiveness of Stan a bit,
we could do the same kind of thing.
- Bob
P.S. This model rocks in Stan, by the way -- it's what I
showed you the fits for. We've got about 30 datasets
for different words, with Mechanical Turk coders identifying
the dictionary meaning for a word in context. We have I=1000, K=7,
N=25,000 and it converges in fewer than 10 steps and mixes
very very well (n_eff is between 1/3 and 2/3 the number of draws!).
It takes about 2 minutes for 500 iterations of a chain, and
even runs on Becky's ancient 2GB 32-bit Windows machine. The discrete
form of the model is a bit of a strain on JAGS -- Becky couldn't
get it in memory and it took hours to run without good
mixing on my machine -- I haven't done the fair comparison
of using the marginalized out model in JAGS. The obvious
next step is to do something hierarchical over the 30
data sets -- I already fit hierarchical models for theta within
datasets. We also want to fit item effects, so I'm going to
move it all over to a more multilevel-like logistic structure
as Andrew and Jennifer first suggested about four years ago!
(I've done this in JAGS already, but didn't have enough data to
pin down the difficulty per item -- that's why we collected so
much data, 25 labels per word instance.)
Matt Hoffman
Oct 26, 2012, 8:44:16 PM10/26/12
to stan...@googlegroups.com
I like the Gibbs block idea. But I'm not sure that this example is
quite correct. This is a valid Gibbs sampler:
> x[1], y[1], z[1] init
>
> x[2] | y[1], z[1]
> y[2] | x[2], z[1]
> z[2] | x[2], y[2]
>
> x[3] | y[2], z[2]
This is not:
> x[2] | y[1], z[1]
> y[2] | x[1], z[1]
> z[2] | x[1], y[1]
It's not just a question of mixing speed, it's a question of
correctness. Consider two binary variables x and y, with contingency
table
p(x=0, y=0) = 0.25
p(x=1, y=0) = 0.25
p(x=0, y=1) = 0.5
p(x=1, y=1) = 0
suppose x[1] = 0, y[1] = 0. Then if we resample x and y given their
previous values simultaneously, then there's a positive probability
(1/3) of arriving in the 0-probability state x=1, y=1. Clearly that's
not ok.
I'm sure that allowing users to implement their own Gibbs samplers
will lead to all kinds of nasty bugs getting posted to stan-users. But
that certainly doesn't mean we shouldn't do it!
Matt
quite correct. This is a valid Gibbs sampler:
> x[1], y[1], z[1] init
>
> x[2] | y[1], z[1]
> y[2] | x[2], z[1]
> z[2] | x[2], y[2]
>
> x[3] | y[2], z[2]
> x[2] | y[1], z[1]
> y[2] | x[1], z[1]
> z[2] | x[1], y[1]
correctness. Consider two binary variables x and y, with contingency
table
p(x=0, y=0) = 0.25
p(x=1, y=0) = 0.25
p(x=0, y=1) = 0.5
p(x=1, y=1) = 0
suppose x[1] = 0, y[1] = 0. Then if we resample x and y given their
previous values simultaneously, then there's a positive probability
(1/3) of arriving in the 0-probability state x=1, y=1. Clearly that's
not ok.
I'm sure that allowing users to implement their own Gibbs samplers
will lead to all kinds of nasty bugs getting posted to stan-users. But
that certainly doesn't mean we shouldn't do it!
Matt
Bob Carpenter
Oct 27, 2012, 5:25:33 PM10/27/12
to stan...@googlegroups.com
On 10/26/12 8:44 PM, Matt Hoffman wrote:
> I like the Gibbs block idea. But I'm not sure that this example is
> quite correct. This is a valid Gibbs sampler:
>
>> x[1], y[1], z[1] init
>>
>> x[2] | y[1], z[1]
>> y[2] | x[2], z[1]
>> z[2] | x[2], y[2]
>>
>> x[3] | y[2], z[2]
>
> This is not:
>
>> x[2] | y[1], z[1]
>> y[2] | x[1], z[1]
>> z[2] | x[1], y[1]
the joint
x[2], y[2], z[2] |
to do them all at once, or you could block them as
x[2], y[2] | z[1]
z[2] | x[2], y[2]
...
> I'm sure that allowing users to implement their own Gibbs samplers
> will lead to all kinds of nasty bugs getting posted to stan-users.
> But
> that certainly doesn't mean we shouldn't do it!
optimal Gibbs, although it leans very heavily on the
user's coding.
- Bob
Andrew Gelman
Oct 27, 2012, 9:42:09 PM10/27/12
to stan...@googlegroups.com
Hi--I have a thought. Suppose we go with the current plan and program discrete data by requiring the user to specify the relevant portions of the log posterior. Then there is the potential for inconsistency, if the user mistypes or misunderstands the model. Here's a way we could check this. When the density ratio is calculated (for Metrop accept/reject or for slice or whatever), we could usually use the shorter fast formula suppled by the user. But every once in awhile we could compute the ratio using the full joint posterior. If these two ratios are ever different (by more than some small round-off value), this would indicate a mistake, and we could print out an error message.
A
> --
>
A
>
Matt Hoffman
Oct 27, 2012, 11:04:41 PM10/27/12
to stan...@googlegroups.com
I really like this idea. We could also have it as a command-line
option (or the R equivalent), similar to the current functionality for
testing gradients. Then, once you're confident your Gibbs sampler's
implemented correctly you can take the water wings
(http://en.wikipedia.org/wiki/Inflatable_armbands) off.
This could be useful even for 100% discrete models like collapsed LDA!
I'm always terrified (often with justification) that I'm implementing
my samplers incorrectly, but I'm also always too lazy to implement a
Cook et al. style validation scheme. It'd be great to know that Stan's
watching my back.
Matt
> --
>
>
option (or the R equivalent), similar to the current functionality for
testing gradients. Then, once you're confident your Gibbs sampler's
implemented correctly you can take the water wings
(http://en.wikipedia.org/wiki/Inflatable_armbands) off.
This could be useful even for 100% discrete models like collapsed LDA!
I'm always terrified (often with justification) that I'm implementing
my samplers incorrectly, but I'm also always too lazy to implement a
Cook et al. style validation scheme. It'd be great to know that Stan's
watching my back.
Matt
>
>
Bob Carpenter
Oct 28, 2012, 1:41:37 AM10/28/12
to stan...@googlegroups.com
Yes, this is a very good idea for all the reasons
Matt and Andrew say.
Somebody (Michael? Andrew?) brought up the idea of
using the full log prob to check the conditionals during
the meeting --- it's in my notebook with stars next
to it.
We have similar tests in many places using auto-dif on
simple definitions to test hand-derived derivatives.
It's also very much like the finite-diff gradient
check Matt originally suggested would be so helpful.
If we go to NASA-style triple redundancy, we could
probably write a paper about testing our derivatives:
by hand, by auto-dif, by finite diffs. Qudruple
redundancy could be had with symbolic derivatives
from sympy or Mathematica.
Speaking of tests, we really need to do the first
thing Michael said we should do for debugging and
allow dumping of each leapfrog step with unconstrained
parameters and gradients (we could also dump constrained
params).
- Bob
Matt and Andrew say.
Somebody (Michael? Andrew?) brought up the idea of
using the full log prob to check the conditionals during
the meeting --- it's in my notebook with stars next
to it.
We have similar tests in many places using auto-dif on
simple definitions to test hand-derived derivatives.
It's also very much like the finite-diff gradient
check Matt originally suggested would be so helpful.
If we go to NASA-style triple redundancy, we could
probably write a paper about testing our derivatives:
by hand, by auto-dif, by finite diffs. Qudruple
redundancy could be had with symbolic derivatives
from sympy or Mathematica.
Speaking of tests, we really need to do the first
thing Michael said we should do for debugging and
allow dumping of each leapfrog step with unconstrained
parameters and gradients (we could also dump constrained
params).
- Bob
Michael Betancourt
Oct 29, 2012, 1:33:04 AM10/29/12
to stan...@googlegroups.com
Yeah, we had discussed adding some means to validate the consistency of the conditional distributions, if we allow users to define them independently of the joint.
Personally I think this is a great idea -- for every freedom we give the user we introduce some means of validation to avoid the inevitable errors. Not only for conditional distributions but also, for example, implementing gradients by hand. If the issue is speed we'd compare the user's implementation against autodiff, and if they were implementing by hand because autodiff couldn't handle something we'd compare against finite differences.
> --
>
>
Personally I think this is a great idea -- for every freedom we give the user we introduce some means of validation to avoid the inevitable errors. Not only for conditional distributions but also, for example, implementing gradients by hand. If the issue is speed we'd compare the user's implementation against autodiff, and if they were implementing by hand because autodiff couldn't handle something we'd compare against finite differences.
>
>
Reply all
Reply to author
Forward
0 new messages