Little mixing after messing with adapt_delta and stepsize_jitter

[Joseph Burling]

Hello,

I'm trying to figure out ways to improve my model below. I'm trying to fit an overdispersed, zero-inflated poisson multilevel model. I get very little mixing, particularly with a few of the sigma parameters below. I've tried decreasing step size by messing with adapt_delta and also playing around with different values of stepsize_jitter. Everyone once in a while I'll have a chain that hits the max tree depth for a while, which drastically increases the time after running some test samples.

Not sure if some parts need to be reparameterized or perhaps some stronger constraints. Any suggestions would be very helpful. Thanks.

data {

    int<lower=1>        N;             // number of total observations

    int<lower=0>        Pixels[N];     // pixel count data (object area from a 640x480 img)

    int<lower=1>        J;             // number of dyads observed

    int<lower=1>        I;             // number of objects used

    int<lower=1>        K;             // number of effects from design matrix

    matrix[N, K]        X;             // population predictor (design) matrix

    matrix[N, K*J]      Z_dyad;        // dyad predictor matrix

    matrix[N, K*I]      Z_obj;         // object predictor matrix

    vector<lower=0>[N]  Exposure;      // exposure (offset)

    real                sig_T;         // sigma truncation = 3 * residual sd using mean (~ 15)

    real                b_T;           // beta truncation = log(10 x image area) (~ 17)

}

transformed data {

  real                 sigma_T;        // truncated cauchy - upper limit

  vector[K]            mu_zeros;       // group means for level 2 betas (all zeros)

  sigma_T <- atan(sig_T);  

  mu_zeros <- rep_vector(0.0, K);      // used for both grouping factors since estimating same K set of effects

}

parameters {

    real<lower=-b_T, upper=b_T>            beta_mu;           // prior mean for all population beta parameters

    real<lower=0, upper=sigma_T>           sigma_pop_unif;    // prior sd for all population beta parameters

    vector<lower=-b_T, upper=b_T> [K]      beta_pop_raw;      // matt trick adjustment for each population beta parameter

    matrix<lower=-b_T, upper=b_T> [K, J]   beta_dyad_raw;     // group effect beta parameters (dyad errors)

    matrix<lower=-b_T, upper=b_T> [K, I]   beta_obj_raw;      // group effect beta parameters (object errors)

    vector<lower=0, upper=sigma_T>[K]      sigma_dyad_unif;   // group effect sd uniform (for cauchy)

    vector<lower=0, upper=sigma_T>[K]      sigma_obj_unif;    // group effect sd uniform (for cauchy)

    corr_matrix[K]                         Rho_dyad;          // dyad corr matrix

    corr_matrix[K]                         Rho_obj;           // obj corr matrix

    real<lower=0, upper=sigma_T>           epsilon_unif;      // sd epsilon (for overdispersion)

    vector<lower=-b_T, upper=b_T> [N]      epsilon_raw;       // matt trick adjust for each residual

    real<lower=0, upper=1>                 theta;             // probability of zero count

}

transformed parameters {

    matrix[K, J]                        mu_d;         // temp matrix dyad (for reshaping matrix to vec)

    matrix[K, I]                        mu_o;         // temp matrix object (for reshaping matrix to vec)

    vector[K]                           sigma_dyad;   // final dyad sd for corr matrix

    vector[K]                           sigma_obj;    // final object sd for corr matrix

    vector<lower=-b_T, upper=b_T>[K*J]  Beta_dyad;    // final dyad parameters

    vector<lower=-b_T, upper=b_T>[K*I]  Beta_object;  // final object parameters

    vector<lower=-b_T, upper=b_T>[K]    Beta_mu;      // final tranformed population Beta means

    for (k in 1:K) {

        sigma_dyad[k] <- tan(sigma_dyad_unif[k]);  // cauchy reparameterize

        sigma_obj[k] <- tan(sigma_obj_unif[k]);

        mu_d[k] <- row(beta_dyad_raw, k) - mean(beta_dyad_raw[k]);

        mu_o[k] <- row(beta_obj_raw, k) - mean(beta_obj_raw[k]);

    }

    Beta_dyad <- to_vector(mu_d);

    Beta_object <- to_vector(mu_o);

    Beta_mu <- beta_mu + tan(sigma_pop_unif) * beta_pop_raw;

}

