a (fairly) painless way to implement discrete sampling?

[Andrew Gelman]

Hi all (especially Bob). During the times we've discussed discrete sampling, one possibility we've talked about is to work with latent continuous variables. We've talked about this as a possibility but I always had the sense that it might be tricky to get the distributions right.

But, thinking about it, I think it's close to trivial. For simplicity, suppose our model has only two parameters: a N-dimensional discrete parameter a (which can take on the values 1,2,…,K) and an M-dimensional continuous parameter b. We can define a new continuous parameter z defined on [0,1], along with K probabilities p (i.e., simplex[K] p).

Thus, we expand our model from:

parameters {
integer<lower=1,upper=K> a[N];
real b[M];
}

to

parameters {
real z[N];
simplex[K] p;
real b[M];
}

Now, z relates to a in a very simple way. In the model, define

real<lower=0,upper=1> cutpoints[K]
cutpoints <- cumulative_sum (p)

Then, y=1 if z<cutpoints[1],
y=2 if z is between cutpoints[1] and cutpoints[2]
. . .
y=K if z>cutpoints[K-1]

This can all be implemented in Stan.

That means, I think, that we could implement discrete variables internally: if we have an integer-valued parameter "a", we just define the corresponding "z" and "p" parameters and set up the model as above.

All this so far only works with finitely-valued discrete parameters, but I think we could do countably-infinite sets too, by just defining a preset Zeno-like sequence of cutpoints.

What do you think?

A

P.S. I tried this out by implementing the standard mixture model using the latent-data formulation and then testing it out on fake data. It seemed to work, but it was very very slow to converge. So that's another concern. But if we can get the model set up, improving convergence might be possible using various Matt and Mike tricks...

[Matt Hoffman]

I think this is valid from a language perspective, but it probably
won't ever be fast from an algorithmic perspective. The problem is
that the posterior will become discontinuous unless a is completely
isolated from the rest of the model. (In which case there's not much
point.) Each time you cross a cutpoint you get an abrupt change in
energy, which basically hoses your chance of maintaining a high
acceptance probability.

That said, I seem to remember Mike talking about an integration scheme
that might be able to handle that sort of problem.

Matt

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

[Andrew Gelman]

One question is: How slow would it be on problems of interest? If it is implementable and slow but not prohibitively slow, it could be a useful step along the way, so that users who really want to implement discrete variables will be able to do so. But of course if Mike could whip up a better algorithm, that would be great.
For some models (for example, or multinomial logit), I could actually imagine the continuous parameterization being not too bad, in that the underlying continuous parameter has a statistical interpretation.
A

[Matt Hoffman]

My guess is that performance will be closer to Metropolis, since the
discontinuities are cutting off any gradient signal from the parts of
the model that depend on the discrete variables. Might still work, but
might not be better than JAGS.

[Andrew Gelman]

Performance comparable to Metropolis/Jags isn't bad, in that at least it will run!
A

[Ben Goodrich]

On Thursday, May 2, 2013 6:56:15 PM UTC-4, Andrew Gelman wrote:

Hi all (especially Bob).  During the times we've discussed discrete sampling, one possibility we've talked about is to work with latent continuous variables.  

I talked about constructions like this in my missing data class. At the moment, it is pretty painful to implement in Stan code because there is no support for passing a vector of lower and upper bounds to the latent parameter, which would be a nice feature even if not facing discrete unknowns.

Ben

[Andrew Gelman]

If you wanted to bound the discrete variables, that could be done directly in the model.

In any case, I'm sure there are examples where this would be slow.  What I like about my proposal is that it could all be done internally in Stan, I think.  That is, instead of us having to recommend that users apply "the continuous variable trick," we could just let them specify discrete parameters and Stan would do the sampling under the hood.

Of course, we'd want to make it clear in the manual that, as much as possible, users should integrate out latent parameters so the program will run faster.

[Ben Goodrich]

Seems hard to do it fully internally. In any event, it would be easier to play around with this idea if you could pass vectors of lower and upper bounds to the declaration of z. Also, models with latent variables determining a nominal are a further difficulty because there would have to be some way to declare the constraint that one of the latents is bigger than the others for each observation but which one is the biggest varies across observations. There isn't an easy way to do that in Stan currently.

Ben

[Ben Goodrich]

On Thu, May 2, 2013 at 10:10 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

