optimizing hierarchical non-linear regression

[Steve Bellan]

Hi all, I'm running a hierarchical non-linear regression with 3D multinormal random effects. I adapted the code from section 5.12 of the STAN manual for this task, hoping to maintain as much speed as possible since the data set has 70K individuals and 450K observations. It's working but, as expected, pretty slow in that it takes ~5 minutes per HMC iteration (it still converges in a reasonable # of iterations). I'm wondering if there are any other obvious tricks to speed things up here. For instance, I get the feeling that the for (n in 1:Nobs) loop could be vectorized to be faster but am not sure whether that's true or, if so, how. Any help is much appreciated!!

Thanks,

Steve

data {
  int<lower=0> Nobs; // # observations
  int<lower=1> Npat; // # individuals
  int<lower=1> NcovsPat; // # individual-level predictors (1, sex)
  int<lower=1> NasymParm; // # parameters for non-linear model
  int<lower=1,upper=Npat> jj[Nobs]; // individual for observation
  vector[Nobs] time;    // time of follow up
  matrix[Npat,NcovsPat] Upat;    // individual-level predictors [1, sex]_jj
  vector[Nobs] y;    // outcome
}
parameters {
  matrix[NasymParm,Npat] z;
  cholesky_factor_corr[NasymParm] CovsPatOmega; // correlation matrix for predictor-coefficients
  vector<lower=0>[NasymParm] tau; // prior scale on covar-var matrix
  matrix[NcovsPat,NasymParm] gamma; // individual-level coefficients (1,sex)X(Asym,R0,lrc)
  real<lower=0> sigma; // obs-level error
}
transformed parameters {
  matrix[Npat,NasymParm] beta; // individual-level predictors, each row is a NasymParm-length vector
  vector[Npat] Asym;      // Asymptote
  vector[Npat] R0;      // Initial value
  vector[Npat] lrc;      // log rate of change
  matrix[NasymParm,NasymParm] Sigma_beta; // reparameterize as cov-var matrix
  beta <- Upat * gamma + (diag_pre_multiply(tau,CovsPatOmega) * z)';
  Asym <- col(beta,1);
  R0 <- col(beta,2);
  lrc <- col(beta,3);
  Sigma_beta <- diag_pre_multiply(tau,CovsPatOmega) * diag_pre_multiply(tau,CovsPatOmega)'; // Recalculate as covar-var matrix
}
model {
  vector[Nobs] Mod_jj;
  for (n in 1:Nobs) {
    Mod_jj[n] <- Asym[jj[n]] + (R0[jj[n]] - Asym[jj[n]]) * exp(-exp(lrc[jj[n]]) * time[n]);
  }
  y ~ normal(Mod_jj, sigma);
  tau ~ cauchy(0,2.5);
  to_vector(z) ~ normal(0,1);
  CovsPatOmega ~ lkj_corr_cholesky(2);
  to_vector(gamma) ~ normal(0,5);
}

[Steve Bellan]

The problem is more severe than I thought since I posted this. It took > 2 hours for 30 iterations and then gave a segfault from C stack overflow at the end (see below). It never generated stan sampling object after all that. Any guidance would be very much appreciated:

...

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Iteration:  6 / 30 [ 20%]  (Warmup)
Iteration:  9 / 30 [ 30%]  (Warmup)
Iteration: 12 / 30 [ 40%]  (Warmup)
Iteration: 15 / 30 [ 50%]  (Warmup)
Iteration: 16 / 30 [ 53%]  (Sampling)
Iteration: 18 / 30 [ 60%]  (Sampling)
Iteration: 21 / 30 [ 70%]  (Sampling)
Iteration: 24 / 30 [ 80%]  (Sampling)
Iteration: 27 / 30 [ 90%]  (Sampling)
Iteration: 30 / 30 [100%]  (Sampling)
#  Elapsed Time: 1466.86 seconds (Warm-up)
#                6163.34 seconds (Sampling)
#                7630.2 seconds (Total)


Error: segfault from C stack overflow

[Daniel Lee]

Hi Steve,

What interface are you using? And what version?

Are you able to share the data? If so, I'll try to reproduce the problem on my end.

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

[Steve Bellan]

That'd be awesome, Daniel! I was actually in the middle editing the last post to have attached analysis in reproducible form but got interrupted. Here it is! The scripts are attached and the data are here since they're too big to attach. They were simulated from the same model as that being fit.

I'm using rstan 2.6.0 on R 3.1.1 on OS X 10.10.2.

The model fits fine to simulated data sets with smaller sizes. It was just when I scaled up from 5K to 70K individuals that the segfault occurred.

[Bob Carpenter]

You started with the process we recommend, so it's great
you've already figured out it's not a bug in the model.

The Stan code's very tight, but I can make a couple
suggestions.

First, move the declaration and definition of
Sigma_beta to the generated_quantities block because you
don't use them in the model. This will be a big savings
if the matrix is large.

Second, the expression

exp(-exp(lrc[jj[n]])

is just an underflow waiting to happen. It'll basically underflow
if lrc > 6.5 or so. Is there a way to express this not on
the exp(-exp()) scale?

It's odd that there'd be a stack overflow. That indicates
recursion going beyond the max depth that C++ can handle. We
wrote everything in Stan iteratively to avoid this issue.

We can run in CmdStan to see if it's an R or Stan issue. Daniel
and I are both on the road now, so I don't have anywhere I can
try a long-running model.

- Bob

> <asympStan.R><asympStan.stan>

[Steve Bellan]

Thanks Bob. Very helpful. I moved Sigma_beta to the generated quantities but it only improved speed on a smaller data (1K individuals, 6K observations) set by about 3%. This isn't too surprising I think since that's only a 3x3 matrix. But still very useful for me to note for future work.

As for exp(- exp(lrc) * time) and underflow issues, I'd prefer to keep lrc (log rate of change) to be sampled on a log scale since, as a rate, that is the more meaningful scale to sample it from a normal distribution (it's multinormally distributed with the other parameters). I plotted the model for lrc =seq(-3, 3, len=9) on both scales to double check that the spacing is indeed problematic on the un-logged scale. As you can see from the attached figure, indeed it is.

But in thinking about this I realized that I have ridiculous priors right now and I wonder if that's part of the problem. I kept

to_vector(gamma) ~ Normal(0,5)

from the manual code. But the first row of my gamma matrix are intercepts for each of the non-linear model parameters. And I know a priori that R0 should be in the range of of 2-6 and Asym should be 5-7. lrc is probably in the range of -2 to 2. So I changed the priors to:

  gamma[1,1] ~ normal(6, 5);
  gamma[1,2] ~ normal(4, 5);
  gamma[1,3] ~ normal(0, 2);
  gamma[2] ~ normal(0,5);

This did not speed up the HMC iterations at all on the smaller test set (1K indiv). But I'm hoping that by narrowing the prior on lrc to normal(0,2) it would reduce the chance of underflow, though I'm not sure if I'm thinking about this clearly.

I'll rerun it on a full-sized (ie size of the eventual real data) data set tonight & report what happens in the morning.

[Bob Carpenter]

> On Jun 26, 2015, at 12:04 AM, Steve Bellan <steve....@gmail.com> wrote:
>
> Thanks Bob. Very helpful. I moved Sigma_beta to the generated quantities but it only improved speed on a smaller data (1K individuals, 6K observations) set by about 3%. This isn't too surprising I think since that's only a 3x3 matrix. But still very useful for me to note for future work.
>

> As for exp(- exp(lrc) * time) and underflow issues, I'd prefer to keep lrc (log rate of change) to be sampled on a log scale since, as a rate, that is the more meaningful scale to sample it.

Sorry --- I wasn't suggesting changing the model, just trying to
find some algebra which would avoid the double exp. I'm not even
sure it's possible.

> I plotted the model for lrc =seq(-3, 3, len=9) on both scales to double check that the spacing is indeed problematic on the un-logged scale. As you can see from the attached figure, indeed it is.
>
> But in thinking about this I realized that I have ridiculous priors right now and I wonder if that's part of the problem.

They can be a big problem, especially in problems with
overdispersion (which includes a lot of latent-parameter
models and the negative binomial).

> I kept
> to_vector(gamma) ~ Normal(0,5)

> from the manual code. But the first row of my gamma matrix are intercepts for each of the non-linear model parameters. And I know a priori that R0 should be in the range of of 2-6 and R0 and Asym should be 5-7. lrc is probably in the range of -2 to 2. So I changed the priors to:

> gamma[1,1] ~ normal(6, 5);
> gamma[1,2] ~ normal(4, 5);
> gamma[1,3] ~ normal(0, 2);
> gamma[2] ~ normal(0,5);
> This did not speed up the HMC iterations at all on the smaller test set (1K indiv). But I'm hoping that by narrowing the prior on lrc to normal(0,2) it would reduce the chance of underflow, though I'm not sure if I'm thinking about this clearly.

What it'll do is keep things closer to 0.

Also, running something more constrained than those Cauchy
priors can help, too --- they allow big outliers. So you might
want to see in the traceplots or with max or something whether
some parameters are going to extreme values.

> I'll rerun it on a full-sized (ie size of the eventual real data) data set tonight & report what happens in the morning.

Good luck!

- Bob

[Steve Bellan]

Update: it's been running for 8+ hours and is on iteration 20 of 50. I'll let it keep going but it seems even slower than before! No idea how moving Sigma_beta or changing the priors could have done that. 

Any additional thoughts or help through the reproducible example when you guys finish traveling would be very much appreciated. 

[Steve Bellan]

It finished 50 iterations at 3.3 hours. My computer apparently went to sleep overnight, so wasn't as long as I thought. However, still suffering from stack overflow even with the new narrower prior on lrc. The C segfault didn't crash R, but when I tried to access the fitted object R had it's own segmentation fault.

...
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Iteration:  5 / 50 [ 10%]  (Warmup)
Iteration: 10 / 50 [ 20%]  (Warmup)
Iteration: 15 / 50 [ 30%]  (Warmup)
Iteration: 20 / 50 [ 40%]  (Warmup)
Iteration: 25 / 50 [ 50%]  (Warmup)
Iteration: 26 / 50 [ 52%]  (Sampling)
Iteration: 30 / 50 [ 60%]  (Sampling)
Iteration: 35 / 50 [ 70%]  (Sampling)
Iteration: 40 / 50 [ 80%]  (Sampling)
Iteration: 45 / 50 [ 90%]  (Sampling)
Iteration: 50 / 50 [100%]  (Sampling)
#  Elapsed Time: 5242.04 seconds (Warm-up)
#                6565.55 seconds (Sampling)
#                11807.6 seconds (Total)

Error: segfault from C stack overflow

> > > > > > > > 11807/3600

[1] 3.279722

> traceplot(fit1, 'gamma', ncol = 3, nrow = 2)

Process R segmentation fault: 11 at Fri Jun 26 12:41:44 2015

[Jiqiang Guo]

The segfault might be due to R's inability to create a data frame with too many columns.  How many parameters (including transformed parameters) do you have?

A similar issue that a segfault happened: https://github.com/stan-dev/rstan/issues/121

Jiqiang 

On Fri, Jun 26, 2015 at 1:41 PM, Steve Bellan <steve....@gmail.com> wrote:

It finished 50 iterations at 3.3 hours. My computer apparently went to sleep overnight, so wasn't as long as I thought. However, still suffering from stack overflow even with the new narrower prior on lrc.

...
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.


Iteration:  5 / 50 [ 10%]  (Warmup)
Iteration: 10 / 50 [ 20%]  (Warmup)
Iteration: 15 / 50 [ 30%]  (Warmup)
Iteration: 20 / 50 [ 40%]  (Warmup)
Iteration: 25 / 50 [ 50%]  (Warmup)
Iteration: 26 / 50 [ 52%]  (Sampling)
Iteration: 30 / 50 [ 60%]  (Sampling)
Iteration: 35 / 50 [ 70%]  (Sampling)
Iteration: 40 / 50 [ 80%]  (Sampling)
Iteration: 45 / 50 [ 90%]  (Sampling)
Iteration: 50 / 50 [100%]  (Sampling)
#  Elapsed Time: 5242.04 seconds (Warm-up)
#                6565.55 seconds (Sampling)
#                11807.6 seconds (Total)


Error: segfault from C stack overflow
> > > > > > > > 11807/3600
[1] 3.279722

[Steve Bellan]

Sounds plausible. I have 70K individuals with 3 random effects each, so that's 210,028 parameters total I believe. It didn't occur to me that it was storing all those random effects from every iteration. Looking at your linked thread, I'll try to specify to save every parameter except those random effects for now.

I'm almost set up to be running rstan on an XSEDE HPC, so I wonder if the additional memory there will help. But if it's an R issue then I guess it shouldn't.

[Andrew Gelman]

Once you have this sort of informative prior, we recommend you se inits in that range too.

On Jun 26, 2015, at 12:04 AM, Steve Bellan <steve....@gmail.com> wrote:

But in thinking about this I realized that I have ridiculous priors right now and I wonder if that's part of the problem. I kept

to_vector(gamma) ~ Normal(0,5)

from the manual code. But the first row of my gamma matrix are intercepts for each of the non-linear model parameters. And I know a priori that R0 should be in the range of of 2-6 and R0 and Asym should be 5-7. lrc is probably in the range of -2 to 2. So I changed the priors to:

[Steve Bellan]

I got 20 iterations to run in an hour without saving all those pesky random effects, no segfault! Thanks for the initial value tip Andrew, I'll make sure to do that as well. Now to try a longer run. Thanks for all the help!

S

[Michael Betancourt]

Assuming reasonable mixing (say 5 effective samples within those 20 samples),

5 effective samples per hour for 200,000+ parameters is pretty awesome!

[Steve Bellan]

Indeed it is pretty impressive! 

Another quick Q: would it be worthwhile to reparameterize my model so that all normal distributions are standard normal? I've already got the Cholesky reparameterization going on for the multinormal, but what about for my gamma priors 

  gamma[1,1] ~ normal(6, 5);
  gamma[1,2] ~ normal(4, 5);
  gamma[1,3] ~ normal(0, 2);
  gamma[2] ~ normal(0,5);

and the following line:

  y ~ normal(Mod_jj, sigma);

it seems like changing these to normal(0,1)'s is too trivial to affect the parameterization of the model & its speed, but I'm not sure if there's something I'm missing here about the virtues of sampling std normals.

[Steve Bellan]

So actually the # of effective samples is pretty bad when I have the full-sized data set (70K indiv X 6 obs/indiv). I've pasted the results of 2 chains of 100 iterations below. Things seem to work a bit better when I reduce the data set to 10K indiv X 6obs/indiv, with far better mixing. I don't understand why the # n_eff / iteration should go down so quickly with increasing data set size. How can I maintain high mixing efficiency at this greater data set size?

70K indiv x 6 obs/indiv with 2 chains of 100 iterations:

                   n_eff  Rhat
gamma[1,1]          3.05  4.88
gamma[1,2]          3.35  2.40
gamma[1,3]          3.34  4.25
gamma[2,1]          7.13  1.19
gamma[2,2]         47.30  1.01
gamma[2,3]         68.80  1.13
Sigma_beta[1,1]     1.01 30.50
Sigma_beta[1,2]     4.26  1.86
Sigma_beta[1,3]     1.03 13.50
Sigma_beta[2,1]     4.26  1.86
Sigma_beta[2,2]     1.06  5.04
Sigma_beta[2,3]     3.14  2.66
Sigma_beta[3,1]     1.03 13.50
Sigma_beta[3,2]     3.14  2.66
Sigma_beta[3,3]     3.08  5.55
sigma               1.16  7.50
tau[1]              1.01 60.00
tau[2]              1.10  4.48
tau[3]              3.40  3.79
CovsPatOmega[1,1] 100.00   NaN
CovsPatOmega[1,2] 100.00   NaN
CovsPatOmega[1,3] 100.00   NaN
CovsPatOmega[2,1]   1.13  5.28
CovsPatOmega[2,2]   1.93  2.86
CovsPatOmega[2,3] 100.00   NaN
CovsPatOmega[3,1]   1.07  4.14
CovsPatOmega[3,2]   3.90  1.89
CovsPatOmega[3,3]   4.95  1.54
lp__                1.14  8.07

(note that the n_eff are non-integer because I've extracted the summary matrix from the stanfit object so it no longer rounds)

10K indiv x 6 obs/indiv with 2 chains of 100 iterations 

  

                  n_eff Rhat
gamma[1,1]          7.61 1.27
gamma[1,2]         14.30 1.16
gamma[1,3]          8.51 1.22
gamma[2,1]        100.00 1.02
gamma[2,2]         45.80 1.03
gamma[2,3]        100.00 1.02
Sigma_beta[1,1]     6.32 1.29
Sigma_beta[1,2]     7.37 1.30
Sigma_beta[1,3]    15.60 1.11
Sigma_beta[2,1]     7.37 1.30
Sigma_beta[2,2]     4.19 1.57
Sigma_beta[2,3]    17.30 1.03
Sigma_beta[3,1]    15.60 1.11
Sigma_beta[3,2]    17.30 1.03
Sigma_beta[3,3]     7.41 1.28
sigma               5.87 1.35
tau[1]              6.80 1.27
tau[2]              4.37 1.54
tau[3]              7.29 1.30
CovsPatOmega[1,1] 100.00  NaN
CovsPatOmega[1,2] 100.00  NaN
CovsPatOmega[1,3] 100.00  NaN
CovsPatOmega[2,1]   4.90 1.53
CovsPatOmega[2,2]   6.26 1.33
CovsPatOmega[2,3] 100.00  NaN
CovsPatOmega[3,1]   4.22 1.20
CovsPatOmega[3,2]   3.57 1.97
CovsPatOmega[3,3]   4.69 1.48
lp__                6.20 1.32

On Saturday, June 27, 2015 at 3:40:43 AM UTC-5, Michael Betancourt wrote:

[Bob Carpenter]

Well, those clearly haven't converged. Are the chains close
and just not mixing, or are they clearly stuck in different
modes? That is, does it look like running longer might help?

I don't know if you're getting any divergences here (typically
indicated by warmup messages, but you can also pull them out
of the fit in RStan and they're printed in CmdStan), but if you
are, you probably want to start with a smaller step size and
up the adapt_delta target above its default of 0.8 (perhaps as
high as 0.99). That'll keep step size down and allow the gradients
to do a better job of mimicking the curvature. But it'll slow
down the iterations (at the benefit of getting higher n_eff per
iteration).

The problem with increasing data is that it can increase posterior
curvature along narrow ridges, making it harder for the Hamiltonian
simulation to find a large step size to make large moves. Sometimes
this is good and you actually get sublinear growth with data. Sometimes
it's bad and constrains the parameters and makes them correlated in
a way that makes it much harder to move in Euclidean space with just
a gradient. So you'll probably see the number of leapfrog steps going up and
the step size going down as you get more data to try to compensate.

- Bob

[Michael Betancourt]

You could be running into the circumstance where the data are
very informative and you want to run a centered parameterization,
as discussed in the HMC for Hierarchical Models paper.

[Krzysztof Sakrejda]

On Monday, June 29, 2015 at 3:13:32 PM UTC-4, Bob Carpenter wrote:

Well, those clearly haven't converged.  Are the chains close
and just not mixing, or are they clearly stuck in different
modes?  That is, does it look like running longer might help?

I don't know if you're getting any divergences here (typically
indicated by warmup messages, but you can also pull them out
of the fit in RStan and they're printed in CmdStan), but if you
are, you probably want to start with a smaller step size and
up the adapt_delta target above its default of 0.8 (perhaps as
high as 0.99).  That'll keep step size down and allow the gradients
to do a better job of mimicking the curvature.  But it'll slow
down the iterations (at the benefit of getting higher n_eff per
iteration).

The problem with increasing data is that it can increase posterior
curvature along narrow ridges, making it harder for the Hamiltonian
simulation to find a large step size to make large moves.  Sometimes
this is good and you actually get sublinear growth with data.  Sometimes
it's bad and constrains the parameters and makes them correlated in
a way that makes it much harder to move in Euclidean space with just
a gradient.  So you'll probably see the number of leapfrog steps going up and
the step size going down as you get more data to try to compensate.

I've also seen this lead to non-convergence when some chains find the narrow ridges and others never do (and explore a much broader parameter space).  You can ususally tell that from trace plots or just looking at means/variances for multiple chains.

Krzysztof
 

[Steve Bellan]

Bob: I've attached trace plots for both the 10K & 70K models. It looks to me like some parameters have the chains coming close but not mixing (gamma, sigma), but others (the correlation/covariance matrix CovsPatsOmega/Sigma_beta) are not coming close in the amount of time the chain runs. So it does seem like there is a thin likelihood ridge for the fixed effects (gamma) but not necessarily for their correlations. I'm not sure what this means practically. The problem is clearly worse for the 70K model.

Michael: based on your advice I've implemented a centered version of the code below.  I went through the 2013 paper but must admit that I still don't quite follow the rationale for non-centering for small data sets and centering for large data sets. It is unclear how much better/worse the centered version does on small-intermediate size data sets thus far. Some parameters have greater n_eff, some have smaller. I suspect it's hard to tell from only 2 chains. I'm in the process of getting this workflow setup on HPC and will then run batch fits with & without centering, & for a range of deltas (.8,.9,.95,.99) with 24 chains each run, and will then report the results.

My centered model is:

data {
  int<lower=0> Nobs; // # observations
  int<lower=1> Npat; // # individuals

  int<lower=1> NcovsPat; // # individuals-level predictors (1, sex)

  int<lower=1> NasymParm; // # parameters for non-linear model

  int<lower=1,upper=Npat> jj[Nobs]; // individuals for observation

  vector[Nobs] time;    // time of follow up

  matrix[Npat,NcovsPat] Upat;    // individuals-level predictors [1, sex]_jj

  vector[Nobs] y;    // outcome
}
parameters {

  corr_matrix[NasymParm] Omega; // correlation matrix for predictor-coefficients

  vector<lower=0>[NasymParm] tau; // prior scale on covar-var matrix

  matrix[NcovsPat,NasymParm] gamma; // individuals-level coefficients (1,sex)X(Asym,R0,lrc)
  vector[NasymParm] beta[Npat];    // each individual's asymptotic parameters

  real<lower=0> sigma; // obs-level error
}

model {
  matrix[NasymParm,NasymParm] Sigma_beta; // covar-var matrix
  vector[Nobs] Mod_jj;
  vector[NasymParm] u_gamma_t[Npat];
  Sigma_beta <- quad_form_diag(Omega,tau);
  tau ~ cauchy(0,2.5);
  Omega ~ lkj_corr(2);

  gamma[1,1] ~ normal(6, 5);
  gamma[1,2] ~ normal(4, 5);
  gamma[1,3] ~ normal(0, 2);
  gamma[2] ~ normal(0,5);

  for(j in 1:Npat)
    u_gamma_t[j] <- (Upat[j] * gamma)';
  beta ~ multi_normal(u_gamma_t, Sigma_beta);

  for (n in 1:Nobs) {

    Mod_jj[n] <- beta[jj[n],1] + (beta[jj[n],2] - beta[jj[n],1]) * exp(-exp(beta[jj[n],3]) * time[n]);
  }
  y ~ normal(Mod_jj, sigma);
}

[Michael Betancourt]

> I've also seen this lead to non-convergence when some chains find the narrow ridges and others never do (and explore a much broader parameter space). You can ususally tell that from trace plots or just looking at means/variances for multiple chains.

Formally, this isn’t an issue with convergence but rather one of mixing.
Every chain should be exploring both the narrow ridges and the
broader bulk.

[Krzysztof Sakrejda]

Good point, I was thinking of the context where the broader bulk has a much worse log posterior density and so is never explored once the ridge is found but it must happen in less extreme situations too.
 

[Bob Carpenter]

You want to be thinking in terms of posterior mass more than
density.

For example, consider beta(0.9,0.9), which has unbounded
density as y approaches 0 or 1, even though roughly 90% of
the mass is in (0.05, 0.95).

- Bob

[Bob Carpenter]

I didn't get the traceplot attachments.

- Bob

[Steve Bellan]

Here they are:

[Bob Carpenter]

The 10K version seems to have converged, but the 70K obviously
hasn't. And variables like tau[1] seem stuck in different modes.

All I can suggest from looking at this is to (a) start with
a smaller step size, (b) target a higher acceptance rate, and (c)
run longer.

If it's truly multimodal it won't help---maybe just plotting the
after-convergence bit of lp__ would show the densities are way off.
Hard to tell with the warmup skewing the scale.

- Bob

> On Jun 30, 2015, at 4:16 PM, Steve Bellan <steve....@gmail.com> wrote:
>
> Here they are:

>
> --
> 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/d/optout.

> <70Ktraceplot.pdf><10Ktraceplot.pdf>

>
>> On Jun 30, 2015, at 1:51 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>>
>> I didn't get the traceplot attachments.
>>
>> - 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/d/optout.
>>
>> --

>> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
>> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/XuHdf1-QuKc/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Steve Bellan]

Thanks Bob. I doubt it would be multimodal since the fits to smaller data set seem to easily identify the posterior, I don’t see how increasing data size could lead to more posterior modes.

I’ll definitely run longer. As for (a) and (b), my impression was that increasing delta (the target acceptance rate) automatically reduces the step size (which is adaptively tuned by NUTS), and that therefore I can do both (a) and (b) just by tuning this parameter. Or were you suggesting to explicitly set the step size through some other way, like switching to an HMC algorithm?

[Krzysztof Sakrejda]

On Tuesday, June 30, 2015 at 2:51:21 PM UTC-4, Bob Carpenter wrote:

You want to be thinking in terms of posterior mass more than
density.  

For example, consider beta(0.9,0.9), which has unbounded
density as y approaches 0 or 1, even though roughly 90% of
the mass is in (0.05, 0.95).

 

Thanks for keeping me on my toes.  Lucky for me my problem was a Weibull-uniform mixture with reasonable parameter values (shape ~3) so I can sleep without checking but it does make it harder to just use the lpd directly when interpreting a generic model .  Krzysztof
 

[Bob Carpenter]

> On Jun 30, 2015, at 4:41 PM, Steve Bellan <steve....@gmail.com> wrote:
>
> Thanks Bob. I doubt it would be multimodal since the fits to smaller data set seem to easily identify the posterior, I don’t see how increasing data size could lead to more posterior modes.

It depends on the data and how it's drawn. What looks
like an "outlier" in a small data set might become a well-delineated
mode in a bigger data set. Especially if one mode's much
higher mass than the other.

> I’ll definitely run longer. As for (a) and (b), my impression was that increasing delta (the target acceptance rate) automatically reduces the step size (which is adaptively tuned by NUTS), and that therefore I can do both (a) and (b) just by tuning this parameter. Or were you suggesting to explicitly set the step size through some other way, like switching to an HMC algorithm?

Not quite. There is an initial step size (a), which then gets
adapted based on the target acceptance rate (b). Changing the
target acceptance rate has no effect on the initial step size.
Usually if you need a higher acceptance rate, you want to start
with a smaller step size so you don't run too far into the tail
during warmup by taking too-big steps. It'll mainly cut down on
warnings and make convergence smoother and less variable across
chains than anything else.

[Steve Bellan]

What value would be reasonable to try for stepsize or how do I decide what values to  try?

I assume this is in control=list(stepsize)

[Luc Coffeng]

The standard initial stepsize is 1.0. What I usually do is manually tweak to a point where I no longer get divergent iterations during sampling. You don't want to go low than necessary as it will make your sampling unnecessarily inefficient (NUTS treedepth goes up when stepsize is lower).