Discrete parameters...again!

[Ole Petter Hansen]

I'm working on models with discrete parameters that can't be marginalized out. For this I'm using BUGS, which seems to be best way to go. And I've seen the posts about the difficulties in adjusting HMC to dealing with discrete parameters. 

Stan has several nice features that I haven't found anywhere else - even without taking the HMC-bit into account (checks on inputs, good debugging features, documentation, big and active community, vectorisation...). So I was wondering if it could be possible to add Metropolis-Hastings and/or Gibbs sampling to Stan? Not as a mix of MH and HMC, but simply plain, old MH. 

So if a parameter is given a discrete distribution, Stan will just become an beefed up, simple sampler. Could well print a big, fat warning of that HMC is turned off. 

If sometime in the future you can adjust HMC to handle discrete stuff that is great too. But for now I think simply moving discrete models into Stan in a simple way will a major leap forward. 

Ole-Petter

[Bob Carpenter]

Yes, it would be possible to add a random-walk Metropolis for
discrete parameters with a lot of work:

* major change in language to allow discrete parameters
* matching changes in code generator plus model interface
* sampler with adaptation
* interface changes to disallow models with discrete parameters
in HMC and to call Metropolis
* doc (reference manual, algorithms, user's guide)
* unit tests
* integration tests

We even have a Metropolis implementation for continuous
parameters on a branch. We could never figure out how to
do adaptation in a robust way so that it could fit the
BUGS or ARM example models. So we gave up.

BUGS/JAGS use Gibbs. That's much much harder, because we
don't have a graphical modeling language like BUGS/JAGS (contrary
to what a lot of people think). So there's no good way to
pull the Markov blanket out for a variable.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Ole Petter Hansen]

Ok, thanks anyway. Had a hope this wouldn't be that much work, given all of the infrastructure already in place in Stan. 

-op

[Bob Carpenter]

It's overall not that much work compared to the other
things we do, but the real bottleneck is that we couldn't
figure out how to do adaptation for random-walk Metropolis
that actually worked in a range of models.

- Bob

> On Feb 2, 2016, at 2:36 PM, Ole Petter Hansen <olepetter...@outlook.com> wrote:
>
> Ok, thanks anyway. Had a hope this wouldn't be that much work, given all of the infrastructure already in place in Stan.
>
> -op
>

[bee]

Is it possible to marginalize out discrete parameters automatically when the counted number of case is less than 1000 or something?

-Bee

[Andrew Gelman]

For discrete parameters we were talking about slice sampling. But I recall the big issue as not the algorithm, but the programming effort required introduce a new parameter type in Stan, keep track of it, etc.

[Ole Petter Hansen]

Ok,well now you know that discrete-parameter people (at least one) would like to move over to Stan - even if sampler is not better than BUGS. 

[Dustin Tran]

I also recall discussion at some point of continuous relaxation stuff with HMC. Right, even with VI it’s easy to do discrete parameters. If only we can resolve the parameter labelling problem!

Dustin

On Feb 3, 2016, at 2:17 AM, Ole Petter Hansen <olepetter...@outlook.com> wrote:

Ok,well now you know that discrete-parameter people (at least one) would like to move over to Stan - even if sampler is not better than BUGS. 

[Bob Carpenter]

The issue is that because Stan does not define a graphical
model (and in fact allows a far more general class of models
to be written), we can't easily pull out the relevant
probabilities.

We can manage all the data types and everything else, but
the real killer is an easy way to pull out the conditional
density of a single discrete parameter or all of them.
And then how to do adaptation. This goes for Gibbs with
slice sampling (adapt step-out size, need conditional of each
discrete parameter separately) or random-walk Metropolis (adapt proposal
sizes, need conditional of all discrete parameters together).
And either way, we need to be able to control lower-bound, upper-bound
and interval-bounded constraints.

And we have no idea whatsoever how it would interact with
the global adaptation we do for HMC. You (Andrew) created
some kind of technique, I think you called it burnin or something,
to deal with throwing away intermediate samples as you
changed discrete regimes.

[Andrew Gelman]

Hi, I remember Bob and I had a long discussion of this in Paris a couple years ago. We came up with some clever ideas that could speed the computation in various sorts of models. My impression was that the real thing that stopped us from doing it was programming effort. If the programming were easy enough (and if it didn’t require taking apart Stan), I think we could do a reasonable slice sampling implementation that would do the job for lots of problems (and run in general, it would just be slow for hard problems)
A

[Michael Betancourt]

I think we could do a reasonable slice sampling implementation that would do the job for lots of problems (and run in general, it would just be slow for hard problems)

Not only slow but also hard to diagnose, with little theoretical

guarantees of accurate estimates.  Ultimately we just don’t

yet know how to effectively sample from complex discrete 

distributions as we do in the continuous case!

[Andrew Gelman]

Yes, this came up too.  But . . . to the extent that we’re competing with Bugs etc that _are_ handing discrete parameters, it would be good if we could have some implementation of them also.

[Bob Carpenter]

Yes, we could alternate slicing discrete and continuous
parameters. It'll be hugely hugely inefficient because
we can't find the Markov blanket for a discrete parameter.
I don't know how to adapt slice sampling, but I guess Andrew
has a handle on that from what he says below.

The drawback is that we'll evaluate the entire density function
for every step of slice (not just the Markov blanket).

What you (Andrew) wanted was

1. some kind of "adiabatic" thing for the phase when the
discrete params change; I never understood why, so I may
have this wrong

What I wanted was

2. a good way to compute the Markov blanket.