I think we could just write the models as is, but we can see.

It is reasonably easy to internally map some things like

binary ~ bernoulli(inv_logit(mu));

to

z ~ logistic(mu, 1)[lower,upper];

But Stan would need a huge amount of self-awareness to map

binary ~ bernoulli(1.0 / (1.0 + exp(-mu));

to

z ~ logistic(mu, 1)[lower,upper];

Ben

[Andrew Gelman]

Ben,

No, what you're saying is different from what I'm saying.  In your first example, I would _keep_ the line

binary ~ bernoulli(inv_logit(mu));

and then just augment it with

 real z[N];
 simplex[K] p;

 real<lower=0,upper=1> cutpoints[K]
 cutpoints <- cumulative_sum ℗

etc.

A

[Ben Goodrich]

On Thursday, May 2, 2013 10:29:59 PM UTC-4, Andrew Gelman wrote:

Ben,

No, what you're saying is different from what I'm saying.  In your first example, I would _keep_ the line

binary ~ bernoulli(inv_logit(mu));

and then just augment it with

 real z[N];
 simplex[K] p;

 real<lower=0,upper=1> cutpoints[K]
 cutpoints <- cumulative_sum ℗

etc.

I get it, but what do you augment

binary ~ bernoulli(1.0 / (1.0 + exp(-mu));

with?

Ben
 

[Andrew Gelman]

I'm assuming "binary" is a scalar here.  Then you'd augment 

binary ~ bernoulli(inv_logit(mu));

with

real<lower=0,upper=1> z;

real<lower=0,upper=1> p;

binary <- int_step (z-p)

[Ben Goodrich]

You are answering your own question rather than the one I was raising. Ideally, you would want to map

binary ~ bernoulli(1.0 / (1.0 + exp(-mu));

to the same thing that you map

binary ~ bernoulli(inv_logit(mu));

to, but that would entail a massive amount of pattern recognition in the Stan parser. So, it wouldn't be feasible to support and internalize all the Stan syntax that would be legal if binary were fixed, which is why I think we would first need to make it possible for the users to do latent variable constructions in their .stan files and that would entail being able to pass vectors of lower and upper bounds. With that infrastructure, it would be possible to internalize a finite number of common constructions like

binary ~ bernoulli(inv_logit(mu));

or

binary ~ bernoulli(Phi(mu));

Ben

[Andrew Gelman]

I would augment

binary ~ bernoulli(1.0 / (1.0 + exp(-mu));

with the same thing:

real<lower=0,upper=1> z;
real<lower=0,upper=1> p;
binary <- int_step (z-p)

This would still work, no?  I'm not expecting Stan to parse "1.0 / (1.0 + exp(-mu))" to "invlogit" (also, I think/hope it would be ok to use "1" rather than "1.0" in this expression!).  But that's the least of our worries.  We make it clear in the manual that it's best to use canned forms where they exist.

[Ben Goodrich]

How would Stan know to use a truncated logistic distribution for the latent variable whose observation mechanism determines binary if the user writes

binary ~ bernoulli(1.0 / (1.0 + exp(-mu));

unless Stan pattern recognizes that as the standard logistic CDF? Or it could be some other CDF or a CDF parameterized a different way, or a quasilikelihood thing where the variance depends on the location instead of being fixed to 1.

I'm just saying it would be really hard to be able to internalize syntax that is arbitrary but legal if binary were fixed, so you would want to make that possible by allowing users to do the latent variable constructions themselves, even if most of the time they were using common constructions that could be internalized.

Ben

[Andrew Gelman]

I don't think Stan would need to use a truncated logistic distribution in this case.  Couldn't it just define z and p and go from there?

[Bob Carpenter]

Technically, I don't think it would be that hard to add
vectors of lower and upper bounds. They're not too hard
to do by hand if they're in vectors, but not something we
probably want users messing with on a regular basis.

- Bob

[Ben Goodrich]

On Thursday, May 2, 2013 11:12:02 PM UTC-4, Andrew Gelman wrote:

I don't think Stan would need to use a truncated logistic distribution in this case.  Couldn't it just define z and p and go from there?

I'm not seeing how you get the right posterior distribution without it. I was thinking for example of a model of a binary outcome that has missingness and estimating the same model that MCMCpack does with MCMCprobit(), except without the Gibbs sampling. Something like (hypothetically)

data {
  int<lower=1> K;
  int<lower=K> N;
  matrix[N,K]  X;
  int y[n];
  int<lower=0,upper=1> R[N];
}
transformed data {
  vector[N] lower;
  vector[N] upper;
  int n_miss;
  n_miss <- N - sum(R);
  for(i in 1:N) {
    if(R[i] == 1) {
      lower[i] <- 0.0;
      upper[i] <- Infinity();
    }
    else {
      lower[i] <- -Infinity();
      upper[i] <- 0.0;
    }
  }
}
parameters {
  vector<lower=lower,upper=upper>[N] z; # does not work now
  vector[K] beta;
}
model {
  z ~ normal(X * beta, 1)[lower,upper];
  beta ~ cauchy(0,5);
  // maybe a likelihood for R too
}

In that case, you could in principle do the internalization that you are talking about to simplify the syntax, but you can't get away from the fact that the likelihood for z is truncated normal, can you? And thus if you were to do such an internalization, Stan would have to be able to deduce what truncated distribution to use, which is doable for a small number of common constructions but hard for arbitrary syntax.

As it stands now, you can do a model like that but it basically entails two changes-of-variables to circumvent the lack of support for passing vectors of bounds for z.

Ben
 

[Ben Goodrich]

On Thursday, May 2, 2013 11:30:07 PM UTC-4, Ben Goodrich wrote:

I meant to write

transformed data {
  vector[N] lower;
  vector[N] upper;
  int n_miss;
  n_miss <- N - sum(R);
  for(i in 1:N) {


    if(R[i] == 1 && y[i] == 1) {


      lower[i] <- 0.0;
      upper[i] <- Infinity();
    }


    else if(R[i] == 1 && y[i] == 0) {


      lower[i] <- -Infinity();
      upper[i] <- 0.0;
    }


    else { // missing y[i]
      lower[i] <- -Infinity();
      upper[i] <-  Infinity();
    }
  }
}


[Andrew Gelman]

Ben

I'll have to take a look, but my first thought is that you may be making it more difficult than you need to be.  You're expressing the model using these truncated distributions (the way that Bugs does latent variables).  I think it will be ok to express the model using untruncated distributions.  As I see it, the full joint distribution is untruncated; the truncation comes conditional on the discrete variable, but that will come out naturally in the model.  I could be missing something, though (of course).

A

[Ben Goodrich]

On Thursday, May 2, 2013 11:27:34 PM UTC-4, Bob Carpenter wrote:

Technically, I don't think it would be that hard to add
vectors of lower and upper bounds.  They're not too hard
to do by hand if they're in vectors, but not something we
probably want users messing with on a regular basis.

Not too hard to do by hand but possibly involving a change-of-variables. The other thing that would be nice but hard that I mentioned in another post is something like

parameters {
  matrix<vector_of_restrictions>[N,3] Z;
}

where vector_of_restrictions indicates which of the 3 elements in each row of Z should be the largest. That is necessary for the latent variable formulation of the multinomial probit model. You can do it now with something like

parameters {
  matrix[N,3] Z;
}
model {
  matrix[N,3] ZZ;
  ZZ <- Z;
  for(i in 1:N) {
    if(y[i] == 1) {
      if(Z[i,2] > Z[i,3]) {
        ZZ[i,1] <- exp(Z[i,1]) + Z[i,2];
      }
      else {
        ZZ[i,1] <- exp(Z[i,1]) + Z[i,3];
      }
      lp__ <- lp__ + Z[i,1]; // Jacobian
    }
   if(y[i] == 2) {
     ...
   }
   if(y[i] == 3) {
     ...
   }
  }
}

but that is a mess, particularly for an arbitrary number of columns in Z instead of 3.

Ben

[Ben Goodrich]

On Thursday, May 2, 2013 11:36:13 PM UTC-4, Andrew Gelman wrote:

Ben

I'll have to take a look, but my first thought is that you may be making it more difficult than you need to be.  You're expressing the model using these truncated distributions (the way that Bugs does latent variables).  I think it will be ok to express the model using untruncated distributions.  As I see it, the full joint distribution is untruncated; the truncation comes conditional on the discrete variable, but that will come out naturally in the model.  I could be missing something, though (of course).

A

I don't think I've ever seen the posterior distribution that MCMCprobit() draws from written with untruncated distributions. Isn't it the same as any other situation like if theta has restricted support and the user does

theta ~ normal(x * beta,1);

In that case, the amount of mass that the restriction on theta precludes depends on the unknown beta so you have to truncate explicitly to the interval of theta's support. In the missing data example I posted, I was just using a vector z instead of a scalar theta and cutting at 0.

Ben
 

[Andrew Gelman]

Yes, but, in my formulation, z has unrestricted support, its uniformly distributed on (0,1).  Conditional on the discrete variable, z is restricted.  But that's an issue only with Gibbs, not with HMC, I think?

[Ben Goodrich]

On Friday, May 3, 2013 12:05:28 AM UTC-4, Andrew Gelman wrote:

Yes, but, in my formulation, z has unrestricted support, its uniformly distributed on (0,1).  

And just let it reject proposals for z outside of (0,1)? We usually tell users not to do that. I guess, going back to your original post, it isn't clear to me what is the functional form by which z is entering the log-posterior. What hypothetical .stan file would implement something like a probit model with missingness on the outcome?

Ben

[Andrew Gelman]

No need to reject z outside of (0,1).  Stan automatically does a logit transformation on parameters bounded at both ends.

[Ben Goodrich]

On Friday, May 3, 2013 12:29:13 AM UTC-4, Andrew Gelman wrote:

No need to reject z outside of (0,1).  Stan automatically does a logit transformation on parameters bounded at both ends.

Okay, that is what I think we have been calling restricted support, like vector<lower=0,upper=1>[N] z does the logit transformation for you. Regardless, I'm still not understanding how z would enter the log-posterior in a probit model with missingness on the outcome.

Ben
 

[Matt Hoffman]

My guess is it'll be significantly worse than JAGS, for a few reasons:
it's throwing away the advantages of HMC, it's not exploiting the
conditional independence structure of the model, and it's imposing an
order on variables that may not have an intrinsic order (e.g., to get
to a=5 from a=1 you have to go through a=2–4, any or all of which may
be arbitrarily unlikely in, say, a mixture model).

But it's worth testing! It may be that there are some scenarios that
it works ok for. (In particular, it's worth testing before we get too
deep into the language issues involved in doing it automatically.)

Matt

[Bob Carpenter]

On 5/2/13 6:56 PM, Andrew Gelman wrote:
> Hi all (especially Bob). During the times we've discussed discrete sampling, one possibility we've talked about is to work with latent continuous variables. We've talked about this as a possibility but I always had the sense that it might be tricky to get the distributions right.
>

> But, thinking about it, I think it's close to trivial. For simplicity, suppose our model has only two parameters: a N-dimensional discrete parameter a (which can take on the values 1,2,�,K) and an M-dimensional continuous parameter b. We can define a new continuous parameter z defined on [0,1], along with K probabilities p (i.e., simplex[K] p).

[Bob Carpenter]

Sorry for the empty message. What I meant was,
"what Matt said". The step functions essentially
cut off the derivative information from propagating
through the step function. So if you use the resulting
discrete parameter, information from its use (such
as as a mixture component), doesn't get back.

- Bob

[Michael Betancourt]

Ha, this is someone close to the approach that Matt had mentioned (although I'm coming from the physics perspective, and I'm trying to maintain the multivariate structure instead of shoving everything into a big simplex as done here).  The discontinuities in the likelihood (and hence posterior) are indeed the problem -- now you can use some tricks to build an integrator capable of handling the discontinuities but it's not really feasible for that integrator to also work for other continuous parameters.  It seems like you're always going to have to Gibbs between continuous and discrete variables.

Mike

On Fri, May 3, 2013 at 9:05 AM, Bob Carpenter <ca...@alias-i.com> wrote:

Sorry for the empty message.  What I meant was,
"what Matt said".  The step functions essentially
cut off the derivative information from propagating
through the step function.  So if you use the resulting
discrete parameter, information from its use (such
as as a mixture component), doesn't get back.

- Bob

On 5/3/13 2:39 AM, Bob Carpenter wrote:

On 5/2/13 6:56 PM, Andrew Gelman wrote:

Hi all (especially Bob).  During the times we've discussed discrete sampling, one possibility we've talked about is to
work with latent continuous variables.  We've talked about this as a possibility but I always had the sense that it
might be tricky to get the distributions right.

But, thinking about it, I think it's close to trivial.  For simplicity, suppose our model has only two parameters:  a

N-dimensional discrete parameter a (which can take on the values 1,2,…,K) and an M-dimensional continuous parameter

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Marcus Brubaker]

Because of the loss of gradient information in the continuous
representation of the discrete parameters, this is effectively doing
something like Gauss-Metropolis in those variables while doing
HMC/NUTS in the regular continuous variables (although the Hamiltonian
in that case would be very badly behaved). The performance of the
sampler would likely be limited by the performance of the Metropolis
portion though, so this will be a huge drag on the sampling. (Plus,
if you're going to loose gradient info anyway, you might as well do
real x; int y; y<- round(x); I have a hard time seeing a huge
advantage in other parameterizations.)

I suspect that the few direct discrete samplers that would be
relatively easy to implement in Stan (Slice, Gibbs, etc) would be more
efficient than this.

Cheers,
Marcus

>> email to stan-dev+u...@googlegroups.com.

>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>>
>

> --
> You received this message because you are subscribed to the Google Groups
> "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an

> email to stan-dev+u...@googlegroups.com.

[Andrew Gelman]

I guess the question is, how slow will it be. I think it would be much easier to implement than the alternatives.

[Nigel/essence]

I'm surprised nobody has mentioned DREAMZ, athough I know several of you know about it.

I have been using DREAMZ in combination with NUTS for some time, works a dream ha ha.

My current work has been looking at SMC/particles/importance sampling/resampling/move with NUTS. Should be done in a few days.

BTW, I find NUTS to be very reliable, it has been used more or less full time for the last few months.

[Bob Carpenter]

What's DREAMZ? I'm guessing DREAM for integers,
but is there a reference?

I tried searching, but only found an earlier e-mail
from you on our list!

- Bob

On 5/4/13 12:41 PM, Nigel/essence wrote:
> I'm surprised nobody has mentioned DREAMZ, athough I know several of you know about it.
>
> I have been using DREAMZ in combination with NUTS for some time, works a dream ha ha.
>
> My current work has been looking at SMC/particles/importance sampling/resampling/move with NUTS. Should be done in a few days.
>
> BTW, I find NUTS to be very reliable, it has been used more or less full time for the last few months.
>
> On Friday, 3 May 2013 21:20:56 UTC+1, Andrew Gelman wrote:
>> I guess the question is, how slow will it be. I think it would be much easier to implement than the alternatives.
>>
>>
>>
>> On May 3, 2013, at 1:53 AM, Matt Hoffman wrote:
>>
>>
>>
>>> My guess is it'll be significantly worse than JAGS, for a few reasons:
>>
>>> it's throwing away the advantages of HMC, it's not exploiting the
>>
>>> conditional independence structure of the model, and it's imposing an
>>
>>> order on variables that may not have an intrinsic order (e.g., to get
>>

>>> to a=5 from a=1 you have to go through a=2�4, any or all of which may

>>>>>>>> But, thinking about it, I think it's close to trivial. For simplicity, suppose our model has only two parameters: a N-dimensional discrete parameter a (which can take on the values 1,2,�,K) and an M-dimensional continuous parameter b. We can define a new continuous parameter z defined on [0,1], along with K probabilities p (i.e., simplex[K] p).

[Nigel/essence]

Slightly wrong, I mean DREAM(D), there is another one DREAM(ZS). They are all Vrugt and Ter Braak etc.

Don't ask me whether they really fulfill all MCMC criteria! Proofs seem sketchy, as always with these authors, but they seem to work which is the main thing, I haven't come across a practical alternative.

Just completed a study, 60 parameters, 20 discrete, some of the discrete have 3 discrete values, some have 50 discrete values. Got a good variation in results, which pleased my client. I apply NUTS to the non-discrete parameters. It's no use trying to get NUTS to jump between discrete - it's the same problem as jumping between modes, it won't happen.

Got to be careful when hitting bounds.....I normally reflect, but reflecting with discrete needs to be treated carefully, to get correct reversible probabilities.

Hydrol. Earth Syst. Sci., 15, 3701–3713, 2011
www.hydrol-earth-syst-sci.net/15/3701/2011/
doi:10.5194/hess-15-3701-2011