model {

    matrix[K, K]  S_d;               // diagonal sd matrix dyad

    matrix[K, K]  S_o;               // diagonal sd matrix object

    matrix[K, K]  L_d;               // chol decomp of corr matrix for dyad

    matrix[K, K]  L_o;               // chol decomp of corr matrix for obj

    vector[N]     epsilon_raw_adj;   // uncentered residuals

    vector[N]     epsilon_i;         // residuals centered at zero

    vector[N]     lambda;            // predicted count

    // implicit: Rho_dyad ~ lkj_corr(1.0)

    // implicit: Rho_obj ~ lkj_corr(1.0)

    // implicit: sigma_pop_unif ~ uniform(0, sigma_T);

    // implicit: sigma_dyad_unif ~ uniform(0, sigma_T);

    // implicit: sigma_obj_unif ~ uniform(0, sigma_T);

    // implicit: epsilon_unif ~ uniform(0, sigma_T);

    // implicit: theta ~ uniform(0, 1);

    epsilon_raw ~ normal(0, 1);

    beta_mu ~ normal(0, 5);

    beta_pop_raw ~ normal(0, 1);

    S_d <- diag_matrix(sigma_dyad);

    L_d <- cholesky_decompose(Rho_dyad);

    S_o <- diag_matrix(sigma_obj);

    L_o <- cholesky_decompose(Rho_obj);

       

    for (j in 1:J) {

        col(beta_dyad_raw, j) ~ multi_normal_cholesky(mu_zeros, S_d * L_d);

    }

    for (i in 1:I) {

        col(beta_obj_raw, i) ~ multi_normal_cholesky(mu_zeros, S_o * L_o);

    }

    epsilon_raw_adj <- tan(epsilon_unif) * epsilon_raw;

    epsilon_i <- epsilon_raw_adj - mean(epsilon_raw_adj);

    lambda <- Exposure .* exp(X * Beta_mu + Z_dyad * Beta_dyad + Z_obj * Beta_object + epsilon_i);

    for (n in 1:N) {

        if (Pixels[n] == 0) {

            increment_log_prob(log_sum_exp(bernoulli_log(1, theta), bernoulli_log(0, theta) + poisson_log(Pixels[n], lambda[n])));

        } else {

            increment_log_prob(bernoulli_log(0, theta) + poisson_log(Pixels[n], lambda[n]));

        }

    }

}

generated quantities {

    real       epsilon;

    epsilon <- tan(epsilon_unif);

}

[Bob Carpenter]

I don't understand the model's parameterization as is, but
maybe someone else will have a clue. In particular, I'm not
sure what's going on with all the centering you're doing,
such as

mu_d[k] <- row(beta_dyad_raw, k) - mean(beta_dyad_raw[k]);

mu_o[k] <- row(beta_obj_raw, k) - mean(beta_obj_raw[k]);

epsilon_i <- epsilon_raw_adj - mean(epsilon_raw_adj);

Another approach I've used for this in large IRT models that
aren't inherently identified is to give the epsilon_raw_adj
a prior centered at 0 (the uniform does this, but it doesn't favor
a centered result).

Does the model at least converge to a reasonable estimate?
And does it converge to the same place from multiple starting
points? If so, does it converge well before the warmup phase
is over? Adaptation needs to spend some time after convergence
before it can estimate the scales properly.

Is there any way to start with a simpler model and work up
to the complex model you posted? That should let you both
assess whether the model's getting reasonable results and
as you add more complexity, put your finger on where the poor
mixing arises.

If there are only a few hundred parameters, you could also
try a dense mass matrix --- it could be that there's not a reasonable
way to scale all the parameters without a rotation that makes sense.

You only have uniform-in-an-interval priors on many of
the parameters --- is the range reasonable given what you know
about the problem? Adding something more peaked in the way of
priors can help with mixing. Also, initialization (assuming
you don't initialize yourself) will be centered at the center
of the intervals specified as constraints.

If you hit the maximum tree depth, that can be for two reasons.
One is that you just turned the step size down so low due to
jitter that it goes out to the max depth without making a U-turn.
The other is that you really need a small step size because of
high curvature in the posterior or because you're way out in
a tail and have to go a long way before turning around.

To make the model more efficient, I'd suggest not making a diagonal
matrix S_d, but doing the scaling in a loop that's *outside* the
multi_normal_cholesky loop.

Another efficiency issue is in the exp for lambda. This can
be problematic numerically, and we provide an alternative.
Instead of

lambda <- Exposure .* exp(X * Beta_mu + Z_dyad * Beta_dyad + Z_obj * Beta_object + epsilon_i);

and then

poisson_log(Pixels[n], lambda[n]))

you can set

lambda <- log(Exposure) + X * ....

and then use:

poisson_log_log(Pixels[n], lambda[n])

The possion_log distribution takes its parameter on the log scale and does
the exponentiation implicitly inside, which is both more efficient and
more arithmetically stable.

