Initialization issues when modeling Binary decision with lapse rate

[Shaorong Yan]

Hi all, 

     Sorry to bother...  

     I recently tried to use rstan to model a binary decision process with the basic structure as follows: 

     

     ...

     prob = (1 - lapse) * 1 / (1+ exp( (a-x) / b ) ) + lapse * 0.5

     

     y ~ binomial( prob )

     ....

 

     I have constrained the lapse rate to be within 0 and .1, and also the probability p. 

     Further complication is that I am also modeling how lapse rate change over trials, i.e. 

 

     ...

     lapse_sum = lapse + k_l * f(trial) 

     prob = (1 - lapse_sum) * 1 / (1+ exp( (a -x) / b ) ) + lapse_sum * 0.5

     ...

     In the model, how I am currently specifying the model is 

     ...

     parameter{

     real  a  ; 

     real  b  ;

     real<lower = 0, upper = 0.1>  lapse ;  

     real k_l ;

     }

      

     transformed parameters{

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

     real>lower = 0, upper = .1> lapse_sum ;  

     prob = (1 - lapse_sum) * 1 / (1+ exp( (a -x) / b ) ) + lapse_sum * 0.5

     }

     model{

     a ~  normal(0,  100)  ;

     b ~  normal(0,  100)   ;

     lapse ~  uniform(0, .1)   ;

     k_l ~ normal(0, .1)

     

     y ~ binomial(p)

 

     }

     

     ...

     However, I kept getting issues with the sampling processes. The model would either

     1) The sampler cannot initialize and return no samples. I got message saying things like: 

     

"Rejecting initial value:"                                                                                                   " Error evaluating the log probability at the initial value."                                                 "Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:"                        

"validate transformed params: prob[k0__] is 1.38673, but must be less than or equal to 1"                                               

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

"but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."   

     I am a bit confused here. I first thought it is an underflow issue, but I thought that would just give me an NAN? Given that the equation is by nature bounded, I am confused how it can produce prob that is larger than 1 or smaller than 0. Is it because my priors are too generous?

     2) When the sampler did get initialized. I sometimes will still get warnings like 

    "1: There were 3000 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.

2: Examine the pairs() plot to diagnose sampling problems" 

     

     When I increased the delta to .9, sometimes it solves the issue and all chains converge some reasonable number that makes intuitive sense. However, sometimes one or more of the chains still get stuck at the constrained boundaries and does not change much once it gets there nor converge to the same place where other chains converge. 

      Has anyone run into similar issues? And if so, what are the ways to solve the issue? When the sampler did move out of the boundaries, they do converge and always converge on numbers that are meaningful and interpretable. So I had the intuition that it could probably be solve by some re-parameterization and correspondingly adopt different sets of priors. 

   

      Thank you so much for your help! 

Best,

Shaorong

    

    

[Daniel Lee]

Hi Shaorong,

It really helps to provide the whole program. The error you're seeing isn't included in the snippet you sent, so it's much harder helping you. 

We try to be gelpful with the error messages as much as possible. There are still a couple of confusing parts, but it would help if you read the message. I'm try to help out a little below 

     1) The sampler cannot initialize and return no samples. I got message saying things like: 

     

"Rejecting initial value:"                                                                                                   " Error evaluating the log probability at the initial value."                                                 "Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:"                        

"validate transformed params: prob[k0__] is 1.38673, but must be less than or equal to 1"              

Here it's telling you that within the "transformed parameters" block the validation failed.

One particular value isn't obeying the constraints you set. 

Here, we can clearly see this message doesn't belong to the core snippet provided. The message is showing that the variable prob is actually an array or vector (or row vector). In the snippet provided, it's declared as a scalar.

I think you need to look at the definition of prob and you'll need to construct it such that the values satisfy the constraints for all valid values of whatever terms go into it.

                                

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

"but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."   

     I am a bit confused here. I first thought it is an underflow issue, but I thought that would just give me an NAN? Given that the equation is by nature bounded, I am confused how it can produce prob that is larger than 1 or smaller than 0. Is it because my priors are too generous?

