Sampling from inverted beta/dirichlet

[Stephen Martin]

Hey everyone,

I'm fitting a fairly weird model (one that I've posted here in the past) that contains mixtures.

Because I want to favor categorization, I prefer sampling mixture weights (per case) from, e.g., a beta or dirichlet where the parameters are < 1. That puts more probability over 0s and 1s, and helps 'encourage' assignment to one cluster over others.

However, I noticed today that when using beta(.2,.2) as a prior, for example, that sampling was *slow* (like, 45 minutes later, and only hitting 10%, slow), and I assume that it's because so MUCH probability is over 0 and 1, and the geometry of that is awful to sample from (the negative log of that would be like an inverted U, constrained at 0 and 1).

Because of that, is it preferred to instead parameterize as log-odds and sample from a wide logistic?

That is, is it preferred to use:
vector[N] pi_logit;
pi_logit ~ logistic(pi_logit_mean, 5); pi_logit_mean ~ logistic(0,5);

as opposed to:
vector[N] pi;
pi ~ dbeta(alpha*scale,beta*scale); alpha = (pi_mean); beta = 1 - pi_mean; scale = .2; pi_mean ~ uniform(0,1)

I assume the logistic is easier to sample from, geometrically, than a beta/dirichlet that is highly pointed near bounded areas (0 and 1). Am I correct in thinking about this?

[Stephen Martin]

I'll try this out tomorrow, but I think I can solve my sampling problem by using a MVN approximation to logits.

Using logits rather than raw probabilities is easy enough when using 2 groups, but harder with > 2 groups.
If dirichlet, like beta, can't be sampled well due to probability near constraints (and thus, basically, walls), I wonder if I can:
Use softmax regression to estimate logits of groups, then place a MVN(logits[n], 5*1.6*Identity) prior on that as an approximation to a dirichlet.

That is, rather than using: pi[n,1:K] ~ dirichlet(somePredicted/Hierarchical parameters), I could use pi_logit[n,1:K] ~ multinormal(somePredicted/Hierarchical parameters, scaling to adjust multinormal akin to using N(0,1.6) as logistic approx).

I'll play with this tomorrow.

[Daniel Lee]

Hi Stephen,

It's actually very hard to follow without the full model. I'd suggest posting a simple example with the full Stan program if you want specific help. 

Questions:

- why not simplex type for pi?

- I don't understand how you can put either a beta or a dirichlet prior on the mixture weights. They apply to different parameter types (real vs vector), so could you explain what you did? (providing code clarifies something like this)

- Stan doesn't have beta, so either you haven't actually tried it or you're giving non-working code. Could you specify what you've tried? And what you've intended to try, but haven't actually done the work?

- does it work with simulated data?

Regarding geometry: one thing that comes up often is that at first glance you might think that the geometry is determined by the parameters and the model, without data. The data + model together determine the shape of the joint probability distribution function as a function of the pareters. That's why we have this centered vs non-centered parameterization. It's the data in combination with the model that makes the geometry work or not. (I'm sure there are models where no matter what data is provided it won't work.) If you want to start talking about the geometry in depth, you'll have to not only provide the program, but also some data.

Hopefully you're able to provide a working Stan program or two (or three?) and then someone can help with your program. Without them, everything's a wild guess.

Daniel

--
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.

[Daniel Lee]

I think you want to look at the simplex type. 

On Sep 30, 2016, at 12:39 AM, Stephen Martin <hwki...@gmail.com> wrote:

I'll try this out tomorrow, but I think I can solve my sampling problem by using a MVN approximation to logits.

Using logits rather than raw probabilities is easy enough when using 2 groups, but harder with > 2 groups.

# of groups doesn't really matter as long as it's not massive.

If dirichlet, like beta, can't be sampled well due to probability near constraints (and thus, basically, walls), I wonder if I can:
Use softmax regression to estimate logits of groups, then place a MVN(logits[n], 5*1.6*Identity) prior on that as an approximation to a dirichlet.

