Re: model to speed up?

[Bob Carpenter]

On 1/8/13 9:35 PM, Jiqiang Guo wrote:
> Attached vdq2.stan is the model, samples1.csv is a sample
> of one chain, vdq2.txt is the print output for 4 chains, the
> other is the data

I've appended just the model and print
below so everyone else could see it.

> I looked at tree_depth, most of them are 9,10,11. (why we have 11, is the default maximum 10?)

Ack, depth 11 is a huge number -- that's 1000
leapfrog steps (2 ^ (depth-1))! It must be an off-by-one
problem on the bound for number of steps in the code.
I'll take a look.

If you are in fact, hitting the upper limit and getting
stuck, you may just be falling back to a random walk
and it might help actually increasing the maximum number
of leapfrog steps, as counterintuitive as that may sound.

It does appear that low variance may be an issue.
For instance, sigma_age has (0.1,1.3) as its
95% posterior interval, whereas sigma_gender is
even worse at (0.1,11). I guess this makes sense
because there's a lot of posterior uncertainty in
a hierarchical scale parameter estimated from four
coefficients.

The only suggestion I have on the model is to
reparameterize the Cauchy using its inverse
CDF on a uniform variate (as described in the latest
version of the manual on git). If it helps with the
posterior geometry of those scale parameters, it
may make a difference. It certainly won't speed up
the leapfrog steps any.

The intercept variable, a, didn't get a prior
and wound up with a 95% interval of (-6.4,5.9).

- Bob




data {
int<lower=1> Nobs;
int<lower=1> Ngender;
int<lower=1> Nage;
int<lower=1> Nrace;
int<lower=1> Neduc;
int<lower=1> Nfaminc;
int<lower=1> Nemploy;
int<lower=1> Nmarstat;
int<lower=1> Nstateid;
int<lower=1,upper=Ngender> gender[Nobs];
int<lower=1,upper=Nage> age4[Nobs];
int<lower=1,upper=Nrace> race4[Nobs];
int<lower=1,upper=Neduc> educ5[Nobs];
int<lower=1,upper=Nfaminc> faminc5[Nobs];
int<lower=1,upper=Nemploy> employ4[Nobs];
int<lower=1,upper=Nmarstat> marstat3[Nobs];
int<lower=1,upper=Nstateid> stateid[Nobs];
int<lower=0,upper=1> rvote[Nobs];
real<lower=0,upper=1> stateobama[Nstateid];
}


parameters {
real a;

real beta_gender[Ngender]; # 1
real<lower=0> sigma_gender;
real beta_age[Nage]; # 2
real<lower=0> sigma_age;
real beta_race[Nrace]; # 3
real<lower=0> sigma_race;
real beta_educ[Neduc]; # 4
real<lower=0> sigma_educ;
real beta_faminc[Nfaminc]; # 5
real<lower=0> sigma_faminc;
real beta_employ[Nemploy]; # 6
real<lower=0> sigma_employ;
real beta_marstat[Nmarstat]; # 7
real<lower=0> sigma_marstat;

real beta_state[Nstateid];
real z_state[Nstateid];
real<lower=0> sigma_beta_state;
real<lower=0> sigma_z_state;
}

model {
real mu[Nobs];
int state_i;
for (n in 1:Nobs) {
state_i <- stateid[n];
mu[n] <- a + beta_gender[gender[n]] + beta_age[age4[n]] + beta_educ[educ5[n]]
+ beta_race[race4[n]] + beta_faminc[faminc5[n]] + beta_marstat[marstat3[n]]
+ beta_employ[employ4[n]]
+ beta_state[state_i] * stateobama[state_i] + z_state[state_i];
}
rvote ~ bernoulli_logit(mu);

// prior
beta_gender ~ normal(0, sigma_gender); # 1
beta_age ~ normal(0, sigma_age); # 2
beta_race ~ normal(0, sigma_race); # 3
beta_educ ~ normal(0, sigma_educ); # 4
beta_faminc ~ normal(0, sigma_faminc); # 5
beta_employ ~ normal(0, sigma_employ); # 6
beta_marstat ~ normal(0, sigma_marstat); # 7

beta_state ~ normal(0, sigma_beta_state);
z_state ~ normal(0, sigma_z_state);
sigma_gender ~ cauchy(0, 5); # 1
sigma_age ~ cauchy(0, 5); # 2
sigma_race ~ cauchy(0, 5); # 3
sigma_educ ~ cauchy(0, 5); # 4
sigma_faminc ~ cauchy(0, 5); # 5
sigma_employ ~ cauchy(0, 5); # 6
sigma_marstat ~ cauchy(0, 5); # 7
sigma_beta_state ~ cauchy(0, 5);
sigma_z_state ~ cauchy(0, 5);
}



> library(rstan)
> options(width = 120)
> post <- read_stan_csv(c('samples1.csv',
+ 'samples2.csv',
+ 'samples3.csv',
+ 'samples4.csv'))
> print(post)
Inference for Stan model: samples1.
4 chains: each with iter=2000; warmup=1000; thin=1; 1000 iterations saved.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a -0.7 0.4 3.9 -6.4 -1.8 -0.8 0.3 5.9 111 1.0
beta_gender[1] 0.3 0.2 2.9 -4.5 -0.2 0.2 0.6 5.2 220 1.0
beta_gender[2] -0.1 0.2 2.9 -4.8 -0.6 -0.2 0.2 4.9 220 1.0
sigma_gender 2.2 0.2 4.3 0.1 0.4 1.0 2.5 11.0 495 1.0
beta_age[1] -0.2 0.0 0.3 -0.7 -0.3 -0.2 -0.1 0.4 1206 1.0
beta_age[2] -0.1 0.0 0.3 -0.6 -0.3 -0.1 0.0 0.4 1176 1.0
beta_age[3] 0.0 0.0 0.3 -0.4 -0.1 0.0 0.2 0.6 1174 1.0
beta_age[4] 0.3 0.0 0.3 -0.2 0.2 0.3 0.4 0.9 1198 1.0
sigma_age 0.4 0.0 0.3 0.1 0.2 0.3 0.5 1.3 1236 1.0
beta_race[1] 1.0 0.0 1.1 -1.3 0.4 1.0 1.5 3.4 1265 1.0
beta_race[2] -1.9 0.0 1.1 -4.1 -2.4 -1.9 -1.3 0.5 1273 1.0
beta_race[3] 0.3 0.0 1.1 -2.0 -0.3 0.3 0.8 2.7 1266 1.0
beta_race[4] 0.7 0.0 1.1 -1.6 0.1 0.7 1.2 3.1 1267 1.0
sigma_race 2.0 0.0 1.3 0.8 1.2 1.7 2.4 5.3 1851 1.0
beta_educ[1] 0.2 0.0 0.3 -0.3 0.1 0.2 0.3 0.7 515 1.0
beta_educ[2] 0.3 0.0 0.2 -0.2 0.1 0.3 0.4 0.7 507 1.0
beta_educ[3] 0.1 0.0 0.2 -0.3 0.0 0.1 0.3 0.6 454 1.0
beta_educ[4] -0.1 0.0 0.2 -0.6 -0.2 -0.1 0.0 0.4 499 1.0
beta_educ[5] -0.5 0.0 0.2 -1.0 -0.6 -0.5 -0.3 0.0 452 1.0
sigma_educ 0.5 0.0 0.3 0.2 0.3 0.4 0.5 1.3 913 1.0
beta_faminc[1] -0.2 0.0 0.1 -0.5 -0.3 -0.2 -0.1 0.1 720 1.0
beta_faminc[2] -0.1 0.0 0.1 -0.4 -0.2 -0.1 0.0 0.2 719 1.0
beta_faminc[3] 0.0 0.0 0.1 -0.2 0.0 0.0 0.1 0.3 721 1.0
beta_faminc[4] 0.1 0.0 0.1 -0.2 0.0 0.1 0.2 0.4 718 1.0
beta_faminc[5] 0.2 0.0 0.1 -0.1 0.1 0.2 0.2 0.5 728 1.0
sigma_faminc 0.2 0.0 0.2 0.1 0.1 0.2 0.3 0.7 1011 1.0
beta_employ[1] 0.0 0.0 0.1 -0.1 0.0 0.0 0.1 0.1 939 1.0
beta_employ[2] 0.0 0.0 0.1 -0.2 0.0 0.0 0.0 0.1 799 1.0
beta_employ[3] 0.0 0.0 0.1 -0.1 0.0 0.0 0.1 0.2 1043 1.0
beta_employ[4] -0.1 0.0 0.1 -0.2 -0.1 -0.1 0.0 0.0 940 1.0
sigma_employ 0.1 0.0 0.1 0.0 0.1 0.1 0.1 0.3 1308 1.0
beta_marstat[1] 0.0 0.3 2.3 -1.8 0.1 0.2 0.4 1.2 80 1.1
beta_marstat[2] -0.2 0.3 2.3 -2.0 -0.1 0.0 0.2 1.0 80 1.1
beta_marstat[3] -0.5 0.3 2.3 -2.3 -0.4 -0.2 -0.1 0.8 80 1.1
sigma_marstat 1.0 0.3 2.6 0.1 0.3 0.5 0.9 4.2 89 1.1
beta_state[1] 0.1 0.0 0.4 -0.6 -0.1 0.1 0.4 1.0 4000 1.0
beta_state[2] 0.6 0.1 0.6 -0.3 0.1 0.4 1.0 2.0 52 1.1
beta_state[3] 0.3 0.0 0.4 -0.5 0.0 0.2 0.6 1.2 232 1.0
beta_state[4] 0.1 0.0 0.4 -0.6 -0.1 0.0 0.3 0.9 4000 1.0
beta_state[5] -0.2 0.0 0.4 -1.0 -0.5 -0.2 0.0 0.4 214 1.0
beta_state[6] -0.1 0.0 0.3 -0.7 -0.3 0.0 0.1 0.6 4000 1.0
beta_state[7] -0.3 0.0 0.4 -1.1 -0.6 -0.3 0.0 0.3 122 1.0
beta_state[8] -0.6 0.1 0.5 -1.6 -1.0 -0.6 -0.2 0.1 37 1.1
beta_state[9] -0.4 0.1 0.4 -1.3 -0.7 -0.3 0.0 0.3 57 1.1
beta_state[10] -0.1 0.0 0.3 -0.7 -0.3 0.0 0.1 0.6 1282 1.0
beta_state[11] 0.3 0.0 0.4 -0.4 0.0 0.2 0.6 1.2 172 1.0
beta_state[12] -0.3 0.0 0.4 -1.1 -0.6 -0.3 0.0 0.3 162 1.0
beta_state[13] -0.3 0.0 0.4 -1.1 -0.5 -0.2 0.0 0.4 91 1.0
beta_state[14] 0.2 0.0 0.4 -0.6 -0.1 0.1 0.5 1.1 386 1.0
beta_state[15] -0.3 0.0 0.4 -1.1 -0.6 -0.3 0.0 0.4 81 1.0
beta_state[16] -0.1 0.0 0.3 -0.8 -0.3 0.0 0.1 0.7 4000 1.0
beta_state[17] 0.0 0.0 0.4 -0.7 -0.2 0.0 0.2 0.8 4000 1.0
beta_state[18] 0.0 0.0 0.4 -0.7 -0.2 0.0 0.2 0.8 4000 1.0
beta_state[19] 0.6 0.1 0.6 -0.3 0.1 0.4 1.0 1.9 84 1.1
beta_state[20] -0.4 0.1 0.4 -1.3 -0.8 -0.4 -0.1 0.2 55 1.1
beta_state[21] -0.1 0.0 0.3 -0.8 -0.3 -0.1 0.1 0.5 490 1.0
beta_state[22] -0.3 0.0 0.4 -1.2 -0.6 -0.3 0.0 0.4 131 1.0
beta_state[23] -0.2 0.0 0.4 -1.0 -0.5 -0.2 0.0 0.5 232 1.0
beta_state[24] -0.3 0.0 0.4 -1.2 -0.6 -0.2 0.0 0.4 115 1.0
beta_state[25] -0.1 0.0 0.3 -0.8 -0.3 0.0 0.1 0.6 4000 1.0
beta_state[26] 0.7 0.1 0.7 -0.2 0.1 0.5 1.3 2.3 37 1.1
beta_state[27] -0.1 0.0 0.4 -0.9 -0.3 -0.1 0.1 0.6 503 1.0
beta_state[28] 0.1 0.0 0.3 -0.6 -0.1 0.1 0.3 0.8 4000 1.0
beta_state[29] 0.1 0.0 0.4 -0.7 -0.2 0.0 0.3 0.9 4000 1.0
beta_state[30] 0.0 0.0 0.4 -0.7 -0.2 0.0 0.2 0.8 4000 1.0
beta_state[31] -0.3 0.0 0.4 -1.0 -0.5 -0.2 0.0 0.4 158 1.0
beta_state[32] -0.2 0.0 0.4 -1.0 -0.5 -0.2 0.0 0.4 207 1.0
beta_state[33] -0.1 0.0 0.3 -0.7 -0.3 -0.1 0.1 0.6 4000 1.0
beta_state[34] -0.1 0.0 0.3 -0.8 -0.3 -0.1 0.1 0.6 884 1.0
beta_state[35] -0.3 0.0 0.4 -1.1 -0.6 -0.3 0.0 0.3 104 1.1
beta_state[36] -0.2 0.0 0.4 -1.0 -0.4 -0.2 0.0 0.5 257 1.0
beta_state[37] 0.2 0.0 0.4 -0.6 -0.1 0.2 0.5 1.2 270 1.0
beta_state[38] -0.3 0.0 0.4 -1.0 -0.5 -0.2 0.0 0.4 171 1.0
beta_state[39] -0.2 0.0 0.4 -1.0 -0.5 -0.2 0.0 0.5 196 1.0
beta_state[40] -0.4 0.1 0.4 -1.3 -0.7 -0.3 -0.1 0.3 44 1.1
beta_state[41] 0.3 0.0 0.4 -0.5 0.0 0.2 0.5 1.2 286 1.0
beta_state[42] 0.0 0.0 0.4 -0.8 -0.2 0.0 0.2 0.7 4000 1.0
beta_state[43] 0.1 0.0 0.4 -0.5 -0.1 0.1 0.3 0.9 4000 1.0
beta_state[44] 0.2 0.0 0.4 -0.5 -0.1 0.1 0.4 1.0 4000 1.0
beta_state[45] 0.5 0.1 0.6 -0.4 0.0 0.3 0.9 1.9 50 1.1
beta_state[46] 0.0 0.0 0.3 -0.7 -0.2 0.0 0.2 0.7 4000 1.0
beta_state[47] -0.6 0.1 0.6 -1.7 -1.1 -0.6 -0.2 0.2 36 1.1
beta_state[48] -0.3 0.0 0.4 -1.1 -0.6 -0.2 0.0 0.4 66 1.0
beta_state[49] -0.2 0.0 0.4 -0.9 -0.4 -0.1 0.0 0.5 304 1.0
beta_state[50] 0.0 0.0 0.4 -0.8 -0.2 0.0 0.2 0.7 4000 1.0
beta_state[51] 0.4 0.0 0.5 -0.5 0.0 0.2 0.7 1.6 132 1.0
z_state[1] 0.2 0.0 0.2 -0.2 0.0 0.1 0.3 0.5 529 1.0
z_state[2] 0.5 0.0 0.3 0.0 0.3 0.6 0.8 1.0 47 1.1
z_state[3] 0.3 0.0 0.2 -0.1 0.1 0.3 0.4 0.7 61 1.0
z_state[4] 0.1 0.0 0.2 -0.3 0.0 0.1 0.2 0.4 630 1.0
z_state[5] -0.1 0.0 0.2 -0.5 -0.3 -0.2 0.0 0.2 426 1.0
z_state[6] 0.0 0.0 0.2 -0.4 -0.1 0.0 0.1 0.3 4000 1.0
z_state[7] -0.2 0.0 0.2 -0.6 -0.4 -0.2 -0.1 0.2 217 1.0
z_state[8] -0.4 0.0 0.4 -1.1 -0.7 -0.4 -0.1 0.2 92 1.1
z_state[9] -0.3 0.0 0.2 -0.7 -0.5 -0.3 -0.1 0.2 144 1.0
z_state[10] 0.0 0.0 0.2 -0.3 -0.1 0.0 0.1 0.3 4000 1.0
z_state[11] 0.3 0.0 0.2 -0.1 0.1 0.3 0.5 0.7 71 1.0
z_state[12] -0.2 0.0 0.2 -0.7 -0.4 -0.2 0.0 0.2 201 1.0
z_state[13] -0.2 0.0 0.2 -0.6 -0.3 -0.2 0.0 0.2 373 1.0
z_state[14] 0.2 0.0 0.2 -0.1 0.1 0.2 0.4 0.6 231 1.0
z_state[15] -0.2 0.0 0.2 -0.6 -0.4 -0.2 0.0 0.2 295 1.0
z_state[16] 0.0 0.0 0.2 -0.4 -0.1 0.0 0.1 0.3 4000 1.0
z_state[17] 0.1 0.0 0.2 -0.3 0.0 0.1 0.2 0.4 4000 1.0
z_state[18] 0.1 0.0 0.2 -0.3 0.0 0.1 0.2 0.4 4000 1.0
z_state[19] 0.5 0.0 0.3 0.0 0.3 0.6 0.8 1.0 46 1.1
z_state[20] -0.3 0.0 0.2 -0.7 -0.5 -0.3 -0.1 0.1 111 1.1
z_state[21] -0.1 0.0 0.2 -0.4 -0.2 -0.1 0.1 0.3 4000 1.0
z_state[22] -0.2 0.0 0.2 -0.6 -0.4 -0.2 -0.1 0.2 240 1.0
z_state[23] -0.1 0.0 0.2 -0.5 -0.3 -0.1 0.0 0.2 4000 1.0
z_state[24] -0.2 0.0 0.2 -0.6 -0.4 -0.2 -0.1 0.2 193 1.0
z_state[25] 0.0 0.0 0.2 -0.4 -0.1 0.0 0.1 0.3 4000 1.0
z_state[26] 0.7 0.1 0.4 0.0 0.4 0.7 0.9 1.2 35 1.1
z_state[27] -0.1 0.0 0.2 -0.4 -0.2 -0.1 0.0 0.3 4000 1.0
z_state[28] 0.1 0.0 0.2 -0.3 0.0 0.1 0.2 0.5 539 1.0
z_state[29] 0.1 0.0 0.2 -0.3 0.0 0.1 0.2 0.5 4000 1.0
z_state[30] 0.0 0.0 0.2 -0.3 0.0 0.0 0.2 0.4 4000 1.0
z_state[31] -0.2 0.0 0.2 -0.5 -0.3 -0.2 0.0 0.2 4000 1.0
z_state[32] -0.1 0.0 0.2 -0.5 -0.3 -0.1 0.0 0.2 491 1.0
z_state[33] 0.0 0.0 0.2 -0.4 -0.1 0.0 0.1 0.3 4000 1.0
z_state[34] 0.0 0.0 0.2 -0.4 -0.1 0.0 0.1 0.3 4000 1.0
z_state[35] -0.2 0.0 0.2 -0.7 -0.4 -0.2 -0.1 0.2 207 1.0
z_state[36] -0.1 0.0 0.2 -0.5 -0.2 -0.1 0.0 0.2 734 1.0
z_state[37] 0.2 0.0 0.2 -0.1 0.1 0.3 0.4 0.6 165 1.0
z_state[38] -0.2 0.0 0.2 -0.5 -0.3 -0.2 0.0 0.2 301 1.0
z_state[39] -0.1 0.0 0.2 -0.5 -0.3 -0.1 0.0 0.2 532 1.0
z_state[40] -0.3 0.0 0.2 -0.7 -0.5 -0.3 -0.1 0.2 116 1.0
z_state[41] 0.3 0.0 0.2 -0.2 0.1 0.3 0.4 0.6 74 1.0
z_state[42] 0.0 0.0 0.2 -0.4 -0.1 0.0 0.1 0.4 4000 1.0
z_state[43] 0.2 0.0 0.2 -0.2 0.0 0.2 0.3 0.5 317 1.0
z_state[44] 0.2 0.0 0.2 -0.2 0.1 0.2 0.3 0.5 265 1.0
z_state[45] 0.5 0.0 0.3 0.0 0.3 0.6 0.7 0.9 38 1.1
z_state[46] 0.0 0.0 0.2 -0.4 -0.1 0.0 0.1 0.4 4000 1.0
z_state[47] -0.5 0.1 0.3 -1.0 -0.7 -0.5 -0.2 0.1 42 1.1
z_state[48] -0.2 0.0 0.2 -0.6 -0.3 -0.2 -0.1 0.2 298 1.0
z_state[49] -0.1 0.0 0.2 -0.5 -0.2 -0.1 0.0 0.2 861 1.0
z_state[50] 0.0 0.0 0.2 -0.3 -0.1 0.0 0.1 0.3 4000 1.0
z_state[51] 0.4 0.0 0.2 -0.1 0.2 0.4 0.6 0.8 61 1.1
sigma_beta_state 0.5 0.0 0.2 0.1 0.3 0.5 0.7 0.9 26 1.1
sigma_z_state 0.3 0.0 0.1 0.1 0.3 0.3 0.4 0.5 29 1.1
lp__ -32176.0 5.1 32.0 -32215.1 -32196.4 -32184.0 -32165.0 -32090.9 39 1.1