It's not underflow. 

The user is responsible for ensuring transformed variables actually satisfy the constraints. You don't have to worry with parameters. 

The key is to check the equation you're using. If it were me, I'd print all the values before assigning to see what's causing the value to go past 1. 

     2) When the sampler did get initialized. I sometimes will still get warnings like 

    "1: There were 3000 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.

2: Examine the pairs() plot to diagnose sampling problems" 

     

     When I increased the delta to .9, sometimes it solves the issue and all chains converge some reasonable number that makes intuitive sense. However, sometimes one or more of the chains still get stuck at the constrained boundaries and does not change much once it gets there nor converge to the same place where other chains converge. 

I think this is all secondary. You'll need to fix your model first. 

Daniel

      Has anyone run into similar issues? And if so, what are the ways to solve the issue? When the sampler did move out of the boundaries, they do converge and always converge on numbers that are meaningful and interpretable. So I had the intuition that it could probably be solve by some re-parameterization and correspondingly adopt different sets of priors. 

   

      Thank you so much for your help! 

Best,

Shaorong

    

    

--
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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Shaorong Yan]

Hi Daniel, 

  

      Thanks a lot for your reply! I'll first deal with the issues you pointed out and see if that solves the problem. 

Best,

Shaorong

在 2017年4月21日星期五 UTC-4下午6:05:43，Daniel Lee写道：

[Shaorong Yan]

Hi Daniel, 

    So this is the model I specified for the model: 

"data {

    int s ;                         // Number of subjects

    int t ;                         // Number of trials per subject

    int N ;                       // total number of data points  

    int resp[N] ;                // number of total responses

    int subj[N] ;                // Subject number

    real trial[N] ;               // Trial number 

    real vot[N] ;                  //  VOT of each trial

      }

transformed data {

  

}

parameters {

    real a ;                        // Intercept

    real b ;                        // Slope

    real<lower = 0, upper = 0.1> lambda ;   // Lapse rate

    real k_a ;                              // Intercept change

    real k_b ;                              // Slope change

    real k_l ;                              // Lapse change

} 

transformed parameters {

   real<lower = 0, upper = 1> prob[N]  ;      // Probability of choosing b for each trial

   real k_a_sum[N]   ;

   real k_b_sum[N]   ;

   real<lower = 0, upper = 0.1> k_l_sum[N]   ;

  for(i in 1: N){

   k_a_sum[i] = a  + k_a * trial[i]/t;

   k_b_sum[i] = b  +  k_b * trial[i]/t;

   k_l_sum[i] = lambda  +  k_l * trial[i]/t;

   prob[i] = (1 - k_l_sum[i]) / (1 + exp((k_a_sum[i]  - vot[i]) / k_b_sum[i])) + k_l_sum[i] * 0.5 ;

   }

  

}  

  

  

model {

// priors

a ~  normal(20,  100)  ;

b ~  normal(0,  100)   ;

lambda ~  uniform(0, .1)   ;

k_a ~  normal(0,  1)    ;

k_b ~  normal(0,  1)   ;

k_l ~  normal(0,  .1)   ;

// model evaluation

resp ~ binomial(1, prob) ;

  

}

generated quantities {

real Intercept[t] ;

real Slope[t] ;

real Lapse[t] ;

for(i in 1: t){

   Intercept[i] = a  + k_a * trial[i]/t ;

   Slope[i] = b  + k_b * trial[i]/t ;

   Lapse[i] = lambda + k_l * trial[i]/t;

   }

}

"

Given how prob is specified, it should be within 0 and 1 as long as k_l_sum is larger than zero and smaller than 1. But I still sometimes get the error message I attached in my previous mail: 

""Rejecting initial value:"                                                                                        " Error evaluating the log probability at the initial value."                                     "Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:"                        

"validate transformed params: prob[k0__] is 1.38673, but must be less than or equal to 1"   "

     I got the impression that I need to change the parameterization or change the initial values. Is that where the problem is? 

     And if that's what I need to do, is there any recommended principled way to do so that works best when modeling with rstan?

     

      Thank you so much! 

Best,

Shaorong 

在 2017年4月21日星期五 UTC-4下午6:05:43，Daniel Lee写道：

Hi Shaorong,

[Daniel Lee]

The easiest way to figure it this out is to put a print statement with all the variables (lambda, k_l, trial[I], and t) before the assignment. 

I can't tell where it's going to fail by eye. 

You shave to find when that is violated. Perhaps the data isn't right 

But I still sometimes get the error message I attached in my previous mail: 

""Rejecting initial value:"                                                                                        " Error evaluating the log probability at the initial value."                                     "Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:"                        

"validate transformed params: prob[k0__] is 1.38673, but must be less than or equal to 1"   "

     I got the impression that I need to change the parameterization or change the initial values. Is that where the problem is? 

No. I think it'll be an error with the data or a programming error. 

     And if that's what I need to do, is there any recommended principled way to do so that works best when modeling with rstan?

This is basically debugging (holds for deterministic things too). Locate the first place where assumptions are violated. Fix, repeat. 

[Shaorong Yan]

I checked the model and it looks alright. 

I noticed that at the end of error message it says:

" Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model."                                 

[600] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                      

error occurred during calling the sampler; sampling not done

”

That’s why I thought it is a initialization issue. I guess I’ll check my model again and read the relevant parts of the stan-reference document to try figure this out. 

Thanks a lot for your help!

Shaorong

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

[Bob Carpenter]

The most likely failure here looks to me to be the
constraint you put on k_l_sum --- given that it's defined
to have lower=0 and upper=0.1 and has a term that is
unconstrained, you aren't going to have much luck initializing.
You want to make sure that everything is well defined for
every set of parameters that satisfy the declared constraints.

As Daniel says, putting print() [or reject()] statements in will let
you verify where things are going wrong.

You can speed this up a lot by bracketing (trial[i] / t) because otherwise
arithmetic is left associative and you incur a couple autodiff operations
multiplying parameters by constraints --- even better save trial[i] / t
in transformed parameters.

- Bob

[Shaorong Yan]

k. Thank you so much for the advice! I am now initializing in more constrained space so that k_l_sum does not go out of bound all the time.
I’ll look into the bracketing issue you pointed out too.

Thank you so much for your help!


Best,
Shaorong

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

[Bob Carpenter]

That's still going to be problematic.

Stan expects that if a set of parameter values satisfies
the declared constraints then it has a non-zero density
(i.e., a finite log density).

What you will be doing is chopping off part of the sample
space. This can be problematic for normalization and it
also messes with the posterior geometry (providing hard
constraints) through which Stan won't be able to sample
effectively.

We almost always recommend soft constraints on parameters
rather than trying to put hard constraints on derived
quantities. What you're doing doesn't have a straight
generative story from parameters to data, so it winds up
being problematic both statistically and computationally.

- Bob

[Shaorong Yan]

I see. So the right way is to parameterize in a way that the constrained parameters are drawn from a bounded distribution rather than truncated distribution. Is that what you mean?

Best,
Shaorong

[Bob Carpenter]

You have parameters:

>>>>> real<lower = 0, upper = 0.1> lambda ; // Lapse rate

>>>>> real k_l ;

and transformed parameters

>>>> real<lower = 0, upper = 0.1> k_l_sum[N] ;

>>>>
>>>> k_l_sum[i] = lambda + k_l * trial[i]/t;

If you make everything a vector you can do this more
efficiently in one line as

vector[N] k_l_sum = lambda + k_l * (trial / t);

But the real problem is that because k_l isn't bounded
this can run off of the constraint bounds.

You can put hard lower bounds on k_l and lambda if
they both have to be positive. Then you can put
softer upper bounds through a prior. If the fit
then runs past where you think it should, it indicates
a problem with the model's fit to the data. Adding
a hard upper bound in that constraint is just going to
cause problems as the probability mass tries to pile up
at the boundary.

- Bob

[Shaorong Yan]

I see what you mean now. Thanks a lot for answering my questions!

Best,
Shaorong