- 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.
> For more options, visit https://groups.google.com/groups/opt_out.

[Joseph Burling]

First of all, thank you for the quick and thorough response.

I don't understand the model's parameterization as is, but
maybe someone else will have a clue.  In particular, I'm not
sure what's going on with all the centering you're doing,
such as

   mu_d[k] <- row(beta_dyad_raw, k) - mean(beta_dyad_raw[k]);

   mu_o[k] <- row(beta_obj_raw, k) - mean(beta_obj_raw[k]);

   epsilon_i <- epsilon_raw_adj - mean(epsilon_raw_adj);

I did the centering this way because I was trying to get mean population beta parameters (Beta_mu) and the above are the group level betas centered at zero. I know I had a prior of zero to begin with but often one groups' parameters might get stuck far from zero and switch with the population beta parameters.

 

Another approach I've used for this in large IRT models that
aren't inherently identified is to give the epsilon_raw_adj
a prior centered at 0 (the uniform does this, but it doesn't favor
a centered result).

Does the model at least converge to a reasonable estimate?
And does it converge to the same place from multiple starting
points?  If so, does it converge well before the warmup phase
is over?  Adaptation needs to spend some time after convergence
before it can estimate the scales properly.

Yes, the model eventually converges to a reasonable estimate but it can take many hours, and it gets around that point before warmup ends. When printing the diagnostic file, adaption ends at 100 rows. I've initialized the chains to start at several different reasonable points (this was necessary, otherwise auto initialization fails). Some chains get stuck and some don't.

Is there any way to start with a simpler model and work up
to the complex model you posted?  That should let you both
assess whether the model's getting reasonable results and
as you add more complexity, put your finger on where the poor
mixing arises.

I've had less issues when only estimating a varying slope and intercept for each group. Once I start adding more varying effects for each group things start to slow down. 

 

If there are only a few hundred parameters, you could also
try a dense mass matrix --- it could be that there's not a reasonable
way to scale all the parameters without a rotation that makes sense.

You only have uniform-in-an-interval priors on many of
the parameters --- is the range reasonable given what you know
about the problem?  Adding something more peaked in the way of
priors can help with mixing.  Also, initialization (assuming
you don't initialize yourself) will be centered at the center
of the intervals specified as constraints.

I've fit a similar model with the lme4 package in order to obtain wide, reasonable constraints. A lot of the uniform priors with the constraints are used for a truncated half cauchy. If I were to use something more peaked, do you mean something like a beta and reparameterizing the cauchy differently like so, tan(pi()/2 * beta) ? Or did you have something else in mind.

 

If you hit the maximum tree depth, that can be for two reasons.
One is that you just turned the step size down so low due to
jitter that it goes out to the max depth without making a U-turn.
The other is that you really need a small step size because of
high curvature in the posterior or because you're way out in
a tail and have to go a long way before turning around.

To make the model more efficient, I'd suggest not making a diagonal
matrix S_d, but doing the scaling in a loop that's *outside* the
multi_normal_cholesky loop.

Sorry, can you elaborate a little further on this?

 

Another efficiency issue is in the exp for lambda.  This can
be problematic numerically, and we provide an alternative.
Instead of

   lambda <- Exposure .* exp(X * Beta_mu + Z_dyad * Beta_dyad + Z_obj * Beta_object + epsilon_i);

and then

   poisson_log(Pixels[n], lambda[n]))

you can set

   lambda <- log(Exposure) + X * ....

and then use:

   poisson_log_log(Pixels[n], lambda[n])

The possion_log distribution takes its parameter on the log scale and does
the exponentiation implicitly inside, which is both more efficient and
more arithmetically stable.

Thanks again, I'll try this out along with some of the other suggestions you've made. I really appreciate the help.

[Bob Carpenter]

On 1/20/14, 3:59 PM, Joseph Burling wrote:
>
> First of all, thank you for the quick and thorough response.
>
>
> I don't understand the model's parameterization as is, but
> maybe someone else will have a clue. In particular, I'm not
> sure what's going on with all the centering you're doing,
> such as
>
> mu_d[k] <- row(beta_dyad_raw, k) - mean(beta_dyad_raw[k]);
>
> mu_o[k] <- row(beta_obj_raw, k) - mean(beta_obj_raw[k]);
>
> epsilon_i <- epsilon_raw_adj - mean(epsilon_raw_adj);
>
>
>
> I did the centering this way because I was trying to get mean population beta parameters (Beta_mu) and the above are the
> group level betas centered at zero. I know I had a prior of zero to begin with but often one groups' parameters might
> get stuck far from zero and switch with the population beta parameters.

