element-wise exponentiation of vectors (i.e. pow() for vectors

[glae...@gmail.com]

Dear all,

in my model I need to do an element-wise exponentiation of a row vector and that doesn't seem to be possible in Stan. Here's a snippet of my model code (full model is attached).

model {

  row_vector<lower=0,upper=1>[3] a[maxTrials+1]; // attention weights

  vector<lower=0,upper=1>[maxTrials] prob; // Prob(chosenAction)

  ...

  for (sub in 1:nSubjects) {

    ...

    for (t in 1:nTrials) {

      prob[t] <- dot_product(M[sub,t],a[t]^d[sub]) / ( sum(a[t]^d[sub]) + machine_precision() );

      ...

    }

  }

}

a[t] is a trial-specific row_vector that should be raised to the power of d[sub], which is a subject-specific model parameter before taking the dot_product with M[sub,t] (another vector, part of the data)

I keep getting this compiler error:

TRANSLATING MODEL 'rpdf' FROM Stan CODE TO C++ CODE NOW.

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 

  failed to parse Stan model 'rpdf' with error message:

EXPECTATION FAILURE - DIAGNOSTIC(S) FROM PARSER:

LOCATION OF PARSING ERROR (line = 120, position = 43):

PARSED:

model {

  /* Should these declarations be in the transformed_data block? */

  row_vector<lower=0,upper=1>[3] a[maxTrials+1]; // attention weights

  row_vector<lower=0,upper=1>[3] s[maxTrials]; // feedback signal

  vector<lower=0,upper=1>[maxTrials] prob; // Prob(chosenAction)

  r_mu ~ beta(A,B);

  r_prime_kappa ~ beta(A,B);

  p_mu ~ beta(A,B);

  p_prime_kappa ~ beta(A,B);

  d_prime_mu ~ beta(1,1);

  d_prime_kappa ~ beta(A,B);

  f_prime_mu ~ beta(A,B);

  f_prime_kappa ~ beta(A,B);

  for (sub in 1:nSubjects) {

    r[sub] ~ beta(r_a,r_b) T[L,U];

    p[sub] ~ beta(p_a,p_b) T[L,U];

    d_prime[sub] ~ beta(d_a,d_b) T[L,U];

    f_prime[sub] ~ beta(f_a,f_b) T[L,U];

    a[1] <- initA;    

    for (t in 1:nTrials[sub]) {

      prob[t] <- dot_product(M[sub,t],a[t]

EXPECTED: ")" BUT FOUND: 

^d[sub]) / ( sum(

I looked through the Stan User Manual (Section Built-in function), but could find the formal definition of the ^ operator. Shouldn't that be in there or did I miss it somewhere along the way?

I tried to change the line above to using the pow() function: 

prob[t] <- dot_product(M[sub,t],pow(a[t],d[sub])) / ( sum(pow(a[t],d[sub])) + machine_precision() );

but found out through a compiler error and the User Manual that pow() only take scalars (real) as inputs

Of course, one could break this exponentiation of a vector into it's elements by using a loop, but seem rather clumsy. Does anyone have an idea of how to run a power function on a vector?

Thanks a lot,

Jan

[Bob Carpenter]

You're right on both counts. There's no element-wise
exponentiation of row vectors. And we missed "operator^" in
the index. I'll add it to the next release of the manual:

https://github.com/stan-dev/stan/issues/786#issuecomment-51396580

You'll find the more informal doc in the v2.4 manual for ^ on p. 181
under "exponentiation" and in the table of operators. It doesn't
apply to vectors at all at the moment.

You can just assign to a local variable in a loop, though --- it's no slower in Stan,
just more painful to write. And in your case, you wouldn't want to
do

> prob[t] <- dot_product(M[sub,t],a[t]^d[sub]) / ( sum(a[t]^d[sub]) + machine_precision() );

anyway, because that'd require a[t]^d[sub] to be computed twice,
with twice as much work required for the gradients.

For efficiency, it helps to assign complicated subexpressions to
local variables and reuse them. In this case, the best way to write
this is as:

{
vector[K] ad;
real sum_ad;
sum_ad <- 0;
for (int k in 1:K) {
ad[k] <- a[t,k] ^ d[sub,k];
sum_ad_pp <- sum(ad) + machine_precision();
...
prob[t] <- dot_product(M[sub,t],ad[k]) / sum_ad_pp;
}

Even though it's ugly, there's no way around having to define ad in
a local variable like this.

- 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.
> <rpdf.stan>

[Jan Glaescher]

Dear Bob,

thanks again for your quick response and the suggestions for loop to solve the element-wise vector exponentiation.

Cheers,
Jan

Bob Carpenter

August 6, 2014 at 11:07 PM

glae...@gmail.com

August 6, 2014 at 10:34 PM

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

[Jan Glaescher]

Dear Bob,

I tested your suggestion, but in fact the exponentiation with the ^ operator in your loop didn't work. I had to use pow(a[sub,t],d[sub]) function to get past the compiler errors. Not sure, why that is because a has the type

row_vector<lower=0,upper=1>[3] a[maxTrials+1];

and therefore I thought that a[t,k] should be a scalar (real). d[sub] is an individual model parameters sampled from the overarching beta distribution. Actually, d_prime[sub] is sampled from the beta distribution in the model block at the beginning of the subjects loop and d[sub] is then rescaled to a range of 0-5. The latter is done in a loop in the transformed parameter block (see attached model code). Is that OK or could this be the problem?

Nevertheless, with the changes above the model compiled fine (with a few warnings) and started the sampling then terminated with the error message:
SAMPLING FOR MODEL 'rpdf' NOW (CHAIN 1).
Initialization 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 am not sure what Stan is trying to tell me. I am creating chain-specific initial values for each parameter in R by drawing a random number from same uniform beta distribution (truncated to 0.001 and 0.999). Following your previous advice I only set the lower and upper values for the beta distribution to 0 and 1, (but the errors also occur when I truncate to 0.001 and 0.999). Finally, I am not sure, how to approach the reparameterization.

Complete error log and model code are attached.

Thanks a lot for your time and effort.
Jan

Bob Carpenter

August 6, 2014 at 11:07 PM

You're right on both counts.  There's no element-wise
exponentiation of row vectors.  And we missed "operator^" in
the index.  I'll add it to the next release of the manual:

  https://github.com/stan-dev/stan/issues/786#issuecomment-51396580

glae...@gmail.com

August 6, 2014 at 10:34 PM

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

[Bob Carpenter]

You should be able to use (x ^ y) wherever you can use pow(x,y).

If you find a case where that doesn't work, could you send us the

non-compiling model?  It'd help us make sure there's not a bug somewhere

in our parser code.

Failed initialization is due to the fact that the initial values

you are providing have 0 probability (-infinity on the log scale).

The only thing I can suggest is making sure everything's well defined

at the point it's used.  You can use print-based debugging to

see what's going wrong.  If you print lp__ (ignore the warnings)

after each statement, you should be able to see where things go wrong.

- Bob

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+un...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

<error.log><rpdf.stan>

[Jan Glaescher]

Dear Bob,

attached is the Stan model which causes the compile error at the point I am using ^ I am also attaching the latest error.log corresponding to this. Do you need the data or the R script to verify? It works with rpdf.stan model (also attached), in which I just change a[t,k] ^ d[sub] to pow(a[t,k],d[sub]) and similar lines below that. FYI, I am using rstan on R 3.1.1 on a Macbook Pro.

I managed to get past the initialization failure error. You were right that some of my initial values were not well specified. However, now I am getting the following run-time error:

SAMPLING FOR MODEL 'rpdf' NOW (CHAIN 1).

Error : Rejecting user-specified initialization because of vanishing density.

error occurred during calling the sampler; sampling not done

Could you tell me what this means?

It would be tremendously helpful for the beginners in Stan, if the User Manual contained a section with the most common errors, their explanation and suggestions of how to address them. I know that is maybe asking too much, because you developers have already done so much, but maybe this could be a community effort (I would certainly be willing to contribute.)

Thanks again.
Jan

Bob Carpenter

August 8, 2014 at 1:42 AM

Jan Glaescher

August 8, 2014 at 1:06 AM

Dear all,

in my model I need to do an element-wise exponentiation of a row vector and that doesn't seem to be possible in Stan. Here's a snippet of my model code (full model is attached).

model {
  row_vector<lower=0,upper=1>[3] a[maxTrials+1]; // attention weights
  vector<lower=0,upper=1>[maxTrials] prob; // Prob(chosenAction)
  ...
  for (sub in 1:nSubjects) {
    ...
    for (t in 1:nTrials) {
      prob[t] <- dot_product(M[sub,t],a[t]^d[sub]) / ( sum(a[t]^d[sub]) + machine_precision() );
      ...
    }
  }
}

a[t] is a trial-specific row_vector that should be raised to the power of d[sub], which is a subject-specific model parameter before taking the dot_product with M[sub,t] (another vector, part of the data)

I keep getting this compiler error:

TRANSLATING MODEL 'rpdf' FROM Stan CODE TO C++ CODE NOW.
Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'rpdf' with error message:
EXPECTATION FAILURE - DIAGNOSTIC(S) FROM PARSER:



LOCATION OF PARSING ERROR (line = 120, position = 43):

PARSED:

model {

  /* Should these declarations be in the transformed_data block? */
  row_vector<lower=0,upper=1>[3] a[maxTrials+1]; // attention weights
  row_vector<lower=0,upper=1>[3] s[maxTrials]; // feedback signal
  vector<lower=0,upper=1>[maxTrials] prob; // Prob(chosenAction)

  r_mu ~ beta(A,B);
  r_prime_kappa ~ beta(A,B);

  p_mu ~ beta(A,B);
  p_prime_kappa ~ beta(A,B);

  d_prime_mu ~ beta(1,1);
  d_prime_kappa ~ beta(A,B);

  f_prime_mu ~ beta(A,B);
  f_prime_kappa ~ beta(A,B);


  for (sub in 1:nSubjects) {

    r[sub] ~ beta(r_a,r_b) T[L,U];
    p[sub] ~ beta(p_a,p_b) T[L,U];
    d_prime[sub] ~ beta(d_a,d_b) T[L,U];
    f_prime[sub] ~ beta(f_a,f_b) T[L,U];

    a[1] <- initA;    

    for (t in 1:nTrials[sub]) {
      prob[t] <- dot_product(M[sub,t],a[t]

EXPECTED: ")" BUT FOUND: 

^d[sub]) / ( sum(

I looked through the Stan User Manual (Section Built-in function), but could find the formal definition of the ^ operator. Shouldn't that be in there or did I miss it somewhere along the way?

I tried to change the line above to using the pow() function: 

prob[t] <- dot_product(M[sub,t],pow(a[t],d[sub])) / ( sum(pow(a[t],d[sub])) + machine_precision() );

but found out through a compiler error and the User Manual that pow() only take scalars (real) as inputs

Of course, one could break this exponentiation of a vector into it's elements by using a loop, but seem rather clumsy. Does anyone have an idea of how to run a power function on a vector?

Thanks a lot,
Jan

-- 
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.
<rpdf.stan>

[Christian Roy]

Hi Jan,

By experience the error message "Rejecting user-specified initialization because of vanishing density" is due to bad initialization values. I triggered this error message in the past by creating starting values that created -inf values on the log scale or on the logit scale. 

Your model seems to use the beta distribution. I never used the beta distribution with Stan but based on my experience with JAGS and WinBUGS I would guess that your initialization values a probably too close to zero for some parameters. You may want to try starting values that are away from zero to see if you can get the model over the bump. 

HTH

Chris

----------------------------

[Jan Glaescher]

Hi Chris,

thanks for the suggestion. I change my R code and sample the initial values from a beta(30,30) distribution now. Almost all values are between 0.4-0.6, none is below 0.2 or above 0.8. I ran it 3 times, but I am still getting the same error message on every run.

Is there a particular reason, why you don't use beta distributions with Stan? Are they somehow implemented differently than in BUGS/JAGS? Most of my parameters have a natural boundaries of 0 and 1 and I like the flexibility of the beta distribution in that it can accommodate very skewed data. Or would you suggest to use a different distribution?

Thanks a lot,
Jan

Christian Roy

August 8, 2014 at 5:29 PM

Hi Jan,

By experience the error message "Rejecting user-specified initialization because of vanishing density" is due to bad initialization values. I triggered this error message in the past by creating starting values that created -inf values on the log scale or on the logit scale. 

Your model seems to use the beta distribution. I never used the beta distribution with Stan but based on my experience with JAGS and WinBUGS I would guess that your initialization values a probably too close to zero for some parameters. You may want to try starting values that are away from zero to see if you can get the model over the bump. 

HTH

Chris

----------------------------

  

--
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/Q7zADp0Ho9E/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Jan Glaescher

August 8, 2014 at 4:52 PM

[Bob Carpenter]

Thanks.

That's definitely a parsing bug.  The exponentiation operator currently

has the wrong precedence --- it should bind less tightly than indexing. 

I created a high-priority issue to fix it with a simplified model:

https://github.com/stan-dev/stan/issues/834

For now, you'll have to use pow() or wrap the arguments to the ^ operator

in parens; e.g., (a[1])^(b[1]).

- Bob

<rpdf.stan>

-- 
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+un...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

<rpdf_with_exp.stan><error.log>

[Bob Carpenter]

First, I will create an appendix for the manual listing the error

messages and trying to explain them better.  The problem is that they

can often come from multiple sources.  The "vanishing density" and

"informational message" messages are the ones users are most often

confused about.  We'd ideally like to make the messages themselves better

rather than requiring a trip back to the manual, but that's harder than

it sounds --- we've never had any luck explaining the informational message

and users continue to construe it as a bug rather than a warning, no matter

how we word it.  We're taking another pass in the next release.

We do have beta distributions and Stan and would use them wherever

appropriate.  The reason you don't see us coding models the same way

as in JAGS is due to a number of reasons.  First, there's no benefit to conjugate 

priors in Stan models, so we don't use a lot of inverse gamma or Wisharts,

as there are better approaches.  Second, for uniform distributions on bounded

variables, we don't need a prior --- it's already uniform.  So you don't

see beta(1,1) used anywhere as it's redundant.  Third, we often work with

GLMs rather than straight up binomial parameters, and those don't involve

beta densities.  But by all means use them if they make sense for your model.

We also don't use a lot of Dirichlet priors, but that's because they're very

weak in not specifying any covariance structure.

Have you tried my suggestion of trying to diagnose where the log prob

goes to -infinity by using print("lp__=",lp__) statements?  

Your code's too long for me to try to debug by inspection, but you should

also make sure to define every variable before you use it --- the order

of statements matters in Stan.  And make sure every argument is legal to

every function called.

I'd also be careful about your T[L,U].  That's not necessary in Stan

if L and U are constants and the variable is already constrained to lie

between L and U.  You have a problem in that you define

real<lower=0,upper=1> foo;

and then pair that with

 foo ~ beta(alpha,beta) T[0.001,0.999];

This is a very bad thing to do in Stan because the constraint on

the variable foo ([0,1]) does not match the support of the

distribution specified, which is truncated to [0.001,0.999].

- Bob

<rpdf.stan>

-- 
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+un...@googlegroups.com.

[Christian Roy]

Hi Jan,

The only reason I  have never used the beta distribution in Stan is because I used the beta distribution to estimate detection in occupancy modelling or abundance modelling and it's impossible to do those kind of model in Stan (for the moment...maybe forever....but I keep dreaming that it will be possible one day).

I've never worked with the type of model you are using right now but I think the problem is that your starting values are actually defining the parameters of some of your beta distribution. 

For example "r_a" is a parameter for a beta distribution in the model. 

yet "r_a" is the product itself of "r_mu" and "r_kappa"  and r_kappa" itself is defined from "r_prime_kappa".

So even if you try to keep your values away from the extreme you could end up with some kind of highly constrained beta distribution (multiplying a lot of values between 0 and 1 is bound to create some problem at some point). 

Basically I feel like the model will be really tricky to initialize no matter what software you use. 

I have some really painful memories of trying to initialize models like that in JAGS. In those situation I usually started backward. I created a function that would define the most troublesome values and derive the hyper priors from those values.  

 Chris

-------------------

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Jan Glaescher]

Hi Chris,

thanks a lot for sharing your experience and the suggestions. I didn't have a lot of trouble initializing this model in JAGS, but I wonder if the initialization would benefit, if I'd sample the parameter of the group distribution in the model block not from an uninformative beta(1,1), but rather a more symmetrical beta(10,10, or beta(20,20) even though this might entail running the sample for a bit longer until convergence is reached. Any thoughts?

Thanks a lot,
Jan

Christian Roy

August 8, 2014 at 7:23 PM

Hi Jan,

The only reason I  have never used the beta distribution in Stan is because I used the beta distribution to estimate detection in occupancy modelling or abundance modelling and it's impossible to do those kind of model in Stan (for the moment...maybe forever....but I keep dreaming that it will be possible one day).

I've never worked with the type of model you are using right now but I think the problem is that your starting values are actually defining the parameters of some of your beta distribution. 

For example "r_a" is a parameter for a beta distribution in the model. 

yet "r_a" is the product itself of "r_mu" and "r_kappa"  and r_kappa" itself is defined from "r_prime_kappa".

So even if you try to keep your values away from the extreme you could end up with some kind of highly constrained beta distribution (multiplying a lot of values between 0 and 1 is bound to create some problem at some point). 

Basically I feel like the model will be really tricky to initialize no matter what software you use. 

I have some really painful memories of trying to initialize models like that in JAGS. In those situation I usually started backward. I created a function that would define the most troublesome values and derive the hyper priors from those values.  

 Chris

-------------------

--
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/Q7zADp0Ho9E/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Bob Carpenter

August 8, 2014 at 7:09 PM

Jan Glaescher

August 8, 2014 at 5:53 PM

[Bob Carpenter]

On Aug 8, 2014, at 1:23 PM, Christian Roy <christia...@gmail.com> wrote:

Hi Jan,

The only reason I  have never used the beta distribution in Stan is because I used the beta distribution to estimate detection in occupancy modelling or abundance modelling and it's impossible to do those kind of model in Stan (for the moment...maybe forever....but I keep dreaming that it will be possible one day).

Hang on!  It's not impossible to do these kinds of models in Stan.

We've fit all kinds of discrete population models.  For example,

there are some defined in this branch:

feature/issue-451-mark-recapture-egs

You can view it using GitHub's online viewer:

https://github.com/stan-dev/stan/tree/feature/issue-451-mark-recapture-egs/src/models/misc

The examples are in src/models/misc under mark-recapture and occupancy directories.

We also have models that Fraenzi Korner-Nievergelt generalized for us as

part of her book on population models.

What is impossible is doing it with discrete parameters.  But

if you can marginalize out the discrete parameters, it's much more

efficient to sample these models in Stan.  

I need to add some more examples to the manual!

I've never worked with the type of model you are using right now but I think the problem is that your starting values are actually defining the parameters of some of your beta distribution. 

For example "r_a" is a parameter for a beta distribution in the model. 

yet "r_a" is the product itself of "r_mu" and "r_kappa"  and r_kappa" itself is defined from "r_prime_kappa".

This shouldn't be a problem.  The only time you're going to have

a problem is if the variables are out of support or are so extreme

you get underflow or overflow.  That'll only happen with the beta

as the parameters go to 0 or infinity.

So even if you try to keep your values away from the extreme you could end up with some kind of highly constrained beta distribution (multiplying a lot of values between 0 and 1 is bound to create some problem at some point). 

Also, have you tried just letting Stan initialize?  Stan usually

has much better convergence properties than JAGS or BUGS from

extreme initial values. 

- Bob

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.

[Bob Carpenter]

If the initialization works in JAGS, it should work in Stan.

My guess is that there's a bug in your Stan code for the model

that's causing the log prob to be -infinity.  A likely cause is

not having all variables defined.  

I would very very strongly suggest following my earlier advice

and trying to track down where the model goes bad using print()

statements.

- 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+un...@googlegroups.com.

[Christian Roy]

Well...that's very interesting Bob. I'll definitely have to try these models.

And yes usually Stan is quite good at picking it's own starting values. So if you let Stan initialize and you still get an error message there's probably something wrong with the model. 

Chris

----------------------

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Jan Glaescher]

Thanks you both. I do get the "Initialization Failure" error if I let Stan do the inilization. So there must be something wrong with the model.

Bob, should I print lp__ after each line of code or just after each sampling statement? Sorry, I realize these are all beginner's questions ...

- Jan

Bob Carpenter

August 8, 2014 at 8:41 PM

[Bob Carpenter]

On Aug 8, 2014, at 3:02 PM, Jan Glaescher <glae...@gmail.com> wrote:

> Thanks you both. I do get the "Initialization Failure" error if I let Stan do the inilization. So there must be something wrong with the model.
>
> Bob, should I print lp__ after each line of code or just after each sampling statement? Sorry, I realize these are all beginner's questions ...
>

You can try divide and conquer. Print it once in the middle, and if
it's good, look at the second half and if it's bad look at the first half.

You may also wind up wanting to print out some of the other values
to make sure you're passing legal values into all the functions.

The other problem is that we need to enhance the error messages. There's
also a problem on our side in that the initialization just swallows any
informative message and just gives you the inscrutable "vanishing density"
message at the end.

- Bob

[Jan Glaescher]

Dear Bob,

I added a print statement at the end of the model block and re-compiled the model: the log probability jumps in the range of about -2271 to about -2305. None of these values look bad to me (i.e. not NaN or -inf), but the sampling still stop with the failed initialization message after 100 attempts (see attached file lp.txt). I also tried this in the middle of the model block with the same results. I also had a more comprehensive set of print statements that printed current parameter values, but those also seem to be legal.

In the attached file you'll also see some compiler warnings, which sounded harmless to me (but I don't know much about C++ so I maybe wrong).

Is there anything else that I can try in order to get Stan to sample?

Thanks a lot,
Jan

Bob Carpenter

August 8, 2014 at 9:46 PM

You can try divide and conquer. Print it once in the middle, and if
it's good, look at the second half and if it's bad look at the first half.

You may also wind up wanting to print out some of the other values
to make sure you're passing legal values into all the functions.

The other problem is that we need to enhance the error messages. There's
also a problem on our side in that the initialization just swallows any
informative message and just gives you the inscrutable "vanishing density"
message at the end.

- Bob

Jan Glaescher

August 8, 2014 at 9:02 PM

Thanks you both. I do get the "Initialization Failure" error if I let Stan do the inilization. So there must be something wrong with the model.

Bob, should I print lp__ after each line of code or just after each sampling statement? Sorry, I realize these are all beginner's questions ...

- Jan

Bob Carpenter

August 8, 2014 at 8:41 PM

If the initialization works in JAGS, it should work in Stan.

My guess is that there's a bug in your Stan code for the model

that's causing the log prob to be -infinity.  A likely cause is

not having all variables defined.  

I would very very strongly suggest following my earlier advice

and trying to track down where the model goes bad using print()

statements.

- Bob

On Aug 8, 2014, at 2:37 PM, Jan Glaescher <glae...@gmail.com> wrote:

--
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/Q7zADp0Ho9E/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Jan Glaescher

August 8, 2014 at 8:37 PM

Hi Chris,

thanks a lot for sharing your experience and the suggestions. I didn't have a lot of trouble initializing this model in JAGS, but I wonder if the initialization would benefit, if I'd sample the parameter of the group distribution in the model block not from an uninformative beta(1,1), but rather a more symmetrical beta(10,10, or beta(20,20) even though this might entail running the sample for a bit longer until convergence is reached. Any thoughts?

Thanks a lot,
Jan

Christian Roy

August 8, 2014 at 7:23 PM

[Bob Carpenter]

[changed subject]

I thought the vanishing density was when the densities

overflowed or became NaN.  I changed the subject to get

some others' attention.

You should also try to print out the gradients at your initial

point and make sure they're all finite.  I should have suggested

that first.  It's possible that the log probability is OK but the 

gradients are bad. 

If those look OK, then someone else is going to have jump in

or I'll have to dig deeper from your data and model.

- 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+un...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

<lp.txt>

[Jan Gläscher]

Dear Bob, 

thanks for your comments. Can you quickly tells me what the relevant variables are that I should print when trying to assess the gradients?

Thanks,

Jan

On Aug 11, 2014, at 19:34, Bob Carpenter <ca...@alias-i.com> wrote:

[changed subject]

I thought the vanishing density was when the densities

overflowed or became NaN.  I changed the subject to get

some others' attention.

You should also try to print out the gradients at your initial

point and make sure they're all finite.  I should have suggested

that first.  It's possible that the log probability is OK but the 

gradients are bad. 

If those look OK, then someone else is going to have jump in

or I'll have to dig deeper from your data and model.

- Bob

On Aug 11, 2014, at 12:19 PM, Jan Glaescher <glae...@gmail.com> wrote:

Dear Bob,

I added a print statement at the end of the model block and re-compiled the model: the log probability jumps in the range of about -2271 to about -2305. None of these values look bad to me (i.e. not NaN or -inf), but the sampling still stop with the failed initialization message after 100 attempts (see attached file lp.txt). I also tried this in the middle of the model block with the same results. I also had a more comprehensive set of print statements that printed current parameter values, but those also seem to be legal.

In the attached file you'll also see some compiler warnings, which sounded harmless to me (but I don't know much about C++ so I maybe wrong).

Is there anything else that I can try in order to get Stan to sample?

Thanks a lot,
Jan

<postbox-contact.jpg>

Bob Carpenter
August 8, 2014 at 9:46 PM

You can try divide and conquer. Print it once in the middle, and if
it's good, look at the second half and if it's bad look at the first half.

You may also wind up wanting to print out some of the other values
to make sure you're passing legal values into all the functions.

The other problem is that we need to enhance the error messages. There's
also a problem on our side in that the initialization just swallows any
informative message and just gives you the inscrutable "vanishing density"
message at the end.

- Bob

<postbox-contact.jpg>

Jan Glaescher
August 8, 2014 at 9:02 PM
Thanks you both. I do get the "Initialization Failure" error if I let Stan do the inilization. So there must be something wrong with the model.

Bob, should I print lp__ after each line of code or just after each sampling statement? Sorry, I realize these are all beginner's questions ...

- Jan

<postbox-contact.jpg>

Bob Carpenter
August 8, 2014 at 8:41 PM
If the initialization works in JAGS, it should work in Stan.
My guess is that there's a bug in your Stan code for the model
that's causing the log prob to be -infinity.  A likely cause is
not having all variables defined.  

I would very very strongly suggest following my earlier advice
and trying to track down where the model goes bad using print()
statements.

- Bob


On Aug 8, 2014, at 2:37 PM, Jan Glaescher <glae...@gmail.com> wrote:


-- 
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/Q7zADp0Ho9E/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

<postbox-contact.jpg>

Jan Glaescher
August 8, 2014 at 8:37 PM
Hi Chris,

thanks a lot for sharing your experience and the suggestions. I didn't have a lot of trouble initializing this model in JAGS, but I wonder if the initialization would benefit, if I'd sample the parameter of the group distribution in the model block not from an uninformative beta(1,1), but rather a more symmetrical beta(10,10, or beta(20,20) even though this might entail running the sample for a bit longer until convergence is reached. Any thoughts?

Thanks a lot,
Jan 

<compose-unknown-contact.jpg>

Christian Roy
August 8, 2014 at 7:23 PM
Hi Jan,

The only reason I  have never used the beta distribution in Stan is because I used the beta distribution to estimate detection in occupancy modelling or abundance modelling and it's impossible to do those kind of model in Stan (for the moment...maybe forever....but I keep dreaming that it will be possible one day).

I've never worked with the type of model you are using right now but I think the problem is that your starting values are actually defining the parameters of some of your beta distribution. 

For example "r_a" is a parameter for a beta distribution in the model. 

yet "r_a" is the product itself of "r_mu" and "r_kappa"  and r_kappa" itself is defined from "r_prime_kappa".

So even if you try to keep your values away from the extreme you could end up with some kind of highly constrained beta distribution (multiplying a lot of values between 0 and 1 is bound to create some problem at some point). 

Basically I feel like the model will be really tricky to initialize no matter what software you use. 

I have some really painful memories of trying to initialize models like that in JAGS. In those situation I usually started backward. I created a function that would define the most troublesome values and derive the hyper priors from those values.  

 Chris
-------------------




-- 
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/Q7zADp0Ho9E/unsubscribe.
To unsubscribe from this group and all its topics, 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 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+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<lp.txt>

[Bob Carpenter]

I couldn't tell if you're using CmdStan or RStan.
Both have different calls.

If it's CmdStan, just use

% ./model-name diagnose

If it's RStan, add

test_grad=TRUE

option to the stan() call and it'll only print the gradients out
and not try to sample.

- Bob

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

>>> <lp.txt>
>>
>>
>> --
>> 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/Q7zADp0Ho9E/unsubscribe.
>> To unsubscribe from this group and all its topics, 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 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.

[Jan Glaescher]

After compilation it says:

TESTING GRADIENT FOR MODEL 'rpdf' NOW (CHAIN 1).

But it never prints out any information about the gradients. It goes right into initialization (and hence prints out the log probability that I put in the print statement (see attachment).

I did the standard rstan setup (R 3.1.1 on Macbook Pro). Do I need to set an additional debug flag or something else during installation to get the gradient information?

Thanks,
Jan

Bob Carpenter

August 11, 2014 at 8:56 PM

I couldn't tell if you're using CmdStan or RStan.
Both have different calls.

If it's CmdStan, just use

% ./model-name diagnose

If it's RStan, add

test_grad=TRUE

option to the stan() call and it'll only print the gradients out
and not try to sample.

- Bob

Bob Carpenter

August 11, 2014 at 7:34 PM

[changed subject]

I thought the vanishing density was when the densities

overflowed or became NaN.  I changed the subject to get

some others' attention.

You should also try to print out the gradients at your initial

point and make sure they're all finite.  I should have suggested

that first.  It's possible that the log probability is OK but the 

gradients are bad. 

If those look OK, then someone else is going to have jump in

or I'll have to dig deeper from your data and model.

- Bob

On Aug 11, 2014, at 12:19 PM, Jan Glaescher <glae...@gmail.com> wrote:

--
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/Q7zADp0Ho9E/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Bob Carpenter

August 8, 2014 at 9:46 PM

You can try divide and conquer. Print it once in the middle, and if
it's good, look at the second half and if it's bad look at the first half.

You may also wind up wanting to print out some of the other values
to make sure you're passing legal values into all the functions.

The other problem is that we need to enhance the error messages. There's
also a problem on our side in that the initialization just swallows any
informative message and just gives you the inscrutable "vanishing density"
message at the end.

- Bob

Jan Glaescher

August 8, 2014 at 9:02 PM

Thanks you both. I do get the "Initialization Failure" error if I let Stan do the inilization. So there must be something wrong with the model.

Bob, should I print lp__ after each line of code or just after each sampling statement? Sorry, I realize these are all beginner's questions ...

- Jan

Bob Carpenter

August 8, 2014 at 8:41 PM

[Bob Carpenter]

On Aug 11, 2014, at 3:16 PM, Jan Glaescher <glae...@gmail.com> wrote:

> After compilation it says:
>
> TESTING GRADIENT FOR MODEL 'rpdf' NOW (CHAIN 1).
>
> But it never prints out any information about the gradients. It goes right into initialization (and hence prints out the log probability that I put in the print statement (see attachment).
>
> I did the standard rstan setup (R 3.1.1 on Macbook Pro). Do I need to set an additional debug flag or something else during installation to get the gradient information?

No. I think this may just be an undesirable side-effect of
using multiple initializations. Can you try init=0 and
see what that does?

- Bob

[Jan Glaescher]

Still no gradient information, but a new error message:

TESTING GRADIENT FOR MODEL 'rpdf' NOW (CHAIN 1).

log probability: -1879.82
Error : Error during initialization with 0: vanishing density.

error occurred during calling the sampler; sampling not done

Stan model 'rpdf' does not contain samples.

The log probability line is coming from my print statement

Jan

Bob Carpenter

August 11, 2014 at 9:41 PM

No. I think this may just be an undesirable side-effect of
using multiple initializations. Can you try init=0 and
see what that does?

- Bob

Jan Glaescher

August 11, 2014 at 9:16 PM

After compilation it says:

TESTING GRADIENT FOR MODEL 'rpdf' NOW (CHAIN 1).

But it never prints out any information about the gradients. It goes right into initialization (and hence prints out the log probability that I put in the print statement (see attachment).

I did the standard rstan setup (R 3.1.1 on Macbook Pro). Do I need to set an additional debug flag or something else during installation to get the gradient information?

Thanks,
Jan

Bob Carpenter

August 11, 2014 at 8:56 PM

I couldn't tell if you're using CmdStan or RStan.
Both have different calls.

If it's CmdStan, just use

% ./model-name diagnose

If it's RStan, add

test_grad=TRUE

option to the stan() call and it'll only print the gradients out
and not try to sample.

- Bob

Jan Gläscher

August 11, 2014 at 8:47 PM

Dear Bob, 

thanks for your comments. Can you quickly tells me what the relevant variables are that I should print when trying to assess the gradients?

Thanks,

Jan

On Aug 11, 2014, at 19:34, Bob Carpenter <ca...@alias-i.com> wrote:

[glae...@gmail.com]

Dear Bob and others,

as discussed previously, I am posting my model file, the rstan code, a small piece of my data, and an error message.

init=0 alos give me the vanishing density as I posted before. Letting Stan choose it's own initial values also results in a vanishing density error during intialization.

Ias this somehow the same issue that Tamas Papp posted today? Should I provide more print commands of lp__ ?

Thanks a lot,
Jan

--

Jan Glaescher
glae...@gmail.com

[Bob Carpenter]

Yes, the lack of ability to debug by printing lp__ is
related. What Tamas did was print the intermediate values
to see where things went wrong.

I'm working right now on a fix for printing lp__.

If there was supposed to be an attachment, I didn't get it.

- Bob

[glae...@gmail.com]

Ah yes, sorry. Pressed send too early again.

Jan

[glae...@gmail.com]

Dear Bob et al.,

I am reopening this thread from about 8 months ago. Back then I received a “vanishing density” error message during initialization that I wasn’t ab to resolve. In the meantime I continued working with JAGS, where my model compiles and samples without a problem.

Now with Stan 2.6 released I thought I’d give this a try again. The model still won’t initialize, but now I get a new error message (full log is attached):

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

(see attached error_stan26.log)

Afaik, I did the first two recommendations: I initialized all variables in the parameters block (see run_model.r, attached) and I also constrained the underlying values. I also changed the parameterization of variance of the group distribution to a half-cauchy (rpdf_cauchy.stan). All of this had no effect. I still get the above mentioned error message.

I have two questions:
1. Do you or anyone else have another idea that I can try to get Stan to estimate this model.

2. Following, Ahn et al., J Neurosci, Psycho Econ, 2011, I parameterized the group distribution with a beta distribution that was restricted to a uni-modal distribution by limiting the variance parameter. The reason for this choice was that two parameters have natural boundaries between 0 and 1. The two parameters (a,b) of this beta were itself parameterized in terms of mean and variance, which were also drawn from a beta. (see attached code rpdf.stan for the exact implementation).

Do you think that this is a good way to parameterize the model or should I use a different group distribution? Is is possible that the choice of distributions is causing this initialization failure in Stan?

Thanks a lot,
Jan

[Andrew Gelman]

Hi, I’m just coming in from the outside here and have not looked at the model in detail, but this jumped out at me:

real<lower=0.001,upper=0.999>

In general we frown on this sort of hard constraint. Generally we prefer real<lower=0,upper=1> and if you want to keep away from the boundary, just throw in a beta(2,2) prior or something like that.

A

> --

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

> <error_stan26.log><data_small.dump><rpdf_cauchy.stan><rpdf.stan><run_model.r>

[Bob Carpenter]

Thanks for sticking with this and following up. We'll
get this issue sorted out.

I don't understand how you can see the output of
your print statement with the model failing. I'm going
to ping the stan-dev list with a query. I'd thought we
fixed the issue where we were swallowing errors from init
without printing them, too, but apparently not. I'm going
to follow up on the dev list and if that doesn't work, trace
through the CmdStan code.

I don't have time to debug the program for you this weekend, but
should be able to get to it by next weekend (I'm on the road and
need to prepare a class). I'm first going to follow up on the dev
list to see what's up with the interface or look into the code if that
doesn't work.

If you want to try to debug yourself, the thing to do is add more
print statements for intermediate quantities to see what's
going wrong. Or even better, start with a simple model that
works then build up to something more complicated a step
at a time.

To answer your numbered questions:

1. Not yet.

2. No, the choice of priors shouldn't affect initialization unless something
is out of range. The only thing that should prevent initialization is a non-finite
log density or gradient somewhere. I won't have an opinion on priors until I've
understood the model.

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

[glae...@gmail.com]

Dear Andrew and Bob,

thanks for taking up this issue again. I really appreciate the prompt and immediate support thatt you provide on this list.

@Andrew
Yes, this is a left-over from running the model in JAGS. I also wanted to be really sure that I don’t run into extreme values in the beta, when trying to figure our this problem. But I will certainly use natural bounds and vague, but not flat priors when I really get this model to run in Stan.

@Bob
Thanks a lot for getting the devs involved. I am also gone until the middle of next week, so no need to hurry this. In the mean time, I’ll will make the model simpler by eliminating some parameters and possibly also the choice function in the model, which now selected similarly to a softmax.

Are there any other internal variable (besides lp__) that you need for debugging?

Btw, the model is the one from Bishara et al., J Math Psychol, 2010 for the Wisconsin Card Sorting Test. The description of the model is terrible though (within the text), but it works.
doi: 10.1016/j.jmp.2008.10.002
free @ http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2872109/

Thanks a lot,
Jan

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

[glae...@gmail.com]

Dear Bob,

have you anything from the Stan Devs about this issue ("Initialization between (-2, 2) failed after 100 attempts." yet?

In the mean time I can report that reducing the model's complexity by replacing one free parameter after another with a constant until the model does not have any free parameters is still resulting in the same error that I reported before.

Furthermore, to be really sure, I included lots of print statements in the Stan file for debugging purposes to check that every variable that is calculated by Stan is within a valid range(e.g. that probabilities are not accidentally below 0 or above 1. This is indeed the case, so basically I think that the Stan code that I attached earlier is OK and produces valid numbers.

If you have any further idea that I should try to sort this out, please let me know.

Thanks a lot,
Jan


> On Feb 14, 2015, at 3:09 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>

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

[Daniel Lee]

Hi Jan,

It's definitely a bug and I'm trying to track it down. It's pretty deep in the guts of Stan, so it's taking longer to find than I expected. Nothing to do right now. It's happening in between iterations, so you won't be able to track it with print statements, unfortunately.

Daniel

[Bob Carpenter]

Daniel --- is there an issue and/or branch for this yet?
I don't have anything urgent to work on at the moment, so
I could help out.

In the meantime, I'll update the manual for the 2.6.2 release.

- Bob

[glae...@gmail.com]

Dear Daniel and Bob,
thanks for continuing to looking into this. I appreciate it very much. Jan