One thing we discussed for (2) was having users code out the
conditionals directly, just like a stats paper from 10 years
ago. That'd solve the problem at the expense of a lot of work
on the language side and a huge path for mistakes and being
yelled at for having a dumb language.

Andrew thought you'd be able to sort out what you needed
from the code, but it's an undecidable problem in general.
So we'd need subcases. So the other thing we talked about was

3. just restricting discrete models to graphical models.

That would solve (2), but not (1). But now we're talking
about build a graphical modeling library.

- Bob

[Andrew Gelman]

We talked about various options, the simplest of which was to allow the user to tag the program or divide it into blocks to tell it which lines of code needed to be evaluated when updating each discrete parameter. I think that if we felt it was important to allow Stan to work with discrete parameters, we could come up with various reasonable solutions. But we were never convinced it was such a big deal, so we never decided to do it. I _don’t_ think the solution is to go to graphical modeling. I mean, yes, I do want a graphical modeling wrapper for Stan, but I don’t like the idea that Stan internally would use graphical models. To me, that’s too limiting.

[Bob Carpenter]

We seem to have merged two threads along the way.

I don't know what you mean by "internally". A graphical
model with continuous parameters would be compiled down to
a Stan model, which would be run. A graphical mdoel with
discrete parameters would be able to generate all the conditional
densities we need to run Gibbs on top of that. We could
allow other ways to set those conditionals, like just using
the full densities.

- Bob

[Andrew Gelman]

I’d prefer to have the graphcial model in a wrapper that takes the graphcial model, writes the log posterior density, then runs Stan. I’d like “Stan” to just take an objective function (a log-posterior) without any requirement that it be a graphical model.

[Bob Carpenter]

We can't do the efficient sampling with only the
total log density because there's no way to pull out the
Markov blanket for each discrete parameter to compute its
conditionals efficiently.

But if the "objective function" can be extended to include
conditionals for discrete parameters, we could do it efficiently.
(Yes, I know it's the wild West for discrete parameters; I'm just
talking about efficiency relative to other systems like BUGS).

- Bob

[Emmanuel Charpentier]

[ From the perspective of a strictly applied statistics practitioner ]

If I follow you, we should have a compiler for some graphical model description language (let's call it SuperBUGS, for example) which might :

1. extract relevant information not directly representable in Stan terms (e. g. Markov blanket of discrete parameters),
2. write a Stan program encoding the Stan-representable information,
3. run this program
4. read the resultant sample
5. use the non-Stan-representable information to massage this output in a form closer to what is expected from a graphical-model perspective (i. e. the SuperBUGS representation).

Alternatively, one could write a program working directly in the SuperBUGS description (possibly creating a graph of stochastic nodes, as do current BUGS implementations), but using NUTS and other HMC samplers as they happen to fit part of the problem.

In both cases, we would have to subsume BUGS and Stan, in two directions :

• competence (i. e. the class of models representable in SuperBUGS should have both Stan-representable models and BUGS-representable models as subsets),
• performance (i. e. better execution times/CPU/RAM foothold) that should be at least as good as BUGS and Stan for their respective classes of representable models.

The first item is the most important.

In both cases, what the end-user would see would be :
the ease/difficulty of representing his/her model in SuperBUGS (i. e. preparation to step 1 above)
the time/computing resources necessary to go from 1 to 5.

Am I right in my reading ?

[ A question, by the way : do you have examples of Stan-expressible models that *cannot* be represented as graphical models (not necessarily directed acyclic graphs as in BUGS), however awkwardly ? ]

My two cents (and maybe some sense... )

--
Emmanuel Charpentier

[Dustin Tran]

I don’t think we want to internally encode any program as a graphical model. That's more restricting than the current class that Stan supports (the typical example that the PP community likes is programs with recursive functions). If I interpreted Andrew correctly, all he suggested is a wrapper for people to write in some form of syntax that looks like a graphical model, and internally it only outputs the log posterior density.

Dustin

[Bob Carpenter]

Right --- no way to do a superset of Stan with a graphical
model that'll let us do the BUGS-like computations we want.

I see three possible benefits for having a graphical modeling
language.

1. Allow metaprogramming. That is, write programs that build
up graphical models. There's a natural structure of nodes
and links to manipulate. I really have no idea how useful this
would be, so I'd love to see some real examples from, say,
PyMC or NIMBLE, both of which should allow this.

The other use for metaprogramming would be program transformation.
We have some hope of doing that now because we get an
intermediate abstract syntax tree for all the constructs in a Stan
program.

2. Allow forward sampling. Right now, we can't figure out what's
a prior and what's higher level because there's no graphical structure.
We'd like to be able to automatically simulate data sets.

3. Allow easy Markov blanket computations so we can do
conditional sampling for discrete variables, which would allow us
to have missing data, variable selection (probably won't be able
to fit in any software), direct mixture model representations, etc.

The main reason we haven't done this is that the graphical
modeling language would be more restrictive. For example, you
could write the expectations out for an HMM or a CRF directly using
the forward/backward algorithm in Stan for efficient fitting
and code it directly in Stan's language.

- Bob

[Andrew Gelman]

Hi, Emmanuel.  I’d like to think of Stan as the inference engine which can work with any objective function (which is considered as a log posterior density).  What you’re calling Superbugs would be a way to write graphical models and translate them into Stan for fitting and also be able to do things like fake-data simulation and posterior predictive checking.

In answer to your final question:  I’m not sure.  But I do think that in expressing a model it can be convenient to just add terms to the log-posterior.  Expressing in terms of a graphical model can at times be awkward, compared to simply adding penalty terms to the log-posterior.

A

On Feb 10, 2016, at 2:50 PM, Emmanuel Charpentier <emanuel.c...@gmail.com> wrote: