Hierarchical (Conditional) Multinomial Logit Example

[Sanjog Misra]

First, I would like to thank the Stan team for this most generous gift! I can only imagine the time and effort that has gone into this project. So thanks.

I have been playing around with Stan and the results are impressive. I started with a Hierarchical Linear Regression (as a response to Ian's earlier question on imbalanced panels) and the results were great.

(Blog pose here: www.sanjogmisra.com/blog/stan/)

Since then I have been playing around with a "real" example of a Hierarchical Multinomial Logit (Mixed Logit) as applied to a brand choice context (I'm in a business school!) 

I have attached files for those who are interested.  The code works and the chains generated converge very fast (realtive to BUGS/JAGS). Having said that, I have run into some issues that I was hoping others (the Stan team) could help with:

a) The compilation and sampling time is somewhat long (my comparison code in JAGS compiles and samples much faster but doesn't mix as well).  

b) At some points i get a "Non-Fatal" warning that the inverse Wishart draw is bad. 

Warning (non-fatal rejection): Error in function stan::prob::inv_wishart_log<i>(i): W is not positive definite (0)

c) While the chain converges fast in most cases, once in a while the sampling slows to a crawl and the resultant chain shows poor mixing (see attached pdf where chain #1 has poor mixing). Not sure why this happens. Sometimes the acceptance rates drop quite a bit even after the chain seems to have gotten to the "correct" parameter range. 

Any and all other thoughts are welcome. I am in the process of benchmarking Stan with JAGS, BUGS and available code in R but wanted to make sure I get the most efficient version before I do so. My guess is that I could vectorize my code better but cant seem to figure out how. 

Some Notes: 

I am using the -O3 flag.

I have tried using the multinomial distribution (rather than lp__) that doesn't help speed things up.

I have tried the multi_normal_cholesky (doesn't help speed)

tried and found that softmax function is not available in current version (remove from documentation?)

Looking forward to comments, ideas and suggestions.

Sanjog 

[Ben Goodrich]

On Saturday, September 22, 2012 6:47:16 PM UTC-4, Sanjog Misra wrote:

Some Notes: 

I am using the -O3 flag.

I have tried using the multinomial distribution (rather than lp__) that doesn't help speed things up.

I have tried the multi_normal_cholesky (doesn't help speed)

tried and found that softmax function is not available in current version (remove from documentation?)

Looking forward to comments, ideas and suggestions.

 

etamu ~ multi_normal(mu_m_prior, mu_var_prior);

For any multivariate model where the covariance matrix is known (data), many of the density functions currently in Stan are inefficient because they effectively do a matrix inverse every time the log posterior function is evaluated. In this case, it is faster to calculate the inverse of mu_var_prior once (to get the prior precision matrix) and better still the transpose of Cholesky factor (t_L) of the prior precision matrix. Then, you can implement the log multivariate normal yourself with

v <- t_L * (etamu - mu_m_prior);
lp__ <- lp__ - 0.5 * dot_product(v,v);

Note that the determinant part of the multivariate normal is a constant so you don't need to worry about that (although it is different below). Also, in your case mu_m_prior is a zero vector, so you can just drop that.

In this line,

eta[h] ~ multi_normal(etamu,etavar);


etavar is an unknown parameter, so using multi_normal() is okay. However, multi_normal_cholesky() is faster if you parameterize the model in terms of the elements of the Cholesky factor of the covariance matrix (etavar). But, doing so presents a complication if you want to put an inverse Wishart prior on the covariance matrix because you have to handle the Jacobian yourself. Something like

parameters {
  real subdiagonals[(H * (H -1)];
  real<lower=0> diagonals[H];
  // other parameters
}
model {
  matrix[H,H] L_etavar;
  int counter;
  counter <- 1;
  for(i in 1:H) {
    for(j in i:(H - 1)) {
      L_etavar[i,j] <- subdiagonals[counter];
      counter <- counter + 1;
    }
    L_etavar[i,i] <- diagonals[i];
    for(j in (i+1):H) {
      L_etavar[i,j] <- 0.0;
    }
  }

  for(h in 1:H) {
    eta[h] ~ multi_normal_cholesky(etamu,L_etavar);
  }

  L_etavar * L_etavar' ~ inv_wishart((D+K-1), invR);
  // change-of-variables see "Absolute Jacobian Determinant of the Covariance Matrix Inverse Transform" in Stan manual
  for(h in 1:H) {
    lp__ <- lp__ + (H - h + 2) * log(diagonals[h]);
  }

  // rest of model ...
}

At the moment (this could change in the next version of Stan), multi_normal() and multi_normal_cholesky() are not vectorized, which is why you have to loop. This is regrettable because L_etavar gets inverted H times when it only needs to be inverted once. For the time being, you could get around this by parameterizing the model in terms of the Cholesky factor of the precision matrix rather than the covariance matrix. This resembles what we did above, except now you have to account for the determinant part of the (log) multivariate normal

parameters {
  real subdiagonals[(H * (H -1)];
  real<lower=0> diagonals[H];
  // other parameters
}
model {
  vector[H] v;
  matrix[H,H] L_etaprec;
  int counter;
  counter <- 1;
  for(i in 1:H) {
    for(j in i:(H - 1)) {
      L_etaprec[i,j] <- subdiagonals[counter];
      counter <- counter + 1;
    }
    L_etaprec[i,i] <- diagonals[i];
    for(j in (i+1):H) {
      L_etaprec[i,j] <- 0.0;
    }
  }

  for(h in 1:H) {
    v <- L_etaprec' * (eta - etamu);
    lp__ <- lp__ - 0.5 * dot_product(v,v);
  }
  lp__ <- lp__ + H * sum(log(diagonals)); // log-determinant part

  L_etaprec * L_etaprec' ~ wishart((D+K-1), R);
  // change-of-variables see "Absolute Jacobian Determinant of the Covariance Matrix Inverse Transform" in Stan manual
  for(h in 1:H) {
    lp__ <- lp__ + (H - h + 2) * log(diagonals[h]);
  }

  // rest of model ...
}

Also, note that your prior over L_etaprec * L_etaprec' = etaprec is now Wishart, as opposed to having an inverse Wishart prior over L_etavar * L_etavar' = etavar.

Hopefully, that makes sense and hopefully, that will all be less cumbersome soon.

Ben

P.S. Double-check my Stan code for typos / thinkos.



[terrance savitsky]

Why does one have to perform the jacobian increment, manually, when placing  LL' ~ inv_wishart() prior on for L, lower triangular?  Does Stan recognize L as a cholesky factor and therefore not perform the jacobian adjustment to the unrestricted sampling form?  this seem illogical.   how would one have any intuition that they would need to handle the jacobian, themselves?  in particular, for cov_matrix[D] fred;  with fred = LL'; I would not have to perform the adjustment for fred ~ inv_wishart();

--
Thank you, Terrance Savitsky

[Ben Goodrich]

On Saturday, September 22, 2012 9:24:04 PM UTC-4, tds151 wrote:

Why does one have to perform the jacobian increment, manually, when placing  LL' ~ inv_wishart() prior on for L, lower triangular?  Does Stan recognize L as a cholesky factor and therefore not perform the jacobian adjustment to the unrestricted sampling form?  this seem illogical.   how would one have any intuition that they would need to handle the jacobian, themselves?  in particular, for cov_matrix[D] fred;  with fred = LL'; I would not have to perform the adjustment for fred ~ inv_wishart();

This is certainly something that is unique to Stan. I believe the rule is:

1) If the LHS of a ~ is declared in the parameters block, then the user does NOT need to adjust lp__ with a log absolute Jacobian determinant because Stan will do it automatically if necessary. Similarly, if updating lp__ manually instead of using the ~ syntax, if the variable(s) that are contributing that that lp__ update are declared in the parameters block, then the user does NOT need to further adjust lp__ with a log absolute Jacobian determinant because Stan will do it automatically if necessary

2) If the LHS of a ~ is declared in the transformed parameters block or is declared as a local variable in the model block, then the user DOES need to adjust lp__ with a log absolute Jacobian determinant. Similarly, if updating lp__ manually instead of using the ~ syntax, if the variable(s) that are contributing that that lp__ update are declared in the transformed parameters block or the model block, then the user DOES need to further adjust lp__ with a log absolute Jacobian determinant.

When the user does this in the .stan file:

parameters{
  cov_matrix[K] Sigma;
}

The C++ code utilizes a Cholesky factor (L) --- even though L is not mentioned in the .stan file --- and comes out looking like:

Sigma = L * L.transpose();
lp__ += a log absolute Jacobian determinant

In other words, the user is utilizing the Cholesky factor without necessarily knowing it.

However, the C++ code is always utilizing unbounded primitive parameters that correspond to the parameters declared in the .stan file. The only exception is

parameters {
  real mu;
}

Since mu is unrestricted, there is no need to complicate things with a primitive parameter. In contrast,

parameter {
  real<lower=0> sigma;
}

is positive so the C++ will have something like

sigma = exp(sigma_primitive);
lp__ += a log absolute partial derivative

In other words, a Cholesky factor is the primitive parameter underlying a covariance matrix, just like sigma_primitive is the primitive parameter underlying sigma.

See the chapter of the manual on "Transformations of Variables" for more details.

Ben

[Ben Goodrich]

On Saturday, September 22, 2012 9:59:02 PM UTC-4, Ben Goodrich wrote:

2) If the LHS of a ~ is declared in the transformed parameters block or is declared as a local variable in the model block, then the user DOES need to adjust lp__ with a log absolute Jacobian determinant. Similarly, if updating lp__ manually instead of using the ~ syntax, if the variable(s) that are contributing that that lp__ update are declared in the transformed parameters block or the model block, then the user DOES need to further adjust lp__ with a log absolute Jacobian determinant.

P.S. Bob mentioned on stan-dev that we should introduce foo_transform(theta) functions that would adjust lp__ with a log absolute Jacobian determinant for the inverse foo transformation (to theta). Thus, users would need to remember to call foo_transform(theta) in the model block, but would not need to work out what the log absolute Jacobian determinant is.

Hopefully, that will happen soon, but probably not for Stan 1.0.2 .

Ben

[terrance savitsky]

Thanks, Ben.   I now understand the distinction *when* the user needs to include the (log absolute determinant) jacobian adjustment.  that said, i don't understand why restricted parameters declared in the transformed parameters block if one specifies a distribution in the model block requires the adjustment.  If the restricted parameters in the transformed parameters block are a function of parameters in the parameters block (in fashion that obeys previously declared restrictions) and I place a distribution on the transformed parameters in the model block, then the sampling focus shifts to the transformed parameters.   I know one can't declare local parameters as restricted in the models block, but when may do so in the parameters block, right?  i ask these questions b/c if I understand the outlines of how the stan dev team thinks then i may write smarter .stan scripts.

[Ben Goodrich]

On Saturday, September 22, 2012 10:26:43 PM UTC-4, tds151 wrote:

Thanks, Ben.   I now understand the distinction *when* the user needs to include the (log absolute determinant) jacobian adjustment.  that said, i don't understand why restricted parameters declared in the transformed parameters block if one specifies a distribution in the model block requires the adjustment.  If the restricted parameters in the transformed parameters block are a function of parameters in the parameters block (in fashion that obeys previously declared restrictions) and I place a distribution on the transformed parameters in the model block, then the sampling focus shifts to the transformed parameters.   I know one can't declare local parameters as restricted in the models block, but when may do so in the parameters block, right?  i ask these questions b/c if I understand the outlines of how the stan dev team thinks then i may write smarter .stan scripts.

It isn't really true that if you "place a distribution on the transformed parameters in the model block, then the sampling focus shifts to the transformed parameters". In your mind, you might be focusing on the transformed parameters, but from Stan's perspective the sampling is always done on the unbounded, unrestricted, primitive parameters (which the user never even sees). If Stan sampled directly from the transformed parameter space, the sampler (often) might get stuck in a region where the posterior evaluates to zero or is not defined. So, Stan samples from a parameter space where the posterior is defined everywhere and is non-zero everywhere. [In the future, Stan might have an option to bounce around in a transformed parameter subspace, but I doubt that will happen real soon]

Maybe you are thinking of a different example, but let me try to clarify using the AR1 example again. If your .stan file were

data {
  int<lower=1> K;
}
parameters{
  real<lower=0,upper=1> phi;
  real<lower=0> sigma_squared;
}
transformed parameters {
  cov_matrix[K] Sigma;
  // loops to fill in Sigma as a function of phi and sigma_squared
}
model {
  // informative prior on phi
  phi ~ beta(1,2);

  // informative prior on sigma_squared
  sigma_squared ~ exp(1);

  // should NOT have a prior on Sigma
  // should NOT have a Jacobian adjustment

  // some likelihood of the data as a function of Sigma
}

Alternatively --- although this makes little substantive sense to me --- you could do

data {
  int<lower=1> K;
  int<lower=K+1> df;
  matrix[K,K] scale;
}
parameters {
  real<lower=0,upper=1> phi;
  real<lower=0> sigma_squared;
}
transformed parameters {
  cov_matrix[K] Sigma;
  // loops to fill in Sigma as a function of phi and sigma_squared
}
model {
  // prior on Sigma with Jacobian adjustment
  Sigma ~ wishart(df,scale);  
  lp__ <- lp__ + SomeJacobianAdjustment;

  // should NOT have a prior on phi or sigma_squared

  // some likelihood of the data as a function of Sigma
}

In other words, a Wishart prior on the transformed parameter Sigma implies some prior on phi and sigma_squared (technically, it implies some prior on the primitive parameters underlying phi and sigma_squared but Stan automatically adjusts for the change-of-variables to phi and sigma_squared). The reason why I don't like this on substantive grounds is that it is almost completely opaque as to what a Wishart prior on the transformed parameter Sigma implies about phi and sigma_squared. But it is legal in Stan.

Conversely, I do like the first version which has informative priors on phi and sigma_squared (and no prior on Sigma or Jacobian adjustment) because I can easily understand what these two priors entail, even if I might disagree with them. In general, I would say that it is better for a model to put its priors on the parameters that compose a structured covariance matrix than to put a prior on the composed covariance matrix.

I would actually go further and say that --- if no other structure is attractive --- compose a covariance matrix as the product of a diagonal matrix of standard deviations, a correlation matrix, and a diagonal matrix of standard deviations. There are lots of choices of priors for standard deviations and Stan has a nice prior distribution for a correlation matrix called lkj_corr that takes a positive shape parameter, eta. If eta = 1, the prior is jointly uniform over all correlation matrices of order K, and to the extent that eta > 1, the prior has a sharper mode at the identity matrix of order K.

You see a lot of Wishart and inverse Wishart priors in the literature, but that is often because such priors are necessary for a Gibbs sampler. In Hamiltonian Monte Carlo, conjugacy is not relevant, so Stan users might as well use prior distributions that they believe or at least understand. With Wishart and inverse Wishart priors, the expectation is known but that's about it.

Ben

[Sanjog Misra]

Ben,

I took your suggestion and recoded (file attached) using the Cholesky decomposition approach. I had to make some corrections regarding indices etc. bt I think the model is correct. One problem was that the log functions doesnt work on a vector of real<lower=0> type so I had to loop that. Otherwise the model seems to compile fine but gives a segfault when sampling...

Thoughts?

=========

g++ -arch x86_64 -I/Library/Frameworks/R.framework/Resources/include -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG   -I"/Users/sanjogmisra/Library/R/2.15/library/rstan/include//stansrc"  -I"/Users/sanjogmisra/Library/R/2.15/library/RcppEigen/include/"  -I"/Users/sanjogmisra/Library/R/2.15/library/RcppEigen/include//unsupported"  -I"/Users/sanjogmisra/Library/R/2.15/library/rstan/include//stanlib/boost_1.51.0" -I"/Users/sanjogmisra/Library/R/2.15/library/rstan/include" -I/usr/local/include -I"/Users/sanjogmisra/Library/R/2.15/library/Rcpp/include"   -fPIC  -g -O3  -c file910317019b9b.cpp -o file910317019b9b.o

g++ -arch x86_64 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/usr/local/lib -o file910317019b9b.so file910317019b9b.o /Users/sanjogmisra/Library/R/2.15/library/Rcpp/lib/x86_64/libRcpp.a /Users/sanjogmisra/Library/R/2.15/library/rstan/libstan/x86_64/libstan.a -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation

SAMPLING FOR MODEL 'hmnlcode' NOW (CHAIN 1).

 *** caught segfault ***

address 0x0, cause 'unknown'

Traceback:

 1: .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x11bff84d0>,     dll = list(name = "Rcpp", path = "/Users/sanjogmisra/Library/R/2.15/library/Rcpp/libs/x86_64/Rcpp.so",         dynamicLookup = TRUE, handle = <pointer: 0x115c15900>,         info = <pointer: 0x10040adf0>), numParameters = -1L),     <pointer: 0x1181408a0>, <pointer: 0x118178890>, .pointer,     ...)

 2: sampler$call_sampler(args_list[[i]])

 3: .local(object, ...)

 4: sampling(sm, data, pars, chains, iter, warmup, thin, seed, init,     sample_file = sample_file, verbose = verbose, check_data = FALSE,     ...)

 5: sampling(sm, data, pars, chains, iter, warmup, thin, seed, init,     sample_file = sample_file, verbose = verbose, check_data = FALSE,     ...)

 6: stan(model_code = hmnlcode, data = dat, iter = 2000, chains = 2,     verbose = TRUE)

