Slow sampling fitting an auction model

Slow sampling fitting an auction model

Hi all,

I am working on an auction model and the sampling process in Stan is very slow. It took my laptop over 19 hours to fit 4 MCMC chains with 200 iterations for each. Moreover, all transitions after warm up exceeded the maximum treedepth. I’m wondering whether the usage of integral in the modeling dramatically slow down the sampling or not. If yes, how to improve the performance. Otherwise, what can be the potential problem that slow down the sampling.

The auction model I am working on, in brief, models the bid b_{i} of bidder i by the vector of explanatory variables X_i, the number of competitors J_i -1, and the cost c_{i} through

b_i = c_i + \frac{\int_{c_i}^{\bar{c}} [1-\mbox{F}(c; X_i^\top\beta, \sigma^2)]^{J_i - 1} dc}{[1-\mbox{F}(c_i; X_i\beta, \sigma^2)]^{J_i - 1}} + \epsilon_i\;,\; \epsilon_i \sim \mbox{N}(0, 1)\;,

where \mbox{F}(c; X_i^\top\beta, \sigma) is the cdf of a truncated lognormal distribution \mbox{lognormal}(X_i\beta, \sigma^2)\mbox{T}[0, \bar{c}] defined with hyperparameters \beta, \sigma^2 and the upper bound \bar{c}. The cost c_i is hierarcially modeled by the truncated lognormal distribution

c_i \sim \mbox{lognormal}(X_i\beta, \sigma^2)\mbox{T}[0, \bar{c}]\;.

The target is to obtain the posterior inference of cost c_i for all bidders and hyperparameters \beta and \sigma.

The following lists the parameter and model blocks in Stan code. The complete Stan code is in “generative_model.stan” generative_model.stan (3.0 KB) . One can reproduce the simulation by running “auction_proj.R”.auction_proj.R (1.3 KB) , which calls data.R (3.6 KB) to generate data.

parameters {
  vector[K] beta;
  real<lower=0> sigma;
  vector<lower=0, upper = c_max>[N] costs;

model {
  vector[N] inv_hazards;
  // priors
  to_vector(beta) ~ normal(0,1);
  sigma ~ normal(1,0.01);
  for(i in 1:N){
    costs[i] ~ lognormal(X[i] * beta, sigma) T[,c_max]; //dgp cost distribution
  // likelihood
  inv_hazards = getGenerativeTotalAuctionHazard(
                           J, // num bidders in each auction
  (bids - costs) ~ normal(inv_hazards , 1);



Lu Zhang

data.R (3.6 KB) auction_proj.R (1.3 KB) generative_model.stan (3.0 KB)

1 Like

@charlesm93 @bgoodri might you have a few mins to take a look/have any ideas for things to try?

Update, fixed a bug in attachment data.R

The model is likely slow because of the large treedepth. I’m guessing your posterior geometry is relatively flat, so your chain “particle” is given a shove, but never makes a U-turn. Moreover, your model may not be very identifiable, either because your data does not inform your parameters or something is degenerate. Can you share some pairs plot?

Just lookin at the code it’s hard to tell what might be going wrong. I’m not sure if your super tight prior on sigma causes issue. I’d consider fixing sigma or widening the prior.

I’m not too familiar with integrate_1d, and I’m not sure how well its performance has been tested. Since the bonds of your integral are finite, you should be able to rewrite your integral using an ODE solver.

Maybe generate some synthetic data from a model where the integral has a closed form and then substitute that into the model, that would tell you whether or not the problem is caused specifically by the integration.

Hi Charles,

Thank you so much! I couldn’t agree more that the large treedepth is the main reason of the slow sampling. But I doubt this model has a extremely flat posterior geometry. I made up a simulation using the same code except that I replaced the integral with something like c*(1-F)^{J}. I found the sampling for the made-up example was fast. The sampling process for 4 chains with each 1000 iters cost around 200 seconds. Of course, some of the R-hat of the made-up example were large, that probably indicates identifiable problem. But the speed are so different that makes me feel the problem may comes from the integral.

Also, I have tested with different settings, prior with larger variance, prior for the variance of noise, fixing more parameters. All of them are slow. I also tried the marginalized version, i.e., I marginalize out the cost and just focus on the inference of hyperparameters \beta, \sigma. But my model can’t even be initialized at the true value of the parameters.

I think maybe an ODE solver would be helpful. I will have a try and let you know. Thanks again for your help!

Yes, I did something like what you mentioned. I made up a simulation using the same code except that I replaced the integral with something like c*(1-F)^{J}. The sampling process for 4 chains with each 1000 iters cost around 200 seconds. That’s also one of the reason why I think the problem might be caused by the integration.

1 Like

Maybe now replace your closed form with the integration of it’s derivative? If that performs extremely slowly then it has to be the integration step.

I don’t know too much about stan’s internals, but I’m not sure it makes sense for calculating an integral to consume that much extra time – unless the integrand is really wild, it’s not that expensive of a computation in 1 dimension and the gradient shouldn’t be much more difficult than evaluating the integral itself since you can move the gradient inside the integral with Leibniz’s rule.

What is the size of \beta? If \beta is an n-vector then I think it would need to evaluate n integrals to get at the gradient so if n is really large maybe it could get pretty slow.

1 Like

That is a brilliant idea! I just tried it. I replaced

b_i = c_i + c_b(1-F(c_i))^{J - 1} +\epsilon_i


b_i = c_i +\int_{c_i}^{\bar{c}} c_b(J-1)(1-F(u))^{J - 2}f(u) du +\epsilon_i

And here is the comparison of the time estimated by stan:

By the way, the length of \beta is 3.

1 Like

Interesting. I’m not sure about where to go from there tbh. Seems unreasonably slow to me – maybe open another tread on the performance of integration? If the dimension is low, it may be plausible to evaluate the integral offline at the relevant points and then feed that as data to stan and do some interpolation between points – basically a lookup table for the integral? Not sure if that’s a good idea or not. Or maybe reduce the error tolerance for integration?