What do you mean by "sampled"? I'm not understanding your terminology. In Stan, you're specifying the joint log probability distribution function. MCMC is applies to the log joint distribution to generate correlated draws from the posterior. At no point is a prior "sampled" from and used for rejection sampling. 

That is, rather than using: pi[n,1:K] ~ dirichlet(somePredicted/Hierarchical parameters), I could use pi_logit[n,1:K] ~ multinormal(somePredicted/Hierarchical parameters, scaling to adjust multinormal akin to using N(0,1.6) as logistic approx).

I'll play with this tomorrow.

On Thursday, September 29, 2016 at 10:33:14 PM UTC-5, Stephen Martin wrote:

Hey everyone,

I'm fitting a fairly weird model (one that I've posted here in the past) that contains mixtures.

Because I want to favor categorization, I prefer sampling mixture weights (per case) from, e.g., a beta or dirichlet where the parameters are < 1. That puts more probability over 0s and 1s, and helps 'encourage' assignment to one cluster over others.

However, I noticed today that when using beta(.2,.2) as a prior, for example, that sampling was *slow* (like, 45 minutes later, and only hitting 10%, slow), and I assume that it's because so MUCH probability is over 0 and 1, and the geometry of that is awful to sample from (the negative log of that would be like an inverted U, constrained at 0 and 1).

Because of that, is it preferred to instead parameterize as log-odds and sample from a wide logistic?

That is, is it preferred to use:
vector[N] pi_logit;
pi_logit ~ logistic(pi_logit_mean, 5); pi_logit_mean ~ logistic(0,5);

as opposed to:
vector[N] pi;
pi ~ dbeta(alpha*scale,beta*scale); alpha = (pi_mean); beta = 1 - pi_mean; scale = .2; pi_mean ~ uniform(0,1)

I assume the logistic is easier to sample from, geometrically, than a beta/dirichlet that is highly pointed near bounded areas (0 and 1). Am I correct in thinking about this?

--

[Stephen Martin]

Well, the stan model I'm using is fairly massive and novel, and would only obfuscate things if I posted the full thing here. Tomorrow I'll make a smaller, toy example of the aspect of the model that's causing me issues. Really, the part of the model that is currently concerning me is a mixture modeling component.

Just some points
- I do use simplex for pi when K > 2, otherwise I use real<lower=0,upper=1>
- I'm not sure I follow. For mixtures, I assume each person has a probability of belonging to cluster 1 or 2 (for k = 2, vector<lower=0,upper=1>[N] pi), pi ~ beta(alpha, beta); alpha = mean*scale, beta = (1 - mean)*scale; scale = .2; mean ~ uniform). This gives me an overall mean, and allows individuals to individually load on each cluster with some probability pi. The scale essentially inverts the beta and places greater probability on 0 or 1, thus informing the model to try 'harder' to assign people to one cluster or another.
- Stan does have the beta distribution. I've used beta several times. Hell, my model actually works, I'm just trying to extend it to include a softmax regression within it, so you can predict group membership probabilities, but also try and encourage high probability clustering.
- haven't tried simulating data /from/ it, but have fit the model *to* simulated data.

[Stephen Martin]

I understand that stan samples from the joint posterior.

I'm simply stating that my priors are purposefully placing greater probability toward constrained areas to encourage clustering (encouraging probability values of 0 and 1), which likely causes the posterior to be difficult to sample from, as the -log of that creates an extremely steep "hill" of sorts. To avoid that, I can reparameterize in terms of logits (which aren't constrained), and use a prior that is fairly wide so that larger logit values are more probable.
Extending that to a multigroup case is the difficulty. beta -> dirichlet, as logit -> ?; I'm thinking I could use a MVN approximation to a 'multivariate logistic' distribution of sorts.

[Daniel Lee]

On Sep 30, 2016, at 12:52 AM, Stephen Martin <hwki...@gmail.com> wrote:

Well, the stan model I'm using is fairly massive and novel, and would only obfuscate things if I posted the full thing here. Tomorrow I'll make a smaller, toy example of the aspect of the model that's causing me issues. Really, the part of the model that is currently concerning me is a mixture modeling component.

Great. That'll help. Your points below didn't clarify enough for me to follow. 

Just some points
- I do use simplex for pi when K > 2, otherwise I use real<lower=0,upper=1>

So new info (that wasn't present in the original email): K = number of components, N = observations. 

- I'm not sure I follow. For mixtures, I assume each person has a probability of belonging to cluster 1 or 2 (for k = 2, vector<lower=0,upper=1>[N] pi), pi ~ beta(alpha, beta); alpha = mean*scale, beta = (1 - mean)*scale; scale = .2; mean ~ uniform). This gives me an overall mean, and allows individuals to individually load on each cluster with some probability pi. The scale essentially inverts the beta and places greater probability on 0 or 1, thus informing the model to try 'harder' to assign people to one cluster or another.

Is there 1 alpha? Or N?

What's motivation for tying alpha and beta together in that way? And when you say individually, does that really mean you have N means? Or N alphas? Or N betas? 

What is an overall mean? Is it estimated? Transformed data? A weighted average? (sorry for being pedantic, but you're using words that I understand, but I really don't understand how you're using them. Math written carefully is unambiguous. And a Stan program is also unambiguous. But the words alone aren't pinpointing your meaning [to me])  

- Stan does have the beta distribution. I've used beta several times. Hell, my model actually works, I'm just trying to extend it to include a softmax regression within it, so you can predict group membership probabilities, but also try and encourage high probability clustering.
- haven't tried simulating data /from/ it, but have fit the model *to* simulated data.

Sorry, maybe my phone's autocorrect corrected it. In your original email you wrote "dbeta." That doesn't exist in Stan, unless you wrote your own function. 

Looking forward to a Stan program that would get us talking the same language. 

[Stephen Martin]

I'll write up a quick mixture model in stan tomorrow to show the issue because, yes, that would clarify some things.
As for the dbeta thing, it's an old habit. I use beta; dbeta is a carryover from R and I just accidentally typed it.

I tie alpha and beta together that way because it's essentially parameterizing beta in terms of a mean and scale instead of the typical 'counts' parameterization.

[Michael Betancourt]

Please read the manual about constrained parameters such

as the simplex.  We implicitly do a similar transformation to

an unconstrained space already.  If you want to see why the

sampler is slow then use CmdStan with the “diagnostic_file”

option to get traces of the unconstrained parameters.

[Stephen Martin]

Alright. Here is the code:

data {
    int N;
    vector[N] x;
    int K;
    real dirichlet_scale;
}
parameters{
    ordered[K] means;
    vector<lower=0>[K] sds;
//    simplex[K] pi[N];
    vector<lower=0>[K] pi_logit[N];
    vector<lower=0>[K] beta0;
}
transformed parameters{
    simplex[K] pi_mean;
    simplex[K] pi[N];
    pi_mean = softmax(beta0);
    for(n in 1:N){
        pi[n] = softmax(pi_logit[n]);
    }
}
model{
    vector[K] lp_nk[N];
    vector[N] lp_n;
   
    means ~ normal(0,4);
    sds ~ cauchy(0,3);
    beta0 ~ logistic(0,1);
    for(n in 1:N){
        for(k in 1:K){
            pi_logit[n,k] ~ logistic(beta0[k],3);
        }
//        pi[n] ~ dirichlet(pi_mean*dirichlet_scale);
    }
   
    for(n in 1:N){
        for(k in 1:K){
            lp_nk[n,k] = log(pi[n,k]) + normal_lpdf(x[n] | means[k], sds[k]);
        }
        lp_n[n] = log_sum_exp(lp_nk[n]);
    }
target += sum(lp_n);
}

The goal of the model (in this toy example) is to fit a mixture model, but with the ability to predict latent class memberships.
In this particular case, it's an intercept-only model, so each person's "predicted membership" is given by beta0.

You can use softmax regression to predict logits; I use priors to constrain the intercepts to be positive and logistic-distributed, because softmax regression needs some constraining (several equal solutions otherwise).
In the dirichlet method, the logits are transformed to probabilities, and each individual's latent membership probabilities (pi) is assumed to be dirichlet-distributed with the predicted probabilities as the hyperparameter; the scale forces probability to be allocated toward 0s and 1s, and so encourages loading on one cluster. The geometry, however, is extremely weird, because it creates a very steep inverted-U that causes difficulty sampling and divergent transitions (to see what it could do, look in R at curve(dbeta(x,.4*.2,.6*.2,log=TRUE)*-1,0,1) ). I'm guessing the posterior has a ton of probability right at the bounded regions and it's creating difficulties.

Another method keeps using logits.
An individual's logits are assumed to be multivariate-logistic distributed with mean-logits at the predicted logits; because the assumed multi-variate-logistic has no covariance, it's the same as using univariate logistics each with a beta0 mean. To again encourage high probabilities for certain clusters, we let the scale be large (and so one can easily have very high logits).
Using this method, no divergent transitions occur, and the model works.

So it seems like I may have answered my own question --- Sampling from inverted-betas and inverted-dirichlets makes the posterior geometry difficult to sample from, it's slow and causes divergent transitions. Using an alternative method that uses only logit-sampling works fine, and it seems to be because the geometry is more forgiving.

On Thursday, September 29, 2016 at 10:33:14 PM UTC-5, Stephen Martin wrote:

[Bob Carpenter]

Thanks so much for following up with all this detail. That's
exactly how we go about trying to diagnose sampling problems like
this.

Another problem you run into is that if you have K degrees of
freedom going into beta0 and run softmax(beta0), then you have
identifiability problems. You can pin one component of beta0 to
0 for identifiability, but that's awkward because now beta0 is
asymmetric.

If you declare a variable as

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

then beta gets a uniform distribution on (0, 1), which is equivalent to
writing:

real beta_raw;

beta_raw ~ logistic(0, 1);

beta = inv_logit(beta_raw);

In fact, if you work through the Jacobian of the inverse of the logit
transform going on under the hood you get the logistic (inv_logit is
cdf of the unit logitistc and its derivative is the logistic(0, 1) density).

- Bob

[Stephen Martin]

Yeah this seems to generally be an issue with softmax regression. It's not identified, and there's no base group. It's really about how big each one is relative to all others.

I finally finished out the softmax model (with predictors) and it worked reasonably well. Constraining the beta0s to be positive-logistic distributed, and identifying latent classes with mean ordering seemed to work.
The only "problem" is that even with priors on the intercepts, it can have large or small predicted logits --- However, this matters exactly none, because once you convert to probabilities, they are all equivalent.
You're right though --- If I constrained one group to have a certain logit, then the rest of them will be scaled appropriately.
For the purposes of a multinomial logistic model predicting latent memberships though, the assigned cluster probabilities will remain the same. A bit of label switching here and there, but marginally it seems to work.

I generated a mixture of gaussians such that N1=25,N2=75,N3=50; mean_1 = 0, mean_2 = 7, mean_3 = 4; I don't recall the sd's. However, the model recovered the means and individuals' cluster probabilities strikingly well.

[Stephen Martin]

For the sake of posterity, here's my [really messy] implementation of multinomial latent class regression:
Two weirds thing in this is that: using a logistic(0,1) vs a logistic(0,3) for softmax really shouldn't make a difference, since it's all relative. However, using 3 seems to still favor high-probability categorization.
Second, Rather than defining one's latent membership probability as a linear function, I still say it's "predicted" from a linear function with a standard error of 1 logit. Either way works, but allowing individuals to differ from the line can improve clustering accuracy (because it permits two individuals with the same betas and predictors to differ). Too large an error (e.g., 5 standard deviations), and the predictors are meaningless and contribute nothing to clustering. So if you want high accuracy clustering and also to estimate how predictors relate to cluster-membership, it seems like setting an sd=1 is reasonable.

data {
    int N;
    int J;
    int K;
    vector[N] x;
    matrix[N,J] predictor;

}
parameters{
    ordered[K] means;
    vector<lower=0>[K] sds;

    vector<lower=0>[K] pi_logit[N];
    vector<lower=0>[K] beta0;

    matrix[J,K] beta;
}
transformed parameters{
    simplex[K] pi[N];
//    vector[K] pi_logit[N];

    for(n in 1:N){
        pi[n] = softmax(pi_logit[n]);
    }
}
model{
    vector[K] lp_nk[N];
    vector[N] lp_n;
   
    means ~ normal(0,4);
    sds ~ cauchy(0,3);

    beta0 ~ logistic(0,3);
    for(k in 1:K){
        beta[1:J,k] ~ logistic(0,3);

    }
    for(n in 1:N){
        for(k in 1:K){

            pi_logit[n,k] ~ logistic(beta0[k] + predictor[n]*beta[1:J,k],1);

        }
    }
   
    for(n in 1:N){
        for(k in 1:K){
            lp_nk[n,k] = log(pi[n,k]) + normal_lpdf(x[n] | means[k], sds[k]);
        }
        lp_n[n] = log_sum_exp(lp_nk[n]);
    }
target += sum(lp_n);
}

On Thursday, September 29, 2016 at 10:33:14 PM UTC-5, Stephen Martin wrote:

[Stephen Martin]

Yet another update, but also a bit of a new question.

As per Bob's suggestion, I constrained beta0[1] (arbitrarily) to 0, so that the logits have some sort of scale. I then allowed the beta0[2:K] to be negative or positive (rather than constrained to be positive).
This led to divergent transitions and very wrong, albeit consistent estimates of beta0[2:K] and the latent class means/sds.

Honestly, I have NO idea why that would be the case.

To clarify (and review):
One model uses softmax regression, but constrains the intercept-logits (beta0) to be positive and ~logistic(0,3); with this parametization, several combinations of beta0s are possible, but it does not affect the model because once softmax'd, the probabilities are the same.
The other model uses softmax, but constrains beta0[1] to be 0; the others are ~ logistic(0,3) with no constraint. This should yield the same model as above (intuitively, to me, anyway), but sets a scale for which other logits to be compared to (no more multiple combinations of beta0s should be possible). This latter model suffers from div. transitions and incorrect parameter estimates.

No idea why, and I welcome any input about that.

On Thursday, September 29, 2016 at 10:33:14 PM UTC-5, Stephen Martin wrote:

[Stephen Martin]

Ah, actually, I wonder if it's because the priors are *too* wide. By allowing the other intercepts to be fairly large (logits of -3 and 3), it may be telling the model that K > 2 clusters are way more probable, and so it's biasing the model toward one dominated by 2 latent clusters... Interesting.

[Stephen Martin]

Nope, not the case. Just lowered the logits and still tons of div transitions and wrong estimates. No idea why this parameterization causes wrong estimates and div transitions.

[Bob Carpenter]

> On Oct 3, 2016, at 3:28 PM, Stephen Martin <hwki...@gmail.com> wrote:
>
> Yet another update, but also a bit of a new question.
>
> As per Bob's suggestion, I constrained beta0[1] (arbitrarily) to 0, so that the logits have some sort of scale. I then allowed the beta0[2:K] to be negative or positive (rather than constrained to be positive).
> This led to divergent transitions and very wrong, albeit consistent estimates of beta0[2:K] and the latent class means/sds.
>
> Honestly, I have NO idea why that would be the case.

It can come about by having a free-floating parameter that's
not being used anywhere and not given a prior. Make sure that's
not going on, first.

- Bob

>
> To clarify (and review):
> One model uses softmax regression, but constrains the intercept-logits (beta0) to be positive and ~logistic(0,3); with this parametization, several combinations of beta0s are possible, but it does not affect the model because once softmax'd, the probabilities are the same.
> The other model uses softmax, but constrains beta0[1] to be 0; the others are ~ logistic(0,3) with no constraint. This should yield the same model as above (intuitively, to me, anyway), but sets a scale for which other logits to be compared to (no more multiple combinations of beta0s should be possible). This latter model suffers from div. transitions and incorrect parameter estimates.
>
> No idea why, and I welcome any input about that.
>
> On Thursday, September 29, 2016 at 10:33:14 PM UTC-5, Stephen Martin wrote:
> Hey everyone,
>
> I'm fitting a fairly weird model (one that I've posted here in the past) that contains mixtures.
>
> Because I want to favor categorization, I prefer sampling mixture weights (per case) from, e.g., a beta or dirichlet where the parameters are < 1. That puts more probability over 0s and 1s, and helps 'encourage' assignment to one cluster over others.
>
> However, I noticed today that when using beta(.2,.2) as a prior, for example, that sampling was *slow* (like, 45 minutes later, and only hitting 10%, slow), and I assume that it's because so MUCH probability is over 0 and 1, and the geometry of that is awful to sample from (the negative log of that would be like an inverted U, constrained at 0 and 1).
>
> Because of that, is it preferred to instead parameterize as log-odds and sample from a wide logistic?
>
> That is, is it preferred to use:
> vector[N] pi_logit;
> pi_logit ~ logistic(pi_logit_mean, 5); pi_logit_mean ~ logistic(0,5);
>
> as opposed to:
> vector[N] pi;
> pi ~ dbeta(alpha*scale,beta*scale); alpha = (pi_mean); beta = 1 - pi_mean; scale = .2; pi_mean ~ uniform(0,1)
>
> I assume the logistic is easier to sample from, geometrically, than a beta/dirichlet that is highly pointed near bounded areas (0 and 1). Am I correct in thinking about this?
>