[Jiqiang Guo]

I thought it was a problem in RStan, so I have singled out the stan code and data.  It turns out it is from somewhere else. 

For the time being, I just attached the stan and rdump file.  I guess somebody else better than me at finding out the very problem would take a book before I have time. 

Jiqiang

[Ben Goodrich]

@Bob,

The problem seems to be near line 293 of the attached .cpp file

            for (int h = 1; h <= H; ++h) {
                assign(vk, multiply(transpose(L_etaprec),subtract(get_base1(eta,h,"eta",1),etamu)));
                assign(lp__, (lp__ - (0.5 * dot_product(vk,vk))));
            }

For some reason, vk is not getting assigned right, as can be seen in the attached gdb log. Its first two elements are on the order of 10^-317 and accessing its third element causes an assert. But eta, etamu, and L_etaprec seem to have valid data, so I don't know what is going wrong.

Ben

[Jeffrey Arnold]

Ben, 

Could you add your discussion of the rules of when the user needs to manually adjust lp__ to the manual? Unless I am missing something, the only discussion is in Ch 16.2 and it only talks about the case in which the LHS of ~ is an expression or lp__ is adjusted directly. 

And to make sure that I am understanding this correctly, for the following model

parameters {
  real<lower=0> sigma;
  ...
}
transformed parameters {
  real<lower=0> tau;
  tau <- sqrt(sigma);
}
model {
  ...
  tau ~ gamma(0.1, 0.1);
  ...
}