Samples were drawn using NUTS2 at Wed Jan 02 05:46:01 2013.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
>

[Andrew Gelman]

J
I'd suggest starting with Bob's suggestions, then another idea is we could throw in boundary-avoiding priors for the sigma parameters, as a computational device to get things moving faster.
A

> --
>

[Andrew Gelman]

J

Actually, even before the boundary-avoiding priors, you can try running with the variance parameters set to fixed (and reasonable values). If _that_ is slow, then we have a different sort of computational problem. If that goes fast, then we can try gradually relaxing the priors until it starts to go wrong.

I disagree with Bob's suggestion that we should be increasing the max number of leapfrog steps. My intuition is that, if we're taking 1000 steps on this simple hierarchical regression, we need to figure out what's going wrong, rather than trying to solve it using brute force.

A

[Daniel Lee]

On Wed, Jan 9, 2013 at 11:44 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

I disagree with Bob's suggestion that we should be increasing the max number of leapfrog steps.  My intuition is that, if we're taking 1000 steps on this simple hierarchical regression, we need to figure out what's going wrong, rather than trying to solve it using brute force.

Bob's suggestion was for the number of leapfrog steps, not the number of iterations. Perhaps the trajectory isn't getting to the U-turn condition before it maxes out on the number of iterations (given the size of the step). If the model is hitting the boundary here, I think Bob's intuition is right. There may be some other fundamental computational problems with the model, as you point out.

--

[Jiqiang Guo]

I am trying all the options.  

In addition, Yuanjun pointed out under current model specification, the correlations between betas for a factor is strong too much:

> e <- extract(fit, pars = c('beta_gender', 'beta_educ', 'beta_age', 'beta_race'))

> cor(e$beta_gender)

          [,1]      [,2]

[1,] 1.0000000 0.9999785

[2,] 0.9999785 1.0000000

> cor(e$beta_educ)

          [,1]      [,2]      [,3]      [,4]      [,5]

[1,] 1.0000000 0.9648541 0.9649637 0.9629956 0.9599406

[2,] 0.9648541 1.0000000 0.9947731 0.9933266 0.9907213

[3,] 0.9649637 0.9947731 1.0000000 0.9943725 0.9916071

[4,] 0.9629956 0.9933266 0.9943725 1.0000000 0.9912263

[5,] 0.9599406 0.9907213 0.9916071 0.9912263 1.0000000

> cor(e$beta_age)

          [,1]      [,2]      [,3]      [,4]

[1,] 1.0000000 0.9897810 0.9909723 0.9873223

[2,] 0.9897810 1.0000000 0.9947695 0.9911093

[3,] 0.9909723 0.9947695 1.0000000 0.9942685

[4,] 0.9873223 0.9911093 0.9942685 1.0000000

> cor(e$beta_race)

          [,1]      [,2]      [,3]      [,4]

[1,] 1.0000000 0.9986930 0.9992706 0.9992229

[2,] 0.9986930 1.0000000 0.9980335 0.9980242

[3,] 0.9992706 0.9980335 1.0000000 0.9985834

[4,] 0.9992229 0.9980242 0.9985834 1.0000000

--

Jiqiang 

On Wed, Jan 9, 2013 at 11:44 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

[Bob Carpenter]

On 1/9/13 11:44 PM, Andrew Gelman wrote:
> J
>

> Actually, even before the boundary-avoiding priors, you can try running with the variance parameters set to fixed (and reasonable values). If _that_ is slow, then we have a different sort of computational problem. If that goes fast, then we can try gradually relaxing the priors until it starts to go wrong.

To echo Andrew, we should put this advice in the manual!
It's what Andrew's always said and what I've always tried
in estimating these models myself.

> I disagree with Bob's suggestion that we should be increasing the max number of leapfrog steps.

I figured you would :-)

I know you wanted to set the default even lower.

I should've clarified that increasing the number of leapfrog
steps may help with effective sample size per second for
a fixed model. A better strategy may be looking for a different model
or different parameterization of the same model, as Andrew
suggests above (and below).

> My intuition is that, if we're taking 1000 steps on this simple hierarchical regression, we need to figure out what's going wrong, rather than trying to solve it using brute force.

I agree. Looks like it may be the crazy-high
correlation among the coefficients for different
levels.

- Bob

[Andrew Gelman]

Yah, that's weird.

--
 
 

[Yuanjun Gao]

My NDNUTS(non-diagonal mass NUTS) manages to run 2000 iterations in about 17 minutes. After the warm-up adaptation, the leap-frog per iteration is generally 3,4,5 (I mean, 2^3, 2^4, 2^5).

To better compare the result, I am trying to run the default stan, which only manages to run 300 iterations after 30 minutes. So no doubt my sampler is much faster in this case.

However, in my two runs, one of them is working perfectly with every Rhat smaller than 1.1 while the other does not mixing well(with some Rhat close to 2). Surprisingly, the poor mixing happens to parameters well the correlation is not significant. Like z_state.

I feel that the problem of my sampler is that I am not doing enough shrinkage, so that the limited number of warm-up iterations create some "false" correlation. But I don't have a very good intuition on how much shrinkage I should do and how I should split the warm-up phase. I am afraid that by doing more shrinkage, parameters which do have strong correlation will have a poorer performance.

I was wondering how did you come up with the estimation procedure for diagonal NUTS, especially how did you decide to split the warm-up phase into 10, 20, 40, 80, 160 ,...? Is there some theory or just some empirical result?