The real problem with this kind of centering is that the model
winds up being unidentified on the unconstrained scale (the centered
parameters are identified, but not the parameters being sampled, which are just
the transforms of the variables defined in the parameters block).

Andrew recommended the approach you took for BUGS in his regression book
with Jennifer, but it's problematic for Stan because the adaptation
and sampling happens on the unconstrained scale.

> Another approach I've used for this in large IRT models that
> aren't inherently identified is to give the epsilon_raw_adj
> a prior centered at 0 (the uniform does this, but it doesn't favor
> a centered result).
>
> Does the model at least converge to a reasonable estimate?
> And does it converge to the same place from multiple starting
> points? If so, does it converge well before the warmup phase
> is over? Adaptation needs to spend some time after convergence
> before it can estimate the scales properly.
>
>
>
> Yes, the model eventually converges to a reasonable estimate but it can take many hours, and it gets around that point
> before warmup ends. When printing the diagnostic file, adaption ends at 100 rows. I've initialized the chains to start
> at several different reasonable points (this was necessary, otherwise auto initialization fails). Some chains get stuck
> and some don't.

You really want to get to a reasonable point about halfway through adaptation.
That'll make sampling go much faster by allowing better adaptation.

Getting stuck is usually because of the requirement for detailed balance
in the sampler. If you jump to an unlikely place (or start there), it's
hard to get out again.

> Is there any way to start with a simpler model and work up
> to the complex model you posted? That should let you both
> assess whether the model's getting reasonable results and
> as you add more complexity, put your finger on where the poor
> mixing arises.
>
>
> I've had less issues when only estimating a varying slope and intercept for each group. Once I start adding more varying
> effects for each group things start to slow down.
>
> If there are only a few hundred parameters, you could also
> try a dense mass matrix --- it could be that there's not a reasonable
> way to scale all the parameters without a rotation that makes sense.
>
> You only have uniform-in-an-interval priors on many of
> the parameters --- is the range reasonable given what you know
> about the problem? Adding something more peaked in the way of
> priors can help with mixing. Also, initialization (assuming
> you don't initialize yourself) will be centered at the center
> of the intervals specified as constraints.
>
>
> I've fit a similar model with the lme4 package in order to obtain wide, reasonable constraints. A lot of the uniform
> priors with the constraints are used for a truncated half cauchy. If I were to use something more peaked, do you mean
> something like a beta and reparameterizing the cauchy differently like so, tan(pi()/2 * beta) ? Or did you have
> something else in mind.

We're going to be implementing something like blme in Stan this
year. We're having trouble finding time to do it, but it's one
of our higher priorities.

> If you hit the maximum tree depth, that can be for two reasons.
> One is that you just turned the step size down so low due to
> jitter that it goes out to the max depth without making a U-turn.
> The other is that you really need a small step size because of
> high curvature in the posterior or because you're way out in
> a tail and have to go a long way before turning around.
>
> To make the model more efficient, I'd suggest not making a diagonal
> matrix S_d, but doing the scaling in a loop that's *outside* the
> multi_normal_cholesky loop.
>
>
> Sorry, can you elaborate a little further on this?

Sure. Rather than this:

S_d <- diag_matrix(sigma_dyad);
L_d <- cholesky_decompose(Rho_dyad);

for (j in 1:J) {
col(beta_dyad_raw, j) ~ multi_normal_cholesky(mu_zeros, S_d * L_d);
}

it'd be better to do this:

L_d <- cholesky_decompose(Rho_dyad);
for (i in 1:cols(L_d))
for (j in i:cols(L_d))
L_d[i,j] <- L_d[i,j] * sigma_dyad[i];

for (j in 1:J) {

col(beta_dyad_raw, j) ~ multi_normal_cholesky(mu_zeros, L_d);
}

You should check my algebra --- I'm not very good at matrices.

The loop saves a bunch of multiplies over the S_d * L_d, and they're
only done once.

An intermediate solution would be to assign (S_d * L_d) to a temporary
matrix variable outside of the loop.

>
> Another efficiency issue is in the exp for lambda. This can
> be problematic numerically, and we provide an alternative.
> Instead of
>
> lambda <- Exposure .* exp(X * Beta_mu + Z_dyad * Beta_dyad + Z_obj * Beta_object + epsilon_i);
>
> and then
>
> poisson_log(Pixels[n], lambda[n]))
>
> you can set
>
> lambda <- log(Exposure) + X * ....
>
> and then use:
>
> poisson_log_log(Pixels[n], lambda[n])
>
> The possion_log distribution takes its parameter on the log scale and does
> the exponentiation implicitly inside, which is both more efficient and
> more arithmetically stable.
>
>
> Thanks again, I'll try this out along with some of the other suggestions you've made. I really appreciate the help.

