Initializing a hierarchical probit

[Ryan Dew]

Hi all,

Thank you for building such a great tool. I’ve been using rstan with great success over the past year, and until now have been able to sort out my problems on my own. However, I’m now having an issue with something that seems simple, but I can’t figure out what’s going on.

The model is a hierarchical probit model, with simulated data of individuals purchasing a given product. Individuals are simulated to have 20 observations each, where the person either buys or does not buy at each occassion, and the covariates are the price and an indicator for coupon, all set randomly. The simulation is included in the attached code.

I am having problems with the initialization: regardless of whether I set the inits, or let stan do it, I get the same error:

“SAMPLING FOR MODEL 'probit' NOW (CHAIN 1).
Error :

Initialization between (-2, 2) failed after 100 attempts.
 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

error occurred during calling the sampler; sampling not done”

I have tried initializing at exactly the simulated values and get the error. There are a few other things I’ve investigated.

1. The non-hierarchical equivalent has no issues, with the parameter set at the prior mean from the hierarchical version
   
2. The same error presents regardless of the hyperprior specification. I have tried not specifying a hyper prior, set it to be Wishart, Inverse Wishart, and as you can see in the attached, the recommended LKJ formulation from the stan reference manual. Regardless, I get this error.
   
3. I have tried specifying the model in a vectorized form (beta and y vectorized), and nothing changes.
   
4. I’ve tried varying the simulation parameters and the initialization, and nothing changes the error.

I believe I’ve set the constraints correctly, and I don’t think there is anything noteworthy (bad or good) about the parametrization, but undoubtedly I'm wrong on one of these. Please let me know if you have any insight!

Thanks for your help!

Ryan


In case needed, the sessionInfo():

R version 3.1.3 (2015-03-09)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.1 (Yosemite)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods 
[7] base    

other attached packages:
[1] rstan_2.6.0   inline_0.3.13 Rcpp_0.11.5   MASS_7.3-39 

loaded via a namespace (and not attached):
[1] codetools_0.2-10 stats4_3.1.3     tools_3.1.3   

[Bob Carpenter]

I don't see any immediate problem with the constraints.
How big is K? There may be instability with the correlation
matrix or multiplication.

In general, you're better off using a Cholesky factor
parameterization instead of a full-on covariance matrix.

There's an example in the manual, but the basic idea is
that you are much better for for arithmetic stability and
efficiency with:

cholesky_corr[K] L_Omega;
...

matrix[K,K] L_Sigma;
L_Sigma <- diag_pre_multiply(tau, L_Omega);

L_Omega ~ lkj_corr_cholesky(2);


You can then reconstruct Omega as L_Omega * L_Omega'
in the generated quantities if you want to look at it.

- 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.
> <probit_model.txt>

[Ryan Dew]

Hi Bob,

Thanks for the help. This worked perfectly. K was small—only 3—so I didn’t think it would make a big difference to parametrize it like that, but it seems it did!

Best,

Ryan

--
Ryan Dew
RDe...@gsb.columbia.edu

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/quzreoGxNcM/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

I think we'll be making it more stable in 2.7. Ben
found an issue with how we were validating that we have
fixed on a branch. It should still be much more efficient
with the Cholesky version.

- Bob

[Ryan Dew]

I also want to add that I discovered far too late that, after my initial trials with inits and no inits, I had forgotten to re-include the init part of the stan() command, as evident in the code I posted at the start. So for some of my initial value trials, it wasn't actually calling the specified inits. When I legitimately initialize the posted model at the true parameter values, it does work, and also at random values near zero, although the cholesky factorization version is much more stable (at a wide range of starting values) and more efficient. I discovered this only after I had the cholesky version up and running, and it's a pretty embarrassing mistake, but since I've referenced more than a handful of these threads looking for solutions to problems, I wanted to include this endnote in case anyone else has this problem.

Thanks again for the help.