[Stephen Martin]

From what I can tell, all parameters are defined and with a prior. However, the problem is that the model is putting noone on cluster 1, so the cluster 1 parameters are not being estimated from the data and are purely draws from the priors (pi[n,1] is essentially 0). Again, I'm not sure why it is any different from a non-scaled simplex regression. I may play with this some more tomorrow and post the code. I just finished cleaning up the working model to be faster (22 seconds for a multinomial logistic model predicting latent class membership; not too shabby!), so I'd have to undo a lot of that to be parameterized using K-1 logits with one fixed to 0.

[Daniel Lee]

Hi Stephen,

I'm finally getting around to actually looking at the Stan program you've sent. It really does help to have the program. Thanks.

If you have the generated data too, please send that along.

If you want to know how the simplex parameter is parameterized in Stan, check the manual. It describes how we parameterize it with K-1 unconstrained parameters. (Yes, you can build the transformation on your own and it is instructive, but we have it built in and tested.)

Daniel

On Mon, Oct 3, 2016 at 9:58 PM, Stephen Martin <hwki...@gmail.com> wrote:

From what I can tell, all parameters are defined. However, the problem is that the model is putting noone on cluster 1, so the cluster 1 parameters are not being estimated from the data and are purely draws from the priors (pi[n,1] is essentially 0). Again, I'm not sure why it is any different from a non-scaled simplex regression. I may play with this some more tomorrow and post the code. I just finished cleaning up the working model to be faster (22 seconds for a multinomial logistic model predicting latent class membership; not too shabby!), so I'd have to undo a lot of that to be parameterized using K-1 logits with one fixed to 0.

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