Hopefully someone like Michael or Matt or Andrew who understand the
parameterizations and their effects better than me will have other
suggestions.

As Andrew likes to point out, hard questions on mailing lists
get fewer answers because they take a lot longer for the devs to
figure out what's going on.

So if you can break the reparameterization down into a simple question,
you might have more luck getting an answer.

- Bob

[Joseph Burling]

I was trying to look for other examples for the centering problem and come across this example from the ARM book using STAN.

https://github.com/stan-dev/stan/blob/feature/ARM/src/models/ARM/Ch.19/pilots.stan

If the sampling happens on the unconstrained scale which is unidentified here I'm assuming, are there any alternative reparameterizations for the type of model?

 

>         mu_d[k] <- row(beta_dyad_raw, k) - mean(beta_dyad_raw[k]);
>
>         mu_o[k] <- row(beta_obj_raw, k) - mean(beta_obj_raw[k]);
>
>         epsilon_i <- epsilon_raw_adj - mean(epsilon_raw_adj);

[Bob Carpenter]

Could someone who understands identification a bit better
than me help out in figuring out which of the two proposals
below lead to an identified parameterization of a centered
vector on the unconstrained scale?

I'd like to drop a suggestion into the next manual and
also think about providing a centered vector data type
in the future.

Proposal A:

parameters {
simplex[K] beta_ratios;
real<lower=0> beta_scale;
...
transformed parameters {
vector[K] beta;
...
beta <- beta_scale * (beta_ratios - (1.0 / K));

Proposal B:

parameters {
vector[K-1] beta_prefix;
...
}
transformed parameters {
vector[K] beta;
...
for (k in 1:(K-1))
beta[k] <- beta_prefix[k];
beta[K] <- - sum(beta_prefix);
...

Proposal (A) seems identified to me, but there are
K free parameters (K-1 for the simplex and 1 for the scale),
whereas I'm pretty sure there should only be (K-1), as
in proposal (B).

Proposal (B) feels less natural to me because of the asymmetry in
the beta and the resulting difficulty in finding natural priors.

- Bob

[Bob Carpenter]

I posted a new query that hopefully someone who understands
these matters better than me will take up. (And hopefully
explain how to answer in general.)

I don't see how the pilots.stan model relates to the question.

- Bob

[Michael Betancourt]

> Could someone who understands identification a bit better
> than me help out in figuring out which of the two proposals
> below lead to an identified parameterization of a centered
> vector on the unconstrained scale?

Neither and both? The problem is that there are an infinite number of
ways to map between the unconstrained space and the constrained space.
More formally, the unconstrained space is an affine subspace of the
constrained space, i.e. the unconstrained space decomposes into
an origin and a plane.

Now given any transformation we can always rotate and scale to obtain a new,
equally valid transformation. Moreover there will be combinatorial symmetry
in decided which unconstrained point is interpreted as the “origin”.

> Proposal A:
>
> parameters {
> simplex[K] beta_ratios;
> real<lower=0> beta_scale;
> ...
> transformed parameters {
> vector[K] beta;
> ...
> beta <- beta_scale * (beta_ratios - (1.0 / K));

First thing, the simplex won’t work here because it assumes n _positive_ variables
that sum to 1, not n unconstrained variables.

Second thing, the beta_scale is not necessary. It just indexes many of the possible
transformations you could have (by the scaling noted above) but should not be
interpreted as a parameter. In other words, beta_scale is a constant.

> Proposal B:
>
> parameters {
> vector[K-1] beta_prefix;
> ...
> }
> transformed parameters {
> vector[K] beta;
> ...
> for (k in 1:(K-1))
> beta[k] <- beta_prefix[k];
> beta[K] <- - sum(beta_prefix);

> …

The problem here is not that the transformation is invalid, but rather a uniform distribution
on the unconstrained space is not a uniform distribution on the constrained space, and
vice-versa.

Unfortunately there’s no way to define a truly symmetric transformation — one of the
constrained variables will always have to be the origin which will distinguish it from
the others. You can try to marginalize over all the possible assignments, but then
you have to deal with the Jacobian being a sum of terms (that you would have to
compute before taking the log). Not fun.

Anyways, it’s easy to go from constrained to unconstrained, but trying to go from
unconstrained to constrained will always break the symmetry. As far as I can tell.

[Bob Carpenter]

On 1/26/14, 4:35 PM, Michael Betancourt wrote:
>
>> Could someone who understands identification a bit better
>> than me help out in figuring out which of the two proposals
>> below lead to an identified parameterization of a centered
>> vector on the unconstrained scale?
>
> Neither and both? The problem is that there are an infinite number of
> ways to map between the unconstrained space and the constrained space.
> More formally, the unconstrained space is an affine subspace of the
> constrained space, i.e. the unconstrained space decomposes into
> an origin and a plane.
>
> Now given any transformation we can always rotate and scale to obtain a new,
> equally valid transformation. Moreover there will be combinatorial symmetry
> in decided which unconstrained point is interpreted as the “origin”.

By identified here, I really just meant are there different
values of the parameters that lead to the same transformed parameter
vector?

If you are building something like an IRT model and require one
set of coefficients to have mean 0, it identifies the model as a whole.

>> Proposal A:
>>
>> parameters {
>> simplex[K] beta_ratios;
>> real<lower=0> beta_scale;
>> ...
>> transformed parameters {
>> vector[K] beta;
>> ...
>> beta <- beta_scale * (beta_ratios - (1.0 / K));
>
> First thing, the simplex won’t work here because it assumes n _positive_ variables
> that sum to 1, not n unconstrained variables.

Maybe I wasn't clear on the goal here, which is to
define a vector beta of unconstrained values that
sums to 0.

Although the simplex is constrained to positive
values that sums to 1, (beta_ratios - 1.0 / K) is
a set of values that sums to 0 and whose absolute
values sum to 1. It's only going to be non-negative
in the symmetric case.

For example, if you look at the simplex (0.3, 0.7),
then (0.3, 0.7) - 0.5 = (-0.2, 0.2).

> Second thing, the beta_scale is not necessary. It just indexes many of the possible
> transformations you could have (by the scaling noted above) but should not be
> interpreted as a parameter. In other words, beta_scale is a constant.

The (beta_ratios - 1.0 / K) is constrained to have absolute values
sum to 1. By scaling by beta_scale, the absolute value sums to
beta_scale.

Isn't the result that beta can be any vector that sums to 0?

>> Proposal B:
>>
>> parameters {
>> vector[K-1] beta_prefix;
>> ...
>> }
>> transformed parameters {
>> vector[K] beta;
>> ...
>> for (k in 1:(K-1))
>> beta[k] <- beta_prefix[k];
>> beta[K] <- - sum(beta_prefix);
>> …
>
> The problem here is not that the transformation is invalid, but rather a uniform distribution
> on the unconstrained space is not a uniform distribution on the constrained space, and
> vice-versa.

In this particular case, we're going to put a prior on the parameters, not
on the transformed parameters. That's what I meant by it being a bit odd
because of the asymmetry.

>
> Unfortunately there’s no way to define a truly symmetric transformation — one of the
> constrained variables will always have to be the origin which will distinguish it from
> the others. You can try to marginalize over all the possible assignments, but then
> you have to deal with the Jacobian being a sum of terms (that you would have to
> compute before taking the log). Not fun.
>
> Anyways, it’s easy to go from constrained to unconstrained, but trying to go from
> unconstrained to constrained will always break the symmetry. As far as I can tell.

I think that's right unless we want to apply the Jacobian adjustment.
Then we could put the prior on beta and we'd be back to what I meant
by symmetric in that you could put a prior on each of the beta[k] for
k in 1:K.

- Bob

[Michael Betancourt]

By identified here, I really just meant are there different
values of the parameters that lead to the same transformed parameter
vector?

Yes, I know.  Let me try again.

In general we have a constrained space, C, and an unconstrained space, U.

Now in many cases there is a canonical isomorphism between C and U,

for example with the simplex.

In this case, however, there is no canonical transformation; rather there is

an infinite family of equivalent transformations.  Pick any one and you can map

the constrained space to an unconstrained space where the system is

identified.

With beta_scale you’re trying to parameterize the transformations — but

you don’t need to do that.  You only need to pick out one of them.

Although the simplex is constrained to positive
values that sums to 1, (beta_ratios - 1.0 / K) is
a set of values that sums to 0 and whose absolute
values sum to 1.  It's only going to be non-negative
in the symmetric case.

For example, if you look at the simplex (0.3, 0.7),
then (0.3, 0.7) - 0.5 = (-0.2, 0.2).

Yes, they can be negative but only so much.  The resulting beta have

support on (-1/K, infinity), not (-infinty, infinity).

In this particular case, we're going to put a prior on the parameters, not
on the transformed parameters.  That's what I meant by it being a bit odd
because of the asymmetry.

Right.  But because there’s no canonical transformation between the

constrained and unconstrained space this may not be the best idea.

Unless you’re careful to define an asymmetric prior so that the resulting

prior on the constrained scale is symmetric.

[Andrew Gelman]

I hate these models! We can discuss tomorrow.

A

[Bob Carpenter]

On 1/26/14, 6:11 PM, Michael Betancourt wrote:
>> By identified here, I really just meant are there different
>> values of the parameters that lead to the same transformed parameter
>> vector?
>
> Yes, I know. Let me try again.
>
> In general we have a constrained space, C, and an unconstrained space, U.
> Now in many cases there is a canonical isomorphism between C and U,
> for example with the simplex.

Now I have to figure out what a canonical isomorphism is. The quick
Wikipedia lookup's less than helpful:

http://en.wikipedia.org/wiki/Canonical

Canonical homomorphism, canonical isomorphism: an homomorphism that is uniquely defined by its main property, for
example, the canonical homomorphism of the integers into the rational numbers, the canonical homomorphism of a vector
space onto the quotient by a subspace. More precisely, a canonical morphism or canonical homomorphism is a morphism
whose existence and unicity are provided by a universal property.

This is a bit more useful:

http://www.quora.com/Mathematics/What-is-a-canonical-isomorphism

but it's going to take some thought on my part. Ah, the education continues.

> In this case, however, there is no canonical transformation; rather there is
> an infinite family of equivalent transformations. Pick any one and you can map
> the constrained space to an unconstrained space where the system is
> identified.
>
> With beta_scale you’re trying to parameterize the transformations — but
> you don’t need to do that. You only need to pick out one of them.

I think we'll have to take this up at the next meeting
when there's more bandwidth! I don't even see how
(0,inf) <--> (-inf,inf) has a canonical transform. There
are lots of continuous monotonic functions that map one to
the other. For example, log and log^3 both map
(0,inf) to (-inf,inf).

>> Although the simplex is constrained to positive
>> values that sums to 1, (beta_ratios - 1.0 / K) is
>> a set of values that sums to 0 and whose absolute
>> values sum to 1. It's only going to be non-negative
>> in the symmetric case.
>>
>> For example, if you look at the simplex (0.3, 0.7),
>> then (0.3, 0.7) - 0.5 = (-0.2, 0.2).
>
> Yes, they can be negative but only so much. The resulting beta have
> support on (-1/K, infinity), not (-infinty, infinity).

They can only be -1/K negative, whereas they can be up to
(1 - 1/K) positive. But the final values aren't constrained
at all, at least independently.

If the simplex is (0.3,0.7) and beta_scale = 1e6, then
you get (-2e5,2e5) out.

Going the other way, I'm stuck on the algebra in trying to
solve

Given K-vector beta such that

sum(beta) = 0,

solve for beta_simplex and beta_scale such that

beta = (beta_simplex - 1/K) * beta_scale
sum(beta_simplex) = 1
beta_simplex[k] > 0
beta_scale > 0

I don't think there's going to be a unique solution because there
are too many free variables in beta_simplex and beta_scale compared
to beta.

>> In this particular case, we're going to put a prior on the parameters, not
>> on the transformed parameters. That's what I meant by it being a bit odd
>> because of the asymmetry.
>
> Right. But because there’s no canonical transformation between the
> constrained and unconstrained space this may not be the best idea.
> Unless you’re careful to define an asymmetric prior so that the resulting
> prior on the constrained scale is symmetric.

I'm probably missing something else here, because this seems like the same problem
as for the simplex itself. The simplex has K values, but only (K-1)
underlying values --- the last one gets constrained so that the sum
is 1.

We really do want to discuss parameterizations of multi-logit and multi-probit,
because they bring up the same issues and we should be making recommendations
to the users (and to me!).

- Bob

[Michael Betancourt]

> Now I have to figure out what a canonical isomorphism is.

Sorry I just meant a natural one. For example, the tangent space and cotangent
space on a manifold have a canonical isomorphism because they are dual to
each other. You can define many more, but that’s the only one that arises
naturally from the underlying geometric structure.

> I think we'll have to take this up at the next meeting
> when there's more bandwidth! I don't even see how
> (0,inf) <--> (-inf,inf) has a canonical transform. There
> are lots of continuous monotonic functions that map one to
> the other. For example, log and log^3 both map
> (0,inf) to (-inf,inf).

You’re right — you need to impose additional assumptions to make log natural.
A better example is the logit — logit(ax + b) works just as well as logit(x). In
fact, this provides a nice parallel to this case.

For a multivariate system we should be able to choose any logit( a^{T} X b + c),
where a and b are vectors, X a matrix, and c a scalar (i.e. you can rotate,
translate, and then scale the argument). Now a, b, and c aren’t parameters
but rather constants that determine the mapping and its inverse. Same thing
happens here — you can always rotate and scale the unconstrained space
relative to the constrained space.

> I don't think there's going to be a unique solution because there
> are too many free variables in beta_simplex and beta_scale compared
> to beta.

Yeah, this is where specifying the transformation is really important. Some of the
constants cancel out in the constrained -> unconstrained mapping but matter
for the inverse.

> I'm probably missing something else here, because this seems like the same problem
> as for the simplex itself. The simplex has K values, but only (K-1)
> underlying values --- the last one gets constrained so that the sum
> is 1.

The simplex is different because of the positivity constraint. This limits the potential
transformations — you no longer have translation or rotation invariance — I think you
still have scale invariance but it always cancels out. The lack of rotational invariance is
the biggest feature because it prevents the permutation symmetries we see here.

> We really do want to discuss parameterizations of multi-logit and multi-probit,
> because they bring up the same issues and we should be making recommendations
> to the users (and to me!).

See the discussion above — lots of ways to do it, no natural way. The logit/probit case
is nice because the rotation/scaling/translation have intuitive meaning or can be
absorbed into existing parameters.

[Bob Carpenter]

Thanks! That makes a lot more sense now. I'm sure
I'm missing a lot of the subtleties, but to paraphrase
Daniel, baby steps are good if they're in the right direction.

- Bob

[Joseph Burling]

Trying to diagnose this mixing issue and wanted to hear yours thoughts about the trace for a few parameters. Seems like an identifiability problem to me relating back to the centering issue but not sure.

bad mixing example

some mixing but weird spikes

Seems okay to me

[Bob Carpenter]

The first example does seem like identifiabilit and the third
seems OK, but if you don't know the correct answer, it's hard
to say.

The "weird spikes" may be a result of one of the chains (the yellow
one here) finding a better step size and thus exploring more of
the posterior. This may also be due to identifiability.

The only way to know for sure is to simulate fake data where you
know the answer and fit that. Even then, unless there's an
analytic posterior you won't know what the posterior variance
should be.

- Bob

On 1/29/14, 3:54 PM, Joseph Burling wrote:
> Trying to diagnose this mixing issue and wanted to hear yours thoughts about the trace for a few parameters. Seems like
> an identifiability problem to me relating back to the centering issue but not sure.
>
>
> bad mixing example
>

> <https://lh3.googleusercontent.com/-9tUvFOhca3w/UukVjSoIcpI/AAAAAAAAHOU/0dCFFnfSQd0/s1600/Screen+Shot+2014-01-29+at+8.47.36+AM.png>

>
>
> some mixing but weird spikes
>

> <https://lh6.googleusercontent.com/-5-26e7VXU0M/UukVzcu7fuI/AAAAAAAAHOc/wNfuO-H0DF0/s1600/Screen+Shot+2014-01-29+at+8.47.02+AM.png>
>
>
> Seems okay to me
>
> <https://lh5.googleusercontent.com/-utuvJNdX0Ss/UukV6wR4dWI/AAAAAAAAHOk/zoQN3TrfVOg/s1600/Screen+Shot+2014-01-29+at+8.43.57+AM.png>

>
>
>
>
>
>
>
> On Monday, January 20, 2014 9:17:15 AM UTC-6, Bob Carpenter wrote:
>
>
>
> On 1/20/14, 3:59 PM, Joseph Burling wrote:
> >
> > First of all, thank you for the quick and thorough response.
> >
> >
> > I don't understand the model's parameterization as is, but
> > maybe someone else will have a clue. In particular, I'm not
> > sure what's going on with all the centering you're doing,
> > such as
> >
> > mu_d[k] <- row(beta_dyad_raw, k) - mean(beta_dyad_raw[k]);
> >
> > mu_o[k] <- row(beta_obj_raw, k) - mean(beta_obj_raw[k]);
> >
> > epsilon_i <- epsilon_raw_adj - mean(epsilon_raw_adj);
> >
> >
> >
> > I did the centering this way because I was trying to get mean population beta parameters (Beta_mu) and the above
> are the
> > group level betas centered at zero. I know I had a prior of zero to begin with but often one groups' parameters
> might
> > get stuck far from zero and switch with the population beta parameters.
>
> The real problem with this kind of centering is that the model
> winds up being unidentified on the unconstrained scale (the centered
> parameters are identified, but not the parameters being sampled, which are just
> the transforms of the variables defined in the parameters block).
>
> Andrew recommended the approach you took for BUGS in his regression book
> with Jennifer, but it's problematic for Stan because the adaptation
> and sampling happens on the unconstrained scale.
>
>