requires an extra statement in the model block to adjust the log-posterior

lp__ <-  lp__ + log(0.5 * pow(sigma, -0.5));

That wasn't a real model, and in practice I'd put the prior on sigma and not tau (and not use the evil gamma), but I just want to make sure I understand where adjustments to the log posterior is needed.

Also this may be a terrible idea, but I'll throw it out there anyways. Since accumulation seems to be used often in Stan, any chance 

of adding the equivalent of =+ and =-?  I guess that would be <-+ and <-- (the later being definitely suboptimal syntax).

---

Jeffrey Arnold

Department of Political Science

University of Rochester

http://jrnold.me

jeffrey...@gmail.com

jeffrey...@rochester.edu

[Ben Goodrich]

On Sunday, September 23, 2012 8:01:01 PM UTC-4, Jeffrey Arnold wrote:

Ben, 

Could you add your discussion of the rules of when the user needs to manually adjust lp__ to the manual? Unless I am missing something, the only discussion is in Ch 16.2 and it only talks about the case in which the LHS of ~ is an expression or lp__ is adjusted directly. 

And to make sure that I am understanding this correctly, for the following model

parameters {
  real<lower=0> sigma;
  ...
}
transformed parameters {
  real<lower=0> tau;
  tau <- sqrt(sigma);
}
model {
  ...
  tau ~ gamma(0.1, 0.1);
  ...
}

requires an extra statement in the model block to adjust the log-posterior

lp__ <-  lp__ + log(0.5 * pow(sigma, -0.5));

You are correct that this model requires the manual adjustment for the change-of-variables. But I think you have it reversed: What is needed is the (log absolute) derivative of the transformation from tau to sigma. Since sigma = tau^2, d sigma / d tau = 2.0 * tau and

lp__ <- lp__ + log(2.0 * tau); // no abs() necessary because already positive

Chapter 26 of 1.0.1 reference manual covers this in the most detail. If this were a real model that utilized both sigma and tau for some reason, I would declare tau as the parameter and sigma as the transformed parameter (or as a local parameter within the model block) if you wanted to place an informative prior on gamma.

Also this may be a terrible idea, but I'll throw it out there anyways. Since accumulation seems to be used often in Stan, any chance 

of adding the equivalent of =+ and =-?  I guess that would be <-+ and <-- (the later being definitely suboptimal syntax).

Maybe. I think we are trying to avoid adding syntax unless there is an appreciable benefit. In some languages like C++,

x += 42;

is actually faster than

x = x + 42;

because temporaries are avoided, but that is not really a relevant consideration for the Stan parsing language.

I also wouldn't necessarily agree that accumulation needs to be used "often". It is just that we have been getting a lot of questions from people with big multivariate models, and Stan's implementations of the relevant densities work for everything but are optimal for nothing. So, for now, the fastest way to use multivariate normal / t or inverse / plain Wishart distributions is with hand-optimized accumulation, but it won't be that way forever.

Ben

[Jeffrey Arnold]

Thanks for clearing that up; although now I'm feeling pretty stupid for calculating the wrong derivative.   Chapter 26 covers the details of the transformations, but the manual doesn't explicitly mention that parameters defined in the  transformed parameters block on the LHS of ~ need this adjustment to lp__.  However, after writing this, I saw that the example basic_distributions/wishart2x2.stan is an example where the log Jacobian transform is added to lp__.

You're correct, "often" was not the correct adjective. More like "unusual", in which case you wouldn't want syntatic sugar which would make it easier to not think about what you're doing. 

Jeff

[Ben Goodrich]

On Sunday, September 23, 2012 11:02:08 PM UTC-4, Jeffrey Arnold wrote:

Thanks for clearing that up; although now I'm feeling pretty stupid for calculating the wrong derivative.   Chapter 26 covers the details of the transformations, but the manual doesn't explicitly mention that parameters defined in the  transformed parameters block on the LHS of ~ need this adjustment to lp__. However, after writing this, I saw that the example basic_distributions/wishart2x2.stan is an example where the log Jacobian transform is added to lp__.

I think we are going to have to beef up the manual in this context. I think that

https://code.google.com/p/stan/source/browse/src/models/basic_distributions/wishart2x2.stan

is supposed to be mathematically equivalent to

https://code.google.com/p/stan/source/browse/src/models/basic_distributions/wishart.stan

but the former requires the user to adjust for the Jacobian while Stan handles the later behind the scenes because the cov_matrix is declared in the parameters block.

Ben
 

[Ben Goodrich]

On Sunday, September 23, 2012 11:02:08 PM UTC-4, Jeffrey Arnold wrote:

Thanks for clearing that up; although now I'm feeling pretty stupid for calculating the wrong derivative.  

Argh. You were right the first time; I was wrong. Sorry. Everyone vote down my response on Google Groups.

Do-over. You originally said:

parameters {
  real<lower=0> sigma;
  ...
}
transformed parameters {
  real<lower=0> tau;
  tau <- sqrt(sigma);
}
model {
  ...
  tau ~ gamma(0.1, 0.1);
  ...
}

requires an extra statement in the model block to adjust the log-posterior

lp__ <-  lp__ + log(0.5 * pow(sigma, -0.5));

If you are going to use informative priors, Stan makes you apply them to things in the parameters block only. However, since

p(sigma | ...) \propto abs(det( J[sigma -> gamma] )) * p(gamma | ...)

a specification for p(gamma | ...) commits you to a specification for p(sigma | ...) but you have to adjust the log posterior yourself. In this case, the Jacobian is 1x1, so it is just the partial derivative of gamma with respect to sigma (which I totally messed up in my first attempt at a response). And thus

lp__ <-  lp__ + log(0.5 * pow(sigma, -0.5));

is the correct adjustment.
 

Chapter 26 of 1.0.1 reference manual covers this in the most detail. If this were a real model that utilized both sigma and tau for some reason, I would declare tau as the parameter and sigma as the transformed parameter (or as a local parameter within the model block) if you wanted to place an informative prior on gamma.

That is still intended to be a factual statement. Also, as Jeff points out, the simplest non-trivial example of a multivariate change-of-variables is

https://code.google.com/p/stan/source/browse/src/models/basic_distributions/wishart2x2.stan

Ben

[Ben Goodrich]

On Monday, September 24, 2012 12:42:31 AM UTC-4, Ben Goodrich wrote:

On Sunday, September 23, 2012 11:02:08 PM UTC-4, Jeffrey Arnold wrote:

Thanks for clearing that up; although now I'm feeling pretty stupid for calculating the wrong derivative.  

Argh. You were right the first time; I was wrong. Sorry. Everyone vote down my response on Google Groups.

Do-over. You originally said:

parameters {
  real<lower=0> sigma;
  ...
}
transformed parameters {
  real<lower=0> tau;
  tau <- sqrt(sigma);
}
model {
  ...
  tau ~ gamma(0.1, 0.1);
  ...
}

requires an extra statement in the model block to adjust the log-posterior

lp__ <-  lp__ + log(0.5 * pow(sigma, -0.5));

If you are going to use informative priors, Stan makes you apply them to things in the parameters block only. However, since

p(sigma | ...) \propto abs(det( J[sigma -> gamma] )) * p(gamma | ...)

a specification for p(gamma | ...) commits you to a specification for p(sigma | ...) but you have to adjust the log posterior yourself. In this case, the Jacobian is 1x1, so it is just the partial derivative of gamma with respect to sigma (which I totally messed up in my first attempt at a response).

Damn. Yesterday was not a good day. The transformed parameter in question is tau, not gamma, although it was given a gamma distribution in the example. So, this should have said:

If you are going to use informative priors, Stan makes you apply them to things in the parameters block only. However, since

p(sigma | ...) \propto abs(det( J[sigma -> tau] )) * p(tau | ...)

a specification for p(tau | ...) commits you to a specification for p(sigma | ...) but you have to adjust the log posterior yourself. In this case, the Jacobian is 1x1, so it is just the partial derivative of tau with respect to sigma (which I totally messed up twice, so far).

Ben

[Bob Carpenter]

On 9/23/12 7:06 PM, Ben Goodrich wrote:> @Bob,

>
> The problem seems to be near line 293 of the attached .cpp file
>
> |
> for (int h = 1; h <= H; ++h) {
> assign(vk, multiply(transpose(L_etaprec),subtract(get_base1(eta,h,"eta",1),etamu)));
> assign(lp__, (lp__ - (0.5 * dot_product(vk,vk))));
> }
>
> |
>
> For some reason, vk is not getting assigned right, as can be seen in the attached gdb
> log. Its first two elements are on
> the order of 10^-317 and accessing its third element causes an assert. But eta, etamu,
> and L_etaprec seem to have valid
> data, so I don't know what is going wrong.

...snip...

Thanks for the dumps. I don't think this is the problem as
I recreated a standalone unit test which seems to work just
fine (see below and test/agrad/matrix_test.cpp in the repository).

I don't know what to do in the face of large, complex
user models like this one. I just don't have the time to
spend debugging.

One thing users can do to help is to try to isolate
where the problem is and create a small model that causes
the problem. I know this can be hard, but without it,
it's very hard for me to dive into these big models.

Before getting back to the case at hand,
I want to recommend that users (including me when I'm
using Stan) should build models up stepwise from simple to
complex. That lets you start with something that works,
then gradually add new effects, new hierarchical levels,
new priors, etc. A diff between a model that worked and
the one step more complex model that doesn't is invaluable
for debugging.

The print statements coming in Stan 1.0.2 should help in
these models with complex assignments, too.

The seg fault could come about if one of the variables in
the matrix was not initialized (then the .vi_->val will
cause a segfault) or we somehow to tried to index beyond the
size of either v1 or v2 in v1[i] and v2[i]:

Program received signal SIGSEGV, Segmentation fault.
0x000000000040f565 in stan::agrad::(anonymous
namespace)::dot_product_vv_vari::var_dot<Eigen::Block<Eigen::Matrix<stan::agrad::var, -1, -1, 0, -1, -1> const, 1, -1,
false, true>, Eigen::Block<Eigen::Matrix<stan::agrad::var, -1, 1, 0, -1, 1> const, -1, 1, true, true> > (v1=..., v2=...)
at src/stan/agrad/matrix.hpp:672
672 result += v1[i].vi_->val_ * v2[i].vi_->val_;
(gdb) bt

I created an independent unit test in test/agrad/matrix_test.cpp,
in which I recreated the set of ops that fail in the user model,
but made sure everything was initialized properly.

So I'm guessing something's not getting defined somewhere in the
user model. As soon as prints are in for 1.0.2, this kind of
thing should be easier to debug.

Here's the new test case:

TEST(AgradMatrix, UserCase1) {
using std::vector;
using stan::math::multiply;
using stan::math::transpose;
using stan::math::subtract;
using stan::math::get_base1;
using stan::agrad::assign;
using stan::math::dot_product;

size_t H = 2;
size_t DpKm1 = 4;

vector_v vk(DpKm1);
for (size_t k = 0; k < DpKm1; ++k)
vk[k] = (k + 1) * (k + 2);

matrix_v L_etaprec(DpKm1,DpKm1);
for (size_t m = 0; m < DpKm1; ++m)
for (size_t n = 0; n < DpKm1; ++n)
L_etaprec(m,n) = (m + 1) * (n + 1);

vector_d etamu(DpKm1);
for (size_t k = 0; k < DpKm1; ++k)
etamu[k] = 10 + (k * k);

vector<vector_d> eta(H,vector_d(DpKm1));
for (size_t h = 0; h < H; ++h)
for (size_t k = 0; k < DpKm1; ++k)
eta[h][k] = (h + 1) * (k + 10);

var lp__ = 0.0;

for (size_t h = 1; h <= H; ++h) {

assign(vk, multiply(transpose(L_etaprec),
subtract(get_base1(eta,h,"eta",1),
etamu)));
assign(lp__, (lp__ - (0.5 * dot_product(vk,vk))));
}

EXPECT_TRUE(lp__.val() != 0);
}

- Bob

[Bob Carpenter]

Bunch of comments inline.

- Bob

On 9/23/12 8:58 PM, Ben Goodrich wrote:
> On Sunday, September 23, 2012 8:01:01 PM UTC-4, Jeffrey Arnold wrote:
>
> Ben,
>
> Could you add your discussion of the rules of when the user needs to manually adjust lp__ to the manual? Unless I am
> missing something, the only discussion is in Ch 16.2 and it only talks about the case in which the LHS of ~ is an
> expression or lp__ is adjusted directly.

We'll add a change of variables chapter to the
user's guide.

I want to try to make clear when you can transform
without a Jacobian and when's one necessary.

Which specific issue were you concerned about adding?
Just Ben's two-point post?

> And to make sure that I am understanding this correctly, for the following model
>
> parameters {
> real<lower=0> sigma;
> ...
> }
> transformed parameters {
> real<lower=0> tau;
> tau <- sqrt(sigma);
> }
> model {
> ...
> tau ~ gamma(0.1, 0.1);
> ...
> }
>
> requires an extra statement in the model block to adjust the log-posterior
>
> lp__ <- lp__ + log(0.5 * pow(sigma, -0.5));

I'd suggest putting it in the transformed parameters block.

And that is the derivative you need. I got SO confused
with the directionality of these in our first implementation
of Stan, that I can feel everyone's pain dealing with this
very acutely. The transform here (at least in the notation
of Gelman et al.'s BDA and our manual) is from tau to
sigma. The inverse transfrom is from sigma to tau,
and it's the log Jacobian of the inverse that needs to be added.

We'll add some simple univariate, so you should be able to do
this:

transformed Parameters {
real<lower=0> tau;
tau <- sqrt_transform(sigma);
}

or go even further and eliminate the intermediate variable and
just use this:

model {
...
sqrt_transform(sigma) ~ gamma(0.1,0.1);


The sqrt_transform() function bring in the Jacobian
itself behind the scenes. The discussion we've been having is
how to compose these in multivariate settings because we can't easily compute
and then auto-dif Jacobians (we can do one or the other, but
not second-order derivatives).

> Also this may be a terrible idea, but I'll throw it out there anyways. Since accumulation seems to be used often in
> Stan, any chance
> of adding the equivalent of =+ and =-? I guess that would be <-+ and <-- (the later being definitely suboptimal
> syntax).

Yes, that's on the to-do list. The usual syntax
in C, Java, etc. for compound assignment has the operator
first. The statement

a += b;

has the same effect on a as:

a = a + b;

Syntactically and semantically it's no issue to add,
but I'm uncertain what to call it given how ugly
"+<-" looks. But that's what it'll have to be.

An alternative would be:

lp_increment(x);

to stand in for:

lp <- lp + x;

Thoughts?

> Maybe. I think we are trying to avoid adding syntax unless there is an appreciable benefit.

I think there's a benefit to this just as in C++.
It's cleaner and clearer. I just haven't added it
because our to-do list is already a mile long.

> In some languages like C++,
>
> x += 42;
>
> is actually faster than
>
> x = x + 42;
>
> because temporaries are avoided, but that is not really a relevant consideration for the Stan parsing language.

It's relevant insofar as we translate

lp__ <- lp__ + x;

into C++ as:

lp__ = lp__ + x;

I'm guessing the compilers are good enough to figure
out this can be optimized as well as:

lp += x;

> I also wouldn't necessarily agree that accumulation needs to be used "often". It is just that we have been getting a lot
> of questions from people with big multivariate models, and Stan's implementations of the relevant densities work for
> everything but are optimal for nothing. So, for now, the fastest way to use multivariate normal / t or inverse / plain
> Wishart distributions is with hand-optimized accumulation, but it won't be that way forever.

I want to have this available for sums in loops, too.

- Bob

[Jeffrey Arnold]

comments inline

 

On 9/23/12 8:58 PM, Ben Goodrich wrote:
> On Sunday, September 23, 2012 8:01:01 PM UTC-4, Jeffrey Arnold wrote:
>
>     Ben,
>
>     Could you add your discussion of the rules of when the user needs to manually adjust lp__ to the manual? Unless I am
>     missing something, the only discussion is in Ch 16.2 and it only talks about the case in which the LHS of ~ is an
>     expression or lp__ is adjusted directly.

We'll add a change of variables chapter to the
user's guide.

I want to try to make clear when you can transform
without a Jacobian and when's one necessary.

Which specific issue were you concerned about adding?
Just Ben's two-point post?

Just the points made in Ben's two point post with the objective of clearly defining when manual adjustments to the log posterior need to be made.  When I first read Ben's post, I went back through all my models to make sure I hadn't used a sampling statement with any transformed variables. Thankfully, I hadn't. But when writing them I didn't know that it would have been a problem,even though I was aware that if I made manual adjustments to the log posterior or transformed a variable in the model block I would likely have to add the log Jacobian.

 

>     And to make sure that I am understanding this correctly, for the following model
>
>         parameters {
>            real<lower=0> sigma;
>            ...
>         }
>         transformed parameters {
>            real<lower=0> tau;
>            tau <- sqrt(sigma);
>         }
>         model {
>            ...
>            tau ~ gamma(0.1, 0.1);
>            ...
>         }
>
>     requires an extra statement in the model block to adjust the log-posterior
>
>         lp__ <-  lp__ + log(0.5 * pow(sigma, -0.5));

I'd suggest putting it in the transformed parameters block.

Why in the transformed parameters block? Is this for style or efficiency? 

 

And that is the derivative you need.  I got SO confused
with the directionality of these in our first implementation
of Stan, that I can feel everyone's pain dealing with this
very acutely.  The transform here (at least in the notation
of Gelman et al.'s BDA and our manual) is from tau to
sigma.  The inverse transfrom is from sigma to tau,
and it's the log Jacobian of the inverse that needs to be added.

And for sparing me that pain, I don't know if I'll ever be able to thank you guys

enough! :-)

 

>     Also this may be a terrible idea, but I'll throw it out there anyways. Since accumulation seems to be used often in
>     Stan, any chance
>     of adding the equivalent of =+ and =-?  I guess that would be <-+ and <-- (the later being definitely suboptimal
>     syntax).

Yes, that's on the to-do list.  The usual syntax
in C, Java, etc. for compound assignment has the operator
first.  The statement

    a += b;

It seems that there is a p=1 that I'll get something backwards when I comment.  

 

has the same effect on a as:

    a = a + b;

Syntactically and semantically it's no issue to add,
but I'm uncertain what to call it given how ugly
"+<-" looks.  But that's what it'll have to be.

An alternative would be:

   lp_increment(x);

to stand in for:

   lp <- lp + x;

Thoughts?

The lp_increment() syntax may be more flexible.  I think somewhere on the TODO list is a way to automatically calculate likelihoods / deviance. 

From the user's perspective, a reasonable implementation would be to add a logical flag on sampling statements and lp_increment()

in order to increment a likelihood variable in addition to lp__.  But I understand that this feature is nowhere near the top of the queue, and is probably 

a bad idea for reasons that i can't think of right now. 

 

> I also wouldn't necessarily agree that accumulation needs to be used "often". It is just that we have been getting a lot
> of questions from people with big multivariate models, and Stan's implementations of the relevant densities work for
> everything but are optimal for nothing. So, for now, the fastest way to use multivariate normal / t or inverse / plain
> Wishart distributions is with hand-optimized accumulation, but it won't be that way forever.

I want to have this available for sums in loops, too.

I generally have to accumulate a sum in a loop for calculating the log-likelihood in all my models. (that probably also explains my previous comment)

 - Jeff

[Bob Carpenter]

On 9/24/12 6:40 PM, Jeffrey Arnold wrote:
> comments inline
>
> On 9/23/12 8:58 PM, Ben Goodrich wrote:
> > On Sunday, September 23, 2012 8:01:01 PM UTC-4, Jeffrey Arnold wrote:

> Which specific issue were you concerned about adding?
> Just Ben's two-point post?
>
>
> Just the points made in Ben's two point post with the objective of clearly defining when manual adjustments to the log
> posterior need to be made. When I first read Ben's post, I went back through all my models to make sure I hadn't used a
> sampling statement with any transformed variables. Thankfully, I hadn't. But when writing them I didn't know that it
> would have been a problem,even though I was aware that if I made manual adjustments to the log posterior or transformed
> a variable in the model block I would likely have to add the log Jacobian.

OK -- will do.

I need to make sure a warning gets printed if anyone
tries to use a transformed parameter on the left
side of a sampling statement. I think that would go
a long way to helping people out.

> > requires an extra statement in the model block to adjust the log-posterior
> >
> > lp__ <- lp__ + log(0.5 * pow(sigma, -0.5));
>
> I'd suggest putting it in the transformed parameters block.
>
>
> Why in the transformed parameters block? Is this for style or efficiency?

Just style. The Jacobian adjustment only depends on
the transformed parameter declarations. In general,
it's better to use things near where they're declared
for program clarity.

Now, where we fall down on that big time is by requiring
variable declarations in a block at the top!

> An alternative would be:
>
> lp_increment(x);
>
> to stand in for:
>
> lp <- lp + x;
>
> Thoughts?
>
>
> The lp_increment() syntax may be more flexible.

I'm not sure what you mean by flexibility, but
given that

lp_increment(x)
=== lp__ +<- x
=== lp__ <- lp__ + x;

the operator pileup +<- is definitely more general,
but the basic assignment <- is most general.

> I think somewhere on the TODO list is a way to automatically calculate
> likelihoods / deviance.

For that, we'd need to break the model apart into
a prior and likelihood component. That would actually
be pretty easy syntactically and in terms of the
translation, but then it'd be up to the user to get
the right things in the right places.

> From the user's perspective, a reasonable implementation would be to add a logical flag on sampling statements and
> lp_increment()
> in order to increment a likelihood variable in addition to lp__. But I understand that this feature is nowhere near the
> top of the queue, and is probably
> a bad idea for reasons that i can't think of right now.

I like having two blocks, but one way we originally
thought about doing declarations was with flags.

- Bob