Re: [stan-users] Digest for stan-users@googlegroups.com - 15 Messages in 4 Topics
85 views
Skip to first unread message
Alexander Budzier
Sep 10, 2012, 1:46:44 PM9/10/12
to stan-...@googlegroups.com
On 9/10/12, stan-...@googlegroups.com <stan-...@googlegroups.com> wrote:
> =============================================================================
> Today's Topic Summary
> =============================================================================
>
> Group: stan-...@googlegroups.com
> Url: http://groups.google.com/group/stan-users/topics
>
> - random effects models with known covariance matrix [2 Updates]
> http://groups.google.com/group/stan-users/t/1dcbd8683bbbd7f6
> - IRT models PL1, PL2 + multilevel versions [4 Updates]
> http://groups.google.com/group/stan-users/t/53f175a7b6128aba
> - optimization level in RStan & discrete parameters in Stan [5 Updates]
> http://groups.google.com/group/stan-users/t/3884a2491bed2721
> - installing Stan [4 Updates]
> http://groups.google.com/group/stan-users/t/ca40e9dada65ed5e
>
>
> =============================================================================
> Topic: random effects models with known covariance matrix
> Url: http://groups.google.com/group/stan-users/t/1dcbd8683bbbd7f6
> =============================================================================
>
> ---------- 1 of 2 ----------
> From: "f....@gmx.net" <f....@gmx.net>
> Date: Sep 10 09:22AM -0700
> Url: http://groups.google.com/group/stan-users/msg/e1c1c80b85cfc684
>
> Hello,
>
> I wanted to try to use STAN for fitting more complex random effects models,
>
> where the random effects have covariance matrices that are known up to a
> constant.
>
> I tried to fit the basic "animal model" (used for computing breeding values
>
> for selection purposes) to get started. In this model, each individual is
> assigned a random effect, which has the following prior
>
>
> data {
> matrix[N,N] A ;
> ...
> }
>
> parameters {
> vector[N] u;
> real <lower=0> sigmasq_u;
> ...
> }
>
> ...
> model{
> u ~ multi_norm(0 , A * sigmasq_u)
> }
>
> where A is the "numerator relationship matrix" and is known from pedigree,
> theory or from genomic marker data.
>
> Even for a relatively moderate size of N = 500 individuals, does STAN take
> extremely long to get in just 500 iterations, and for N = 1000, STAN will
> use up all of my 16 GB memory.
> This can't be, even if STAN is inverting A over and over again.
>
>
> The code is attached. Any suggestion would be appreciated.
>
>
> Regards,
>
> Frank
>
>
> ---------- 2 of 2 ----------
> From: Marcus Brubaker <marcus.a...@gmail.com>
> Date: Sep 10 12:46PM -0400
> Url: http://groups.google.com/group/stan-users/msg/a16374a3ae894d71
>
> One thing you can try is to use multi_normal_cholesky
>
> I.E., add
>
> transformed data {
> matrix[N,N] A_chol;
> A_chol <- cholesky_decompose(A);
> }
>
> then in the model
>
> u ~ multi_normal_cholesky(0,sqrt(sigmasq_u)*A_chol);
>
> Also, you should be able to use the vectorized version of the normal
> distribution. I.E., instead of the for loop over ys write
>
> y ~ normal(beta0 + u, sqrt(sig2_e));
>
> Cheers,
> Marcus
>
>> Hello,
>
>> I wanted to try to use STAN for fitting more complex random effects
> models,
>> where the random effects have covariance matrices that are known up to a
>> constant.
>
>> I tried to fit the basic "animal model" (used for computing breeding
> values
>
>
>
> =============================================================================
> Topic: IRT models PL1, PL2 + multilevel versions
> Url: http://groups.google.com/group/stan-users/t/53f175a7b6128aba
> =============================================================================
>
> ---------- 1 of 4 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 10 12:05AM -0400
> Url: http://groups.google.com/group/stan-users/msg/d130f525a1b325d7
>
> Item-response theory (IRT) models are very much the kind
> of thing we had in mind when designing Stan -- it's essentially
> a logistic regression with high correlation among
> the parameters in the posterior.
>
> I just wrote four Stan models for the basic IRT
> models PL1 (aka Rasch) with no question discrimination
> parameter and PL2 with the discrimination parameter, along
> with multilevel variants of the kind described in Gelman
> and Hill's regression book starting on p. 314.
>
> You can browse the models online in our source repository:
>
> https://code.google.com/p/stan/source/browse/#git%2Fsrc%2Fmodels%2Fmisc%2Firt
>
> I tested them by generating data for 400 students
> on 100 questions, with 25% of the data missing at
> random. This results in around 30,000 data points
> and 500 (PL1) or 600 (PL2) parameters.
>
> The models take about 5 seconds to compile in g++ at -O3
> optimization.
>
> Here's how long Stan took to warmup for 500 iterations
> then sample for 500 iterations. All models converged
> to appropriate answers across multiple chains.
>
> PL1: 80 seconds (1m 20s)
>
> PL2: 180 seconds (3m 00s)
>
> The hierarchical models took the same amount of time
> as their non-hierarchical counterparts.
>
> The PL2 models adapt to a smaller step size and more
> steps, which is why they take longer (that and having
> 20% more parameters and their attendant differentiation
> and updating).
>
> It only consumes 14MB of memory peak, so there's no
> problem running multiple processes concurrently.
>
> I also ran on much more data. Basic IRT for 1000 students
> with 200 questions each (150K data points) ran
> 1000 warmup and 1000 sampling iterations (twice as many
> as reported above) in 4 minutes in 70MB and converged to
> the right answer. I even ran basic IRT for 10,000 students with 200
> questions each (1.5M data points), and though I lost
> patience testing, I extrapolated about 2 hours for
> 1000 warmup and 1000 sampling iterations. (I cut down
> the later tests when I realized it was converging very
> quickly from random inits.)
>
> There's a bit of variation in the runs mainly due to
> waiting for the step sizes to converge during warmup;
> if the random initialization is far from a mode, especially
> for the hierarchical parameters or intercept terms, then
> it can take a bit longer to move to the mode because of
> the need for reduced step sizes to keep the simulation
> accurate (I never saw more than about 50% more time needed
> than reported above). With good inits it can run a bit faster than
> reported above.
>
> The notebook I tested this on is reasonably beefy:
>
> Mac OS X, Mountain Lion,
> 2.3 GHz i7 CPU
> 1600MHz DDR3 memory
>
> Although I'm using the latest g++ from Xcode, version 4.21,
> the latest from Macports is faster (as much as 20% in
> some of our earlier tests of a version of 4.6 and they
> continue to get faster).
>
> - Bob
>
>
> ---------- 2 of 4 ----------
> From: Emmanuel Charpentier <emanuel.c...@gmail.com>
> Date: Sep 10 01:12AM -0700
> Url: http://groups.google.com/group/stan-users/msg/2f1abac00f4419f5
>
> Dear Bob,
>
> Thank you for this effort. I have skimmed the stan files, and notice a
> couple of things :
>
> 1) You use one more parameter, the mean ability, whih, from a "sampling
> efficiency" point of view, should play the same role that the mean of the
> alpha distribution in Gelman & Hill's textbook.
>
> 2) You give to this mean an informative prior. This would be hard to defend
>
> against a reviewer.
>
> 3) Similarly, you give abilities, discriminations and difficulties somewhat
>
> informative priors. The same remark applies : whereas imposing a standard
> deviation on abilities is just a choice of scale; imposing a small standard
>
> deviation to discriminations and difficulties *is* informative. Hard to do
> when you are exploring a situation where no prior knowledge has been
> formalized.
>
> 4) However, in the IRT2 hierrchical model, you do not state a prior for the
>
> abilities and difficulties. This, as far as I know, is equivalent to giving
>
> them the (improper) prior on (0, \infty), ths making the posterior
> improper, AFAIK. So much for model comparison...
>
> In short, where Gelman & Hill's solution introduced both an additive and a
> multiplicative redundancy (by means of centering and scaling "raw"
> parameters), you keep the additive redundancy and replace the
> multiplicative redundancy by imposing more informative distributions.
>
> I'll try to assess the impact of these choices by trying various values of
> the prior parameters.
>
> Again, thank you !
>
> Emmanuel
> Charpentier
>
>
> ---------- 3 of 4 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 10 11:58AM -0400
> Url: http://groups.google.com/group/stan-users/msg/601f3a3e3a3b8eb4
>
> You can put whatever priors you want on the
> variables. It doesn't affect the computation much.
> It does affect the concentration of the posterior,
> which is very sensitive to the choice of prior in
> these models.
>
> I also didn't provide a prior on the deviation
> parameters, which I should probably do.
>
> On 9/10/12 4:12 AM, Emmanuel Charpentier wrote:
>
>> Thank you for this effort. I have skimmed the stan files, and notice a
>> couple of things :
>
>> 1) You use one more parameter, the mean ability, whih, from a "sampling
>> efficiency" point of view, should play the same
>> role that the mean of the alpha distribution in Gelman & Hill's textbook.
>
> Exactly. You can parameterize it the other
> way, but it's then adding a dependency to all the
> parameters on the mean, which is less flexible than
> just having the mean float.
>
>> 2) You give to this mean an informative prior. This would be hard to
>> defend against a reviewer.
>
> Presumably in the same sense that any informative prior might
> be hard to defend against a reviewer!
>
> I could've added the priors as data as in the S. McKay Curtis
> JSS paper on BUGS item-response models, so that they'd be
> easier for users to play with.
>
> I found when the priors were fatter than the values used
> to simulate the data, the mean posterior estimates drifted
> from their simulated values (in predictable, repeatable ways).
>
> I'd recommend the hierarchical models, which don't have this
> problem, but can also be difficult to defend against reviewers
> in my experience.
>
>> : whereas imposing a standard deviation on abilities is just a choice of
>> scale; imposing a small standard deviation to
>> discriminations and difficulties *is* informative. Hard to do when you are
>> exploring a situation where no prior
>> knowledge has been formalized.
>
> Indeed. As I said above, feel free to change the priors!
>
>> 4) However, in the IRT2 hierrchical model, you do not state a prior for
>> the abilities and difficulties. This, as far as
>> I know, is equivalent to giving them the (improper) prior on (0, \infty),
>> ths making the posterior improper, AFAIK. So
>> much for model comparison...
>
> Exactly -- the scale parameters have improper priors.
> The posterior doesn't behave improperly, but I haven't done
> a formal derivation.
>
> Given how the model behaves without a prior (recovering
> the true simulated values), and given the quantity of
> data for the priors (# of questions, # of students), a
> very weak prior shouldn't have much effect.
>
> It's easy enough to add a prior. Andrew likes using
> weakly informative priors in these situations, but
> I realize other researchers have different approaches.
>
>> In short, where Gelman & Hill's solution introduced both an additive and a
>> multiplicative redundancy (by means of
>> centering and scaling "raw" parameters), you keep the additive redundancy
>> and replace the multiplicative redundancy by
>> imposing more informative distributions.
>
> Gelman and Hill didn't propose a single solution. They discuss
> several methods (setting one value based on prior knowledge,
> setting the mean to 0, then their actual model also sets the
> scale to 1 by using a z-score type normalization).
>
> I was just trying to use the simplest thing possible.
>
>> I'll try to assess the impact of these choices by trying various values of
>> the prior parameters.
>
> Let us know what happens. I don't have that much time to play
> around with IRT parameterizations.
>
> - Bob
>
>
> ---------- 4 of 4 ----------
> From: Emmanuel Charpentier <emanuel.c...@gmail.com>
> Date: Sep 10 09:37AM -0700
> Url: http://groups.google.com/group/stan-users/msg/eb6b950c353f97c0
>
> A couple comments below :
>
> Le lundi 10 septembre 2012 17:58:14 UTC+2, Bob Carpenter a écrit :
>> to simulate the data, the mean posterior estimates drifted
>> from their simulated values (in predictable, repeatable ways).
>
>> This is a bit of a problem, don't you think ? You are analyzing these data
>>
> to get information about unknown parameters from your data. What you state
> means that your process will "drift" (i. e. lose precision) if you do not
> give them a bit pf prior information (obtained by sheer divination in
> this case ?). A bit too circular for some reviewers...
>
> I'd recommend the hierarchical models, which don't have this
>> problem, but can also be difficult to defend against reviewers
>> in my experience.
>
> I went for the hierarchical solution because of this. but your hierarchical
>
> solution also incorporates information in the form of informative
> hyperpriors.
>
>
>> Exactly -- the scale parameters have improper priors.
>> The posterior doesn't behave improperly, but I haven't done
>> a formal derivation.
>
> Since one of my goals is odel comparison, I'd be wary of anything leading
> to improper posteriors...
>
>> the true simulated values), and given the quantity of
>> data for the priors (# of questions, # of students), a
>> very weak prior shouldn't have much effect.
>
> I'd be in serious trouble if I had to *prove* that the resultant
> simulations do not come from an improper posterior. I tried this derivation
>
> and found it a bit too hard for my limited algebraic/analytic abilities
> (I'm definitely better at surgery...). A bit of bibliography led me to the
> recent Jean-Paul Fox's monograph, which did not let me with a useable way
> of simulating from this difficult posterior.
>
>
>> It's easy enough to add a prior. Andrew likes using
>> weakly informative priors in these situations,
>
>
> So do I, because 1) your priors being weak, you can succesfully argue that
> they do not carry much information into the posterior, and 2) being proper,
>
> you can (try to) use their resultant simulations in marginal likelihood
> estimation / Bayes factors / whatever model comparison method du jour you
> wish to use.
>
>
>> > imposing more informative distributions.
>
>> Gelman and Hill didn't propose a single solution. They discuss
>> several methods (setting one value based on prior knowledge,
>
> That's for another problem, i. e. to resolve the sign indetermination of
> the general model. Which does not arise in my case, since I can use subject
>
> matter knowledge to postulate that all my items are correlated the same
> way. Hence the lognormal distribution of the gammas.
>
>
>> of the prior parameters.
>
>> Let us know what happens. I don't have that much time to play
>> around with IRT parameterizations.
>
> Neither do I, alas, but I just have to : "When a researche gotta, he's
> gonna...". I was hoping that Stan's superior sampler could help me with a
> better simulation without having to wait hours (literally...) for
> convergence. But the lack of sampling from discrete distributions (if only
> for the missing data) makes it (currently) unfit for my current purposes.
>
> I will definitely try to compare JAGS and Stan (and possibly R Metropolis
> adaptive samplers such as the one exemplified in the adaptMCMC package,
> which performs surprisingly well with very little tuning) for the
> hierarchical models in the complete-cases case. But not now...
>
> Again, thank you for your attention and your time.
>
> Emmanuel
> Charpentier
>
>
>
> =============================================================================
> Topic: optimization level in RStan & discrete parameters in Stan
> Url: http://groups.google.com/group/stan-users/t/3884a2491bed2721
> =============================================================================
>
> ---------- 1 of 5 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 09 10:56PM -0400
> Url: http://groups.google.com/group/stan-users/msg/b4e8120d25ca54f0
>
> On 9/9/12 10:41 AM, Emmanuel Charpentier wrote:
>> the optimization level
>>
>> <http://code.google.com/p/stan/wiki/RStanGettingStarted#Modify_C++_compiler_optimization_level>
>> a hack.
>
>> OK. I'll look into changing that cleanly (maybe on a per-sesion basis,
>> through R's options()).
>
> I don't think that's possible. If you do figure out a way
> to do it, could you let us know?
>
> As a side benefit, -O3 will make other R apps compiled
> with R CMD SHLIB faster :-)
>
> Here's a thread on stackoverflow on the topic:
>
> http://stackoverflow.com/questions/5789673/r-c-optimization-flag-when-using-the-inline-package
>
>> Aha. I overlooked that in the manual. So currently Stan can't impute
>> *discrete* missing data. Nor assess the
>> distributions of *discrete* parameters (so much for model choice
>> algorithms such as umbrella sampling or hierarchical
>> model choice (à la Carlin & Chib (1995) or otherwise)...).
>
> Right. Unless you can marginalize out the
> discrete parameter, Stan can't handle it.
>
> I'm about to put a compile-time check to flag attempted
> uses of discrete params in the parser, which I forgot to
> do the first time around.
>
>> Any idea of the ETA for a discrete sampler ?
>
> We'll probably write a simple Gibbs for finite
> and slice for infinite sampler. Probably some
> time in the next six or nine months.
>
> - Bob
>
>
> ---------- 2 of 5 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 09 11:11PM -0400
> Url: http://groups.google.com/group/stan-users/msg/a934d4fa365f6c63
>
> n 9/9/12 5:26 AM, Emmanuel Charpentier wrote:
>> never initialized.
>
>> OK. Just a little problem : THERE IS NO theta VARIABLE IN THE MODELS. I
>> suspect that the compiler is goofing in building
>> it's error message.
>
> I see the confusion here.
>
> We should make that error message clearer.
> That's just a generic report from binomial_log, where
> the probability-of-success parameter is called theta.
> We should just get rid of the names altogether and refer
> to the names of the variable, like chance-of-success
> or scale or location.
>
> The problem for us is at the point that warning (and it's
> just a warning) gets reported, we've lost the link to the
> underlying variable in the Stan model. It would be possible
> to fix this by passing arguments representing the names of
> all the variables to the distributions or by trying to catch
> errors in the Stan code, but both are a lot of work and
> would make the code very ugly. Probably worth it in the
> longer term, as this is going to be very frustrating for
> users.
>
> - Bob
>
>
> ---------- 3 of 5 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 09 11:15PM -0400
> Url: http://groups.google.com/group/stan-users/msg/86a71d3d880456cc
>
> On 9/9/12 5:40 AM, Emmanuel Charpentier wrote:
>> initialization (to 0, rather thn sampling from the prior as in OpenBUGS)
>> seems to me more reasonable in many cases. I
>> love the stan() R function range of options : "0", "random" or
>> user-specified values. The best of ll worlds would be a
>> named list of parameters, and a default, which could be either "0" or
>> "random". But that might be treated at the R level.
>
> That's what we're aiming for. It'll take a bit of
> work and there are other fires to put out at the
> moment on basic error checking that we're going to
> address first.
>
> The biggest hurdle for the implementation is that the
> inits are on the unconstrained scale.
>
> I don't think JAGS always initializes to 0. It says
> in the manual that it takes reasonable values.
> People like to use very fat priors even when they
> have a good sense of what inits should be. Maybe
> Andrew has something to say about that -- he likes
> dispersed initializations.
>
> RStan and Stan use the same basic C++ code, but
> they have different command-line controls for marshaling
> the arguments. The C++ needs the data to be given
> to the model's class constructor and the constructor
> computes the transformed data, all of which then remains
> fixed throughout sampling.
>
> So far, we haven't had problems with Stan converging,
> even from pretty far-out inits. If they do get so
> bad that they result in NaN or infinite values, we
> resample.
>
> - Bob
>
>
> ---------- 4 of 5 ----------
> From: Andrew Gelman <gel...@stat.columbia.edu>
> Date: Sep 09 11:19PM -0400
> Url: http://groups.google.com/group/stan-users/msg/11ea72e6ff69d20
>
> On Sep 9, 2012, at 11:15 PM, Bob Carpenter wrote:
>
>> have a good sense of what inits should be. Maybe
>> Andrew has something to say about that -- he likes
>> dispersed initializations.
>
> If you're doing Gibbs, it's a good idea to avoid deterministic inits. For
> example, if you set theta_1,…,theta_J=0, then you might estimate tau=0, etc.
> This is probably less of an issue with Nuts or HMC. Still, I'd generally
> prefer some randomness in the inits (except when we want to do some pure
> comparison as in our speed comparisons with Jags).
>
>
> ---------- 5 of 5 ----------
> From: Emmanuel Charpentier <emanuel.c...@gmail.com>
> Date: Sep 10 01:18AM -0700
> Url: http://groups.google.com/group/stan-users/msg/56b72e4eab6acd8f
>
> Dear Bob,
>
> Le lundi 10 septembre 2012 05:11:30 UTC+2, Bob Carpenter a écrit :
>
>
>> longer term, as this is going to be very frustrating for
>> users.
>
>> Indeed...
>
> Re-reading my code (both R and Stan), I have been unable to identify the
> uninitialized parameters. I'm probably becoming myopic with this code, so
> I'll give it a rest, and try to debug it a few days later.
>
> Thank you again !
>
>
> Emmanuel Charpentier
>
>
>
> =============================================================================
> Topic: installing Stan
> Url: http://groups.google.com/group/stan-users/t/ca40e9dada65ed5e
> =============================================================================
>
> ---------- 1 of 4 ----------
> From: Shravan Vasishth <vasishth...@gmail.com>
> Date: Sep 09 11:15AM -0700
> Url: http://groups.google.com/group/stan-users/msg/85cfd4d792048baf
>
> *Sigh*. I managed to live my life out successfully without ever having to
> write (or even look inside, with some exceptions) a Makefile. Advice taken;
>
> I will get to know Makefiles, if only to avoid making a fool of myself in
> the future.
>
> On Sunday, September 9, 2012 12:07:58 PM UTC+2, Emmanuel Charpentier wrote:
>
>
> ---------- 2 of 4 ----------
> From: Daniel Lee <bea...@alum.mit.edu>
> Date: Sep 09 03:30PM -0400
> Url: http://groups.google.com/group/stan-users/msg/ba63641e7449c99b
>
> Shravan,
>
> Before you get too frustrated, all you have to do is to first change
> directory into the Stan directory, then type at the command line "make
> bin/libstan.a". The static library will be built for you at
> bin/libstan.a.
>
> There really isn't a need to get familiar with make, we (the Stan dev
> team) will try to encapsulate as much of the build pain through tools
> such as make.
>
>
>
> Daniel
>
>
>
>
> On Sun, Sep 9, 2012 at 2:15 PM, Shravan Vasishth
>
>
> ---------- 3 of 4 ----------
> From: Shravan Vasishth <vasishth...@gmail.com>
> Date: Sep 09 12:38PM -0700
> Url: http://groups.google.com/group/stan-users/msg/9f7d446c004f7797
>
> Hi Daniel,
>
> yes, I did follow the instructions without thinking about it the second
> time round! Everything works as expected.
>
> Thanks, and apologies for the pointless thread.
>
> On Sunday, September 9, 2012 9:30:30 PM UTC+2, Daniel Lee wrote:
>
>
> ---------- 4 of 4 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 09 10:47PM -0400
> Url: http://groups.google.com/group/stan-users/msg/51b95586f4deb219
>
> On 9/9/12 3:38 PM, Shravan Vasishth wrote:
>> Hi Daniel,
>
>> yes, I did follow the instructions without thinking about it the second
>> time round! Everything works as expected.
>
>> Thanks, and apologies for the pointless thread.
>
> Please don't feel it's your fault for not puzzling
> out how to run Stan from the command line. This is our fault
> for distributing software that doesn't install and
> run like a regular desktop app. (RStan does run
> like a regular R app.)
>
> We're trying to make it as easy as possible
> to use Stan, so this kind of question and feedback
> is really useful to us. We'll work on making the usage
> instructions clearer. Or ideally, figure out
> how to get it run relocatably (any help on
> the part of users who know how to do this better
> than us would also be greatly appreciated).
>
> In the meantime, please send us questions before
> you get too frustrated. We're more than happy to
> help.
>
> - Bob
>
>
>
>
>
--
Sent from my mobile device
------------------------------
Alexander Budzier
+1 (440) 253-9266
+44 (750) 210-7145
+49 (211) 136-4192
> =============================================================================
> Today's Topic Summary
> =============================================================================
>
> Group: stan-...@googlegroups.com
> Url: http://groups.google.com/group/stan-users/topics
>
> - random effects models with known covariance matrix [2 Updates]
> http://groups.google.com/group/stan-users/t/1dcbd8683bbbd7f6
> - IRT models PL1, PL2 + multilevel versions [4 Updates]
> http://groups.google.com/group/stan-users/t/53f175a7b6128aba
> - optimization level in RStan & discrete parameters in Stan [5 Updates]
> http://groups.google.com/group/stan-users/t/3884a2491bed2721
> - installing Stan [4 Updates]
> http://groups.google.com/group/stan-users/t/ca40e9dada65ed5e
>
>
> =============================================================================
> Topic: random effects models with known covariance matrix
> Url: http://groups.google.com/group/stan-users/t/1dcbd8683bbbd7f6
> =============================================================================
>
> ---------- 1 of 2 ----------
> From: "f....@gmx.net" <f....@gmx.net>
> Date: Sep 10 09:22AM -0700
> Url: http://groups.google.com/group/stan-users/msg/e1c1c80b85cfc684
>
> Hello,
>
> I wanted to try to use STAN for fitting more complex random effects models,
>
> where the random effects have covariance matrices that are known up to a
> constant.
>
> I tried to fit the basic "animal model" (used for computing breeding values
>
> for selection purposes) to get started. In this model, each individual is
> assigned a random effect, which has the following prior
>
>
> data {
> matrix[N,N] A ;
> ...
> }
>
> parameters {
> vector[N] u;
> real <lower=0> sigmasq_u;
> ...
> }
>
> ...
> model{
> u ~ multi_norm(0 , A * sigmasq_u)
> }
>
> where A is the "numerator relationship matrix" and is known from pedigree,
> theory or from genomic marker data.
>
> Even for a relatively moderate size of N = 500 individuals, does STAN take
> extremely long to get in just 500 iterations, and for N = 1000, STAN will
> use up all of my 16 GB memory.
> This can't be, even if STAN is inverting A over and over again.
>
>
> The code is attached. Any suggestion would be appreciated.
>
>
> Regards,
>
> Frank
>
>
> ---------- 2 of 2 ----------
> From: Marcus Brubaker <marcus.a...@gmail.com>
> Date: Sep 10 12:46PM -0400
> Url: http://groups.google.com/group/stan-users/msg/a16374a3ae894d71
>
> One thing you can try is to use multi_normal_cholesky
>
> I.E., add
>
> transformed data {
> matrix[N,N] A_chol;
> A_chol <- cholesky_decompose(A);
> }
>
> then in the model
>
> u ~ multi_normal_cholesky(0,sqrt(sigmasq_u)*A_chol);
>
> Also, you should be able to use the vectorized version of the normal
> distribution. I.E., instead of the for loop over ys write
>
> y ~ normal(beta0 + u, sqrt(sig2_e));
>
> Cheers,
> Marcus
>
>> Hello,
>
>> I wanted to try to use STAN for fitting more complex random effects
> models,
>> where the random effects have covariance matrices that are known up to a
>> constant.
>
>> I tried to fit the basic "animal model" (used for computing breeding
> values
>
>
>
> =============================================================================
> Topic: IRT models PL1, PL2 + multilevel versions
> Url: http://groups.google.com/group/stan-users/t/53f175a7b6128aba
> =============================================================================
>
> ---------- 1 of 4 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 10 12:05AM -0400
> Url: http://groups.google.com/group/stan-users/msg/d130f525a1b325d7
>
> Item-response theory (IRT) models are very much the kind
> of thing we had in mind when designing Stan -- it's essentially
> a logistic regression with high correlation among
> the parameters in the posterior.
>
> I just wrote four Stan models for the basic IRT
> models PL1 (aka Rasch) with no question discrimination
> parameter and PL2 with the discrimination parameter, along
> with multilevel variants of the kind described in Gelman
> and Hill's regression book starting on p. 314.
>
> You can browse the models online in our source repository:
>
> https://code.google.com/p/stan/source/browse/#git%2Fsrc%2Fmodels%2Fmisc%2Firt
>
> I tested them by generating data for 400 students
> on 100 questions, with 25% of the data missing at
> random. This results in around 30,000 data points
> and 500 (PL1) or 600 (PL2) parameters.
>
> The models take about 5 seconds to compile in g++ at -O3
> optimization.
>
> Here's how long Stan took to warmup for 500 iterations
> then sample for 500 iterations. All models converged
> to appropriate answers across multiple chains.
>
> PL1: 80 seconds (1m 20s)
>
> PL2: 180 seconds (3m 00s)
>
> The hierarchical models took the same amount of time
> as their non-hierarchical counterparts.
>
> The PL2 models adapt to a smaller step size and more
> steps, which is why they take longer (that and having
> 20% more parameters and their attendant differentiation
> and updating).
>
> It only consumes 14MB of memory peak, so there's no
> problem running multiple processes concurrently.
>
> I also ran on much more data. Basic IRT for 1000 students
> with 200 questions each (150K data points) ran
> 1000 warmup and 1000 sampling iterations (twice as many
> as reported above) in 4 minutes in 70MB and converged to
> the right answer. I even ran basic IRT for 10,000 students with 200
> questions each (1.5M data points), and though I lost
> patience testing, I extrapolated about 2 hours for
> 1000 warmup and 1000 sampling iterations. (I cut down
> the later tests when I realized it was converging very
> quickly from random inits.)
>
> There's a bit of variation in the runs mainly due to
> waiting for the step sizes to converge during warmup;
> if the random initialization is far from a mode, especially
> for the hierarchical parameters or intercept terms, then
> it can take a bit longer to move to the mode because of
> the need for reduced step sizes to keep the simulation
> accurate (I never saw more than about 50% more time needed
> than reported above). With good inits it can run a bit faster than
> reported above.
>
> The notebook I tested this on is reasonably beefy:
>
> Mac OS X, Mountain Lion,
> 2.3 GHz i7 CPU
> 1600MHz DDR3 memory
>
> Although I'm using the latest g++ from Xcode, version 4.21,
> the latest from Macports is faster (as much as 20% in
> some of our earlier tests of a version of 4.6 and they
> continue to get faster).
>
> - Bob
>
>
> ---------- 2 of 4 ----------
> From: Emmanuel Charpentier <emanuel.c...@gmail.com>
> Date: Sep 10 01:12AM -0700
> Url: http://groups.google.com/group/stan-users/msg/2f1abac00f4419f5
>
> Dear Bob,
>
> Thank you for this effort. I have skimmed the stan files, and notice a
> couple of things :
>
> 1) You use one more parameter, the mean ability, whih, from a "sampling
> efficiency" point of view, should play the same role that the mean of the
> alpha distribution in Gelman & Hill's textbook.
>
> 2) You give to this mean an informative prior. This would be hard to defend
>
> against a reviewer.
>
> 3) Similarly, you give abilities, discriminations and difficulties somewhat
>
> informative priors. The same remark applies : whereas imposing a standard
> deviation on abilities is just a choice of scale; imposing a small standard
>
> deviation to discriminations and difficulties *is* informative. Hard to do
> when you are exploring a situation where no prior knowledge has been
> formalized.
>
> 4) However, in the IRT2 hierrchical model, you do not state a prior for the
>
> abilities and difficulties. This, as far as I know, is equivalent to giving
>
> them the (improper) prior on (0, \infty), ths making the posterior
> improper, AFAIK. So much for model comparison...
>
> In short, where Gelman & Hill's solution introduced both an additive and a
> multiplicative redundancy (by means of centering and scaling "raw"
> parameters), you keep the additive redundancy and replace the
> multiplicative redundancy by imposing more informative distributions.
>
> I'll try to assess the impact of these choices by trying various values of
> the prior parameters.
>
> Again, thank you !
>
> Emmanuel
> Charpentier
>
>
> ---------- 3 of 4 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 10 11:58AM -0400
> Url: http://groups.google.com/group/stan-users/msg/601f3a3e3a3b8eb4
>
> You can put whatever priors you want on the
> variables. It doesn't affect the computation much.
> It does affect the concentration of the posterior,
> which is very sensitive to the choice of prior in
> these models.
>
> I also didn't provide a prior on the deviation
> parameters, which I should probably do.
>
> On 9/10/12 4:12 AM, Emmanuel Charpentier wrote:
>
>> Thank you for this effort. I have skimmed the stan files, and notice a
>> couple of things :
>
>> 1) You use one more parameter, the mean ability, whih, from a "sampling
>> efficiency" point of view, should play the same
>> role that the mean of the alpha distribution in Gelman & Hill's textbook.
>
> Exactly. You can parameterize it the other
> way, but it's then adding a dependency to all the
> parameters on the mean, which is less flexible than
> just having the mean float.
>
>> 2) You give to this mean an informative prior. This would be hard to
>> defend against a reviewer.
>
> Presumably in the same sense that any informative prior might
> be hard to defend against a reviewer!
>
> I could've added the priors as data as in the S. McKay Curtis
> JSS paper on BUGS item-response models, so that they'd be
> easier for users to play with.
>
> I found when the priors were fatter than the values used
> to simulate the data, the mean posterior estimates drifted
> from their simulated values (in predictable, repeatable ways).
>
> I'd recommend the hierarchical models, which don't have this
> problem, but can also be difficult to defend against reviewers
> in my experience.
>
>> : whereas imposing a standard deviation on abilities is just a choice of
>> scale; imposing a small standard deviation to
>> discriminations and difficulties *is* informative. Hard to do when you are
>> exploring a situation where no prior
>> knowledge has been formalized.
>
> Indeed. As I said above, feel free to change the priors!
>
>> 4) However, in the IRT2 hierrchical model, you do not state a prior for
>> the abilities and difficulties. This, as far as
>> I know, is equivalent to giving them the (improper) prior on (0, \infty),
>> ths making the posterior improper, AFAIK. So
>> much for model comparison...
>
> Exactly -- the scale parameters have improper priors.
> The posterior doesn't behave improperly, but I haven't done
> a formal derivation.
>
> Given how the model behaves without a prior (recovering
> the true simulated values), and given the quantity of
> data for the priors (# of questions, # of students), a
> very weak prior shouldn't have much effect.
>
> It's easy enough to add a prior. Andrew likes using
> weakly informative priors in these situations, but
> I realize other researchers have different approaches.
>
>> In short, where Gelman & Hill's solution introduced both an additive and a
>> multiplicative redundancy (by means of
>> centering and scaling "raw" parameters), you keep the additive redundancy
>> and replace the multiplicative redundancy by
>> imposing more informative distributions.
>
> Gelman and Hill didn't propose a single solution. They discuss
> several methods (setting one value based on prior knowledge,
> setting the mean to 0, then their actual model also sets the
> scale to 1 by using a z-score type normalization).
>
> I was just trying to use the simplest thing possible.
>
>> I'll try to assess the impact of these choices by trying various values of
>> the prior parameters.
>
> Let us know what happens. I don't have that much time to play
> around with IRT parameterizations.
>
> - Bob
>
>
> ---------- 4 of 4 ----------
> From: Emmanuel Charpentier <emanuel.c...@gmail.com>
> Date: Sep 10 09:37AM -0700
> Url: http://groups.google.com/group/stan-users/msg/eb6b950c353f97c0
>
> A couple comments below :
>
> Le lundi 10 septembre 2012 17:58:14 UTC+2, Bob Carpenter a écrit :
>> to simulate the data, the mean posterior estimates drifted
>> from their simulated values (in predictable, repeatable ways).
>
>> This is a bit of a problem, don't you think ? You are analyzing these data
>>
> to get information about unknown parameters from your data. What you state
> means that your process will "drift" (i. e. lose precision) if you do not
> give them a bit pf prior information (obtained by sheer divination in
> this case ?). A bit too circular for some reviewers...
>
> I'd recommend the hierarchical models, which don't have this
>> problem, but can also be difficult to defend against reviewers
>> in my experience.
>
> I went for the hierarchical solution because of this. but your hierarchical
>
> solution also incorporates information in the form of informative
> hyperpriors.
>
>
>> Exactly -- the scale parameters have improper priors.
>> The posterior doesn't behave improperly, but I haven't done
>> a formal derivation.
>
> Since one of my goals is odel comparison, I'd be wary of anything leading
> to improper posteriors...
>
>> the true simulated values), and given the quantity of
>> data for the priors (# of questions, # of students), a
>> very weak prior shouldn't have much effect.
>
> I'd be in serious trouble if I had to *prove* that the resultant
> simulations do not come from an improper posterior. I tried this derivation
>
> and found it a bit too hard for my limited algebraic/analytic abilities
> (I'm definitely better at surgery...). A bit of bibliography led me to the
> recent Jean-Paul Fox's monograph, which did not let me with a useable way
> of simulating from this difficult posterior.
>
>
>> It's easy enough to add a prior. Andrew likes using
>> weakly informative priors in these situations,
>
>
> So do I, because 1) your priors being weak, you can succesfully argue that
> they do not carry much information into the posterior, and 2) being proper,
>
> you can (try to) use their resultant simulations in marginal likelihood
> estimation / Bayes factors / whatever model comparison method du jour you
> wish to use.
>
>
>> > imposing more informative distributions.
>
>> Gelman and Hill didn't propose a single solution. They discuss
>> several methods (setting one value based on prior knowledge,
>
> That's for another problem, i. e. to resolve the sign indetermination of
> the general model. Which does not arise in my case, since I can use subject
>
> matter knowledge to postulate that all my items are correlated the same
> way. Hence the lognormal distribution of the gammas.
>
>
>> of the prior parameters.
>
>> Let us know what happens. I don't have that much time to play
>> around with IRT parameterizations.
>
> Neither do I, alas, but I just have to : "When a researche gotta, he's
> gonna...". I was hoping that Stan's superior sampler could help me with a
> better simulation without having to wait hours (literally...) for
> convergence. But the lack of sampling from discrete distributions (if only
> for the missing data) makes it (currently) unfit for my current purposes.
>
> I will definitely try to compare JAGS and Stan (and possibly R Metropolis
> adaptive samplers such as the one exemplified in the adaptMCMC package,
> which performs surprisingly well with very little tuning) for the
> hierarchical models in the complete-cases case. But not now...
>
> Again, thank you for your attention and your time.
>
> Emmanuel
> Charpentier
>
>
>
> =============================================================================
> Topic: optimization level in RStan & discrete parameters in Stan
> Url: http://groups.google.com/group/stan-users/t/3884a2491bed2721
> =============================================================================
>
> ---------- 1 of 5 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 09 10:56PM -0400
> Url: http://groups.google.com/group/stan-users/msg/b4e8120d25ca54f0
>
> On 9/9/12 10:41 AM, Emmanuel Charpentier wrote:
>> the optimization level
>>
>> <http://code.google.com/p/stan/wiki/RStanGettingStarted#Modify_C++_compiler_optimization_level>
>> a hack.
>
>> OK. I'll look into changing that cleanly (maybe on a per-sesion basis,
>> through R's options()).
>
> I don't think that's possible. If you do figure out a way
> to do it, could you let us know?
>
> As a side benefit, -O3 will make other R apps compiled
> with R CMD SHLIB faster :-)
>
> Here's a thread on stackoverflow on the topic:
>
> http://stackoverflow.com/questions/5789673/r-c-optimization-flag-when-using-the-inline-package
>
>> Aha. I overlooked that in the manual. So currently Stan can't impute
>> *discrete* missing data. Nor assess the
>> distributions of *discrete* parameters (so much for model choice
>> algorithms such as umbrella sampling or hierarchical
>> model choice (à la Carlin & Chib (1995) or otherwise)...).
>
> Right. Unless you can marginalize out the
> discrete parameter, Stan can't handle it.
>
> I'm about to put a compile-time check to flag attempted
> uses of discrete params in the parser, which I forgot to
> do the first time around.
>
>> Any idea of the ETA for a discrete sampler ?
>
> We'll probably write a simple Gibbs for finite
> and slice for infinite sampler. Probably some
> time in the next six or nine months.
>
> - Bob
>
>
> ---------- 2 of 5 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 09 11:11PM -0400
> Url: http://groups.google.com/group/stan-users/msg/a934d4fa365f6c63
>
> n 9/9/12 5:26 AM, Emmanuel Charpentier wrote:
>> never initialized.
>
>> OK. Just a little problem : THERE IS NO theta VARIABLE IN THE MODELS. I
>> suspect that the compiler is goofing in building
>> it's error message.
>
> I see the confusion here.
>
> We should make that error message clearer.
> That's just a generic report from binomial_log, where
> the probability-of-success parameter is called theta.
> We should just get rid of the names altogether and refer
> to the names of the variable, like chance-of-success
> or scale or location.
>
> The problem for us is at the point that warning (and it's
> just a warning) gets reported, we've lost the link to the
> underlying variable in the Stan model. It would be possible
> to fix this by passing arguments representing the names of
> all the variables to the distributions or by trying to catch
> errors in the Stan code, but both are a lot of work and
> would make the code very ugly. Probably worth it in the
> longer term, as this is going to be very frustrating for
> users.
>
> - Bob
>
>
> ---------- 3 of 5 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 09 11:15PM -0400
> Url: http://groups.google.com/group/stan-users/msg/86a71d3d880456cc
>
> On 9/9/12 5:40 AM, Emmanuel Charpentier wrote:
>> initialization (to 0, rather thn sampling from the prior as in OpenBUGS)
>> seems to me more reasonable in many cases. I
>> love the stan() R function range of options : "0", "random" or
>> user-specified values. The best of ll worlds would be a
>> named list of parameters, and a default, which could be either "0" or
>> "random". But that might be treated at the R level.
>
> That's what we're aiming for. It'll take a bit of
> work and there are other fires to put out at the
> moment on basic error checking that we're going to
> address first.
>
> The biggest hurdle for the implementation is that the
> inits are on the unconstrained scale.
>
> I don't think JAGS always initializes to 0. It says
> in the manual that it takes reasonable values.
> People like to use very fat priors even when they
> have a good sense of what inits should be. Maybe
> Andrew has something to say about that -- he likes
> dispersed initializations.
>
> RStan and Stan use the same basic C++ code, but
> they have different command-line controls for marshaling
> the arguments. The C++ needs the data to be given
> to the model's class constructor and the constructor
> computes the transformed data, all of which then remains
> fixed throughout sampling.
>
> So far, we haven't had problems with Stan converging,
> even from pretty far-out inits. If they do get so
> bad that they result in NaN or infinite values, we
> resample.
>
> - Bob
>
>
> ---------- 4 of 5 ----------
> From: Andrew Gelman <gel...@stat.columbia.edu>
> Date: Sep 09 11:19PM -0400
> Url: http://groups.google.com/group/stan-users/msg/11ea72e6ff69d20
>
> On Sep 9, 2012, at 11:15 PM, Bob Carpenter wrote:
>
>> have a good sense of what inits should be. Maybe
>> Andrew has something to say about that -- he likes
>> dispersed initializations.
>
> If you're doing Gibbs, it's a good idea to avoid deterministic inits. For
> example, if you set theta_1,…,theta_J=0, then you might estimate tau=0, etc.
> This is probably less of an issue with Nuts or HMC. Still, I'd generally
> prefer some randomness in the inits (except when we want to do some pure
> comparison as in our speed comparisons with Jags).
>
>
> ---------- 5 of 5 ----------
> From: Emmanuel Charpentier <emanuel.c...@gmail.com>
> Date: Sep 10 01:18AM -0700
> Url: http://groups.google.com/group/stan-users/msg/56b72e4eab6acd8f
>
> Dear Bob,
>
> Le lundi 10 septembre 2012 05:11:30 UTC+2, Bob Carpenter a écrit :
>
>
>> longer term, as this is going to be very frustrating for
>> users.
>
>> Indeed...
>
> Re-reading my code (both R and Stan), I have been unable to identify the
> uninitialized parameters. I'm probably becoming myopic with this code, so
> I'll give it a rest, and try to debug it a few days later.
>
> Thank you again !
>
>
> Emmanuel Charpentier
>
>
>
> =============================================================================
> Topic: installing Stan
> Url: http://groups.google.com/group/stan-users/t/ca40e9dada65ed5e
> =============================================================================
>
> ---------- 1 of 4 ----------
> From: Shravan Vasishth <vasishth...@gmail.com>
> Date: Sep 09 11:15AM -0700
> Url: http://groups.google.com/group/stan-users/msg/85cfd4d792048baf
>
> *Sigh*. I managed to live my life out successfully without ever having to
> write (or even look inside, with some exceptions) a Makefile. Advice taken;
>
> I will get to know Makefiles, if only to avoid making a fool of myself in
> the future.
>
> On Sunday, September 9, 2012 12:07:58 PM UTC+2, Emmanuel Charpentier wrote:
>
>
> ---------- 2 of 4 ----------
> From: Daniel Lee <bea...@alum.mit.edu>
> Date: Sep 09 03:30PM -0400
> Url: http://groups.google.com/group/stan-users/msg/ba63641e7449c99b
>
> Shravan,
>
> Before you get too frustrated, all you have to do is to first change
> directory into the Stan directory, then type at the command line "make
> bin/libstan.a". The static library will be built for you at
> bin/libstan.a.
>
> There really isn't a need to get familiar with make, we (the Stan dev
> team) will try to encapsulate as much of the build pain through tools
> such as make.
>
>
>
> Daniel
>
>
>
>
> On Sun, Sep 9, 2012 at 2:15 PM, Shravan Vasishth
>
>
> ---------- 3 of 4 ----------
> From: Shravan Vasishth <vasishth...@gmail.com>
> Date: Sep 09 12:38PM -0700
> Url: http://groups.google.com/group/stan-users/msg/9f7d446c004f7797
>
> Hi Daniel,
>
> yes, I did follow the instructions without thinking about it the second
> time round! Everything works as expected.
>
> Thanks, and apologies for the pointless thread.
>
> On Sunday, September 9, 2012 9:30:30 PM UTC+2, Daniel Lee wrote:
>
>
> ---------- 4 of 4 ----------
> From: Bob Carpenter <ca...@alias-i.com>
> Date: Sep 09 10:47PM -0400
> Url: http://groups.google.com/group/stan-users/msg/51b95586f4deb219
>
> On 9/9/12 3:38 PM, Shravan Vasishth wrote:
>> Hi Daniel,
>
>> yes, I did follow the instructions without thinking about it the second
>> time round! Everything works as expected.
>
>> Thanks, and apologies for the pointless thread.
>
> Please don't feel it's your fault for not puzzling
> out how to run Stan from the command line. This is our fault
> for distributing software that doesn't install and
> run like a regular desktop app. (RStan does run
> like a regular R app.)
>
> We're trying to make it as easy as possible
> to use Stan, so this kind of question and feedback
> is really useful to us. We'll work on making the usage
> instructions clearer. Or ideally, figure out
> how to get it run relocatably (any help on
> the part of users who know how to do this better
> than us would also be greatly appreciated).
>
> In the meantime, please send us questions before
> you get too frustrated. We're more than happy to
> help.
>
> - Bob
>
>
>
>
>
--
Sent from my mobile device
------------------------------
Alexander Budzier
+1 (440) 253-9266
+44 (750) 210-7145
+49 (211) 136-4192
Bob Carpenter
Sep 10, 2012, 2:01:09 PM9/10/12
to stan-...@googlegroups.com
Please don't just respond to digests without editing
out irrelevant bits of the messages. I didn't even see
anything in this response, but if there was something,
could you repost it?
Thanks,
- Bob
On 9/10/12 1:46 PM, Alexander Budzier wrote:
> On 9/10/12, stan-...@googlegroups.com <stan-...@googlegroups.com> wrote:
>> =============================================================================
>> Today's Topic Summary
>> =============================================================================
>> ...
out irrelevant bits of the messages. I didn't even see
anything in this response, but if there was something,
could you repost it?
Thanks,
- Bob
On 9/10/12 1:46 PM, Alexander Budzier wrote:
> On 9/10/12, stan-...@googlegroups.com <stan-...@googlegroups.com> wrote:
>> =============================================================================
>> Today's Topic Summary
>> =============================================================================
Reply all
Reply to author
Forward
0 new messages