Also I am curious about the fact that there can be some parameters with n_eff to be 4000. How can a sampler work so well and how to stan calculate n_eff?

Thank you!

[Yuanjun Gao]

I also feel that the non-identifiability is one of the main reasons for running slow. In this model. the intercept a, and the beta's are non-identifiable, which results in the high correlation between the beta's in the same category, as is shown by Jiqiang.

When I was doing the BDA assignments, I will usually drop one of beta's in order to get an identifiable model. For example, I will drop beta_sex for male and only use beta_sex for female. I do know that by doing this I introduce some kind of asymmetry, but at least it runs fast...

If this kind of non-identifiability is desirable in order to preserve the symmetric form, I definitely believe that a good estimation of cov matrix will improve NUTS greatly. Also RMHMC may be powerful in this case.

[Andrew Gelman]

Y

e should discuss the nonidentifiabiltiy. In many cases, nonidentifiability does _not_ slow down the program, all you have to do is monitor convergence on various derived quantities. See page 465 of my book with Jennifer for some discussion of this point. But it sounds like you're saying that, in this example, it's actually causing Nuts to take more steps?

In any case, we really don't need to monitor R.hat on the coefficients, only on the derived quantities.

For example, suppose you have a model, y[i] ~ N (mu + tau*alpha[j[i]], sigma), with 3 groups, thus three parameters alpha[j] ~ N(0,1). In this case, no need to monitor convergence of mu and alpha. Instead, monitor convergence of mu_adj and alpha_adj, where mu_adj = mu + mean(alpha), and alpha_adj[j] = alpha[j] - mean(alpha) for each j. To be clear: mean(alpha) is _not_ the posterior mean. It is the average of the alpha's, that is, (alpha[1]+alpha[2]+alpha[3])/3.

A

> --
>

[Yuanjun Gao]

Dear Professor Gelman,

Sorry for the late response, yes I do feel that the non-identifiability problem is making NUTS slow. The high-correlation will usually force NUTS to take more steps. I am not sure if this is the case for all the HMC method or is it just for NUTS.

Best,
Yuanjun