Flat Trace Plot

[Daniel Berman]

Hi All,


I noticed that upon repeatedly running a model on the same data set, I was getting massively larger variations in the parameter estimates. I decided to look at the values that it was generating at each step in R by calling the stored parameter on the saved data (hierarchical$parameter). When I looked at this, I realized that all the values in a single chain were the same.
What is going on here? Why is it doing this. Here is my code:

data{
  int<lower=1> N;
 
  int items_n_pos;    #integer indicating where the positive data starts
  int items_n_neg;    #integer indicating where the negative data starts
  int items_n_mix;    #integer indicating where the mixed values data starts
  int rawdata[items_n_mix]; #1s and 0s indicating subject's choice. Rawdata equals 1 means the subject choose option B.

  #The prospect data should contain four columns:
  #Column 1: Value1
  #Column 2: probability of Value1
  #Column 3: Value2
  #Column 4: probability of Value4
  matrix[items_n_pos,4] prospects_a_pos;  #Option A, positive data.
  matrix[items_n_neg-items_n_pos,4] prospects_a_neg; #Option A, negative data
  matrix[items_n_mix-items_n_neg,4] prospects_a_mix; #Option A, mixed data
  matrix[items_n_pos,4] prospects_b_pos;  #Option B, positive data
  matrix[items_n_neg-items_n_pos,4] prospects_b_neg; #Option B, negative data
  matrix[items_n_mix-items_n_neg,4] prospects_b_mix; #Option B, mixed data
}

#Transformed data breaks the prospects into individual vectors so it's easier to manipulate.
transformed data{
  vector[items_n_pos] x_a_pos;
  vector[items_n_pos] y_a_pos;
  vector[items_n_pos] px_a_pos;
  vector[items_n_pos] py_a_pos;

  vector[items_n_pos] x_b_pos;
  vector[items_n_pos] y_b_pos;
  vector[items_n_pos] px_b_pos;
  vector[items_n_pos] py_b_pos;

  vector[items_n_neg-items_n_pos] x_a_neg;
  vector[items_n_neg-items_n_pos] y_a_neg;
  vector[items_n_neg-items_n_pos] px_a_neg;
  vector[items_n_neg-items_n_pos] py_a_neg;
 
  vector[items_n_neg-items_n_pos] x_b_neg;
  vector[items_n_neg-items_n_pos] y_b_neg;
  vector[items_n_neg-items_n_pos] px_b_neg;
  vector[items_n_neg-items_n_pos] py_b_neg;

  vector[items_n_mix-items_n_neg] x_a_mix;
  vector[items_n_mix-items_n_neg] y_a_mix;
  vector[items_n_mix-items_n_neg] px_a_mix;
  vector[items_n_mix-items_n_neg] py_a_mix;

  vector[items_n_mix-items_n_neg] x_b_mix;
  vector[items_n_mix-items_n_neg] y_b_mix;
  vector[items_n_mix-items_n_neg] px_b_mix;
  vector[items_n_mix-items_n_neg] py_b_mix;
 
  #positive data
  for (ii in 1:items_n_pos){
    x_a_pos[ii]<-prospects_a_pos[ii,1];
    y_a_pos[ii]<-prospects_a_pos[ii,3];
    px_a_pos[ii]<-prospects_a_pos[ii,2];
    py_a_pos[ii]<-prospects_a_pos[ii,4];
   
    x_b_pos[ii]<-prospects_b_pos[ii,1];
    y_b_pos[ii]<-prospects_b_pos[ii,3];
    px_b_pos[ii]<-prospects_b_pos[ii,2];
    py_b_pos[ii]<-prospects_b_pos[ii,4];
  }

  #negative data
  for (ii in 1:(items_n_neg-items_n_pos)){
    x_a_neg[ii]<-prospects_a_neg[ii,1];
    y_a_neg[ii]<-prospects_a_neg[ii,3];
    px_a_neg[ii]<-prospects_a_neg[ii,2];
    py_a_neg[ii]<-prospects_a_neg[ii,4];
   
    x_b_neg[ii]<-prospects_b_neg[ii,1];
    y_b_neg[ii]<-prospects_b_neg[ii,3];
    px_b_neg[ii]<-prospects_b_neg[ii,2];
    py_b_neg[ii]<-prospects_b_neg[ii,4];
  }

  #mixed data
  for (ii in 1:(items_n_mix-items_n_neg)){
    x_a_mix[ii]<-prospects_a_mix[ii,1];
    y_a_mix[ii]<-prospects_a_mix[ii,3];
    px_a_mix[ii]<-prospects_a_mix[ii,2];
    py_a_mix[ii]<-prospects_a_mix[ii,4];
 
    x_b_mix[ii]<-prospects_b_mix[ii,1];
    y_b_mix[ii]<-prospects_b_mix[ii,3];
    px_b_mix[ii]<-prospects_b_mix[ii,2];
    py_b_mix[ii]<-prospects_b_mix[ii,4];
  }
}

#parameters of interest
parameters{
  real<lower=0,upper=1> alpha;
  real<lower=0,upper=1> beta;
  real<lower=0,upper=1> delta;
  real<lower=0,upper=1> gamma;
  real<lower=1> lambda;
  real<lower=0> luce;
}


model {
#Declaring Variables
  vector[items_n_pos] v_x_a_pos;
  vector[items_n_pos] v_y_a_pos;
  vector[items_n_pos] v_x_b_pos;
  vector[items_n_pos] v_y_b_pos;

  vector[items_n_neg-items_n_pos] v_x_a_neg;
  vector[items_n_neg-items_n_pos] v_y_a_neg;
  vector[items_n_neg-items_n_pos] v_x_b_neg;
  vector[items_n_neg-items_n_pos] v_y_b_neg;
 
  vector[items_n_mix-items_n_neg] v_x_a_mix;
  vector[items_n_mix-items_n_neg] v_y_a_mix;
  vector[items_n_mix-items_n_neg] v_x_b_mix;
  vector[items_n_mix-items_n_neg] v_y_b_mix;
 
  vector[items_n_pos] w_x_a_pos;
  vector[items_n_pos] w_x_b_pos;
  vector[items_n_pos] w_y_a_pos;
  vector[items_n_pos] w_y_b_pos;
 
  vector[items_n_neg-items_n_pos] w_x_a_neg;
  vector[items_n_neg-items_n_pos] w_x_b_neg;
  vector[items_n_neg-items_n_pos] w_y_a_neg;
  vector[items_n_neg-items_n_pos] w_y_b_neg;

  vector[items_n_mix-items_n_neg] w_x_a_mix;
  vector[items_n_mix-items_n_neg] w_x_b_mix;
  vector[items_n_mix-items_n_neg] w_y_a_mix;
  vector[items_n_mix-items_n_neg] w_y_b_mix;
 
  vector[items_n_pos] vf_a_pos;
  vector[items_n_pos] vf_b_pos;

  vector[items_n_neg-items_n_pos] vf_a_neg;
  vector[items_n_neg-items_n_pos] vf_b_neg;

  vector[items_n_mix-items_n_neg] vf_a_mix;
  vector[items_n_mix-items_n_neg] vf_b_mix;
 
  vector[items_n_mix] vf_a;
  vector[items_n_mix] vf_b;
 
  vector[items_n_mix] binval;
 
#Actual code
  #Loops through positive data.
  for (ii in 1:items_n_pos){
    #Applies value function to data
    v_x_a_pos[ii]<-pow(x_a_pos[ii],alpha);
    v_y_a_pos[ii]<-pow(y_a_pos[ii],alpha);
    v_x_b_pos[ii]<-pow(x_b_pos[ii],alpha);
    v_y_b_pos[ii]<-pow(y_b_pos[ii],alpha);
   
    #Applied probability weighting function to data
    w_x_a_pos[ii]<-pow(px_a_pos[ii],gamma)/pow((pow(px_a_pos[ii],gamma)+pow((1-px_a_pos[ii]),gamma)),1/gamma);
    w_x_b_pos[ii]<-pow(px_b_pos[ii],gamma)/pow((pow(px_b_pos[ii],gamma)+pow((1-px_b_pos[ii]),gamma)),1/gamma);
    w_y_a_pos[ii]<-pow(py_a_pos[ii],gamma)/pow((pow(py_a_pos[ii],gamma)+pow((1-py_a_pos[ii]),gamma)),1/gamma);
    w_y_b_pos[ii]<-pow(py_b_pos[ii],gamma)/pow((pow(py_b_pos[ii],gamma)+pow((1-py_b_pos[ii]),gamma)),1/gamma);
   
    #combines the weighted probabilites with the value function to form weighted prospects
    vf_a_pos[ii]<- v_x_a_pos[ii]*w_x_a_pos[ii]+v_y_a_pos[ii]*w_y_a_pos[ii];
    vf_b_pos[ii]<- v_x_b_pos[ii]*w_x_b_pos[ii]+v_y_b_pos[ii]*w_y_b_pos[ii];
  }
 
  #Loops through negative data
  for (ii in 1:(items_n_neg-items_n_pos)){
   
    #Applies value function to data
    v_x_a_neg[ii]<- -lambda*pow(fabs(x_a_neg[ii]),beta);
    v_y_a_neg[ii]<- -lambda*pow(fabs(y_a_neg[ii]),beta);
    v_x_b_neg[ii]<- -lambda*pow(fabs(x_b_neg[ii]),beta);
    v_y_b_neg[ii]<- -lambda*pow(fabs(y_b_neg[ii]),beta);
   
    #Applied probability weighting function to data
    w_x_a_neg[ii]<-pow(px_a_neg[ii],delta)/pow((pow(px_a_neg[ii],delta)+pow((1-px_a_neg[ii]),delta)),1/delta);
    w_x_b_neg[ii]<-pow(px_b_neg[ii],delta)/pow((pow(px_b_neg[ii],delta)+pow((1-px_b_neg[ii]),delta)),1/delta);
    w_y_a_neg[ii]<-pow(py_a_neg[ii],delta)/pow((pow(py_a_neg[ii],delta)+pow((1-py_a_neg[ii]),delta)),1/delta);
    w_y_b_neg[ii]<-pow(py_b_neg[ii],delta)/pow((pow(py_b_neg[ii],delta)+pow((1-py_b_neg[ii]),delta)),1/delta);

    #combines the weighted probabilites with the value function to form weighted prospects
    vf_a_neg[ii]<- v_x_a_neg[ii]*w_x_a_neg[ii]+v_y_a_neg[ii]*w_y_a_neg[ii];
    vf_b_neg[ii]<- v_x_b_neg[ii]*w_x_b_neg[ii]+v_y_b_neg[ii]*w_y_b_neg[ii];
  }

  #Loops through mixed data
  for (ii in 1:(items_n_mix-items_n_neg)){

    #Applies value function to data
    v_x_a_mix[ii]<-pow(x_a_mix[ii],alpha);
    v_y_a_mix[ii]<- -lambda*pow(fabs(y_a_mix[ii]),beta);
    v_x_b_mix[ii]<-pow(x_b_mix[ii],alpha);
    v_y_b_mix[ii]<- -lambda*pow(fabs(y_b_mix[ii]),beta);
   
    #Applied probability weighting function to data
    w_x_a_mix[ii]<-pow(px_a_mix[ii],gamma)/pow((pow(px_a_mix[ii],gamma)+pow((1-px_a_mix[ii]),gamma)),1/gamma);
    w_x_b_mix[ii]<-pow(px_b_mix[ii],delta)/pow((pow(px_b_mix[ii],delta)+pow((1-px_b_mix[ii]),delta)),1/delta);
    w_y_a_mix[ii]<-pow(py_a_mix[ii],gamma)/pow((pow(py_a_mix[ii],gamma)+pow((1-py_a_mix[ii]),gamma)),1/gamma);
    w_y_b_mix[ii]<-pow(py_b_mix[ii],delta)/pow((pow(py_b_mix[ii],delta)+pow((1-py_b_mix[ii]),delta)),1/delta);
   
    #combines the weighted probabilites with the value function to form weighted prospects
    vf_a_mix[ii]<- v_x_a_mix[ii]*w_x_a_mix[ii]+v_y_a_mix[ii]*w_y_a_mix[ii];
    vf_b_mix[ii]<- v_x_b_mix[ii]*w_x_b_mix[ii]+v_y_b_mix[ii]*w_y_b_mix[ii];
  }
 
  #These three loops combine the weighted prospects into one variable.
  for (j in 1:items_n_pos){
    vf_a[j]<-vf_a_pos[j];
    vf_b[j]<-vf_b_pos[j];
  }
  for (j in 1:(items_n_neg-items_n_pos)){
    vf_a[j+items_n_pos]<-vf_a_neg[j];
    vf_b[j+items_n_pos]<-vf_b_neg[j];
  }
  for (j in 1:(items_n_mix-items_n_neg)){
      vf_a[j+items_n_neg]<-vf_a_mix[j];
    vf_b[j+items_n_neg]<-vf_b_mix[j];
  }

  #Applies the Luce Choice rule to find the probability of choosing option B over A.
  for (i in 1:items_n_mix){
    binval[i]<-exp(luce*vf_b[i])/(exp(luce*vf_a[i])+exp(luce*vf_b[i]));
  }

  #priors for the variables of interest
  alpha~uniform(0,1);
  beta~uniform(0,1);
  gamma~uniform(0,1);
  delta~uniform(0,1);
  lambda~uniform(1.01,5);
  luce~uniform(.01,1.5);

  #likelihood function. Takes rawdata and applies a bernouli distribution for the probability of picking option B over A.
  #Rawdata equals 1 means the subject choose option B.
  for (ii in 1:items_n_mix){
    rawdata[ii]~bernoulli(binval[ii]);
  }
}

[Bob Carpenter]

First, you really want to remove all that cut-and-paste and
replace with either or both of arrays and functions. It's very
hard to see what's going on in all that code.

MCMC can get stuck if it can't find a way to jump from where it's
at to a new parameter values with reasonable densities. So you
might try initializing to something reasonable.

Otherwise, you can look at the initialization that's causing problem
and run diagnostics and see what the gradients look like, which
can help diagnose (just be warned they're on the unconstrained scale
and are only indexed by number, not by name --- something we need
to improve on [naming, not being on the uncontsrained scale]).

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

[Daniel Berman]

At the moment I'm not defining initial values. If that's a significant problem I can address that. However, I'm not sure how I would run those diagnostics. How would you recommend doing that?

[Daniel Berman]

I checked with a number of initial values and I consistently get flat trace plots. Additionally, I added print statements after I defined the Priors for the parameters and I can see it change.

How do i determine the gradient? My function is not continuous and has different parameters for positive and negative values.

I apologize for my poor programming style. I am new to programming and working on improving the structure of the code. However, the problem itself is complicated and requires a lot of cut-and-pasting.

[Bob Carpenter]

Overall, my recommendation would be to start with an easier model
Well one problem you're going to have is that your parameter declarations

real<lower=1> lambda;
real<lower=0> luce;

are not consistent with their support:

lambda~uniform(1.01,5);
luce~uniform(.01,1.5);

If you really insist on those uniform priors, which we don't
recommend (see the manual chapter on regression for our recommended
priors and motivation and references), then you need to declare the
variables as

real<lower=1.01, upper=5> lambda;
real<lower=0.01, upper=1.5> luce;

Rather than writing this:

> w_x_a_pos[ii]<-pow(px_a_pos[ii],gamma)/pow((pow(px_a_pos[ii],gamma)+pow((1-px_a_pos[ii]),gamma)),1/gamma);
> w_x_b_pos[ii]<-pow(px_b_pos[ii],gamma)/pow((pow(px_b_pos[ii],gamma)+pow((1-px_b_pos[ii]),gamma)),1/gamma);
> w_y_a_pos[ii]<-pow(py_a_pos[ii],gamma)/pow((pow(py_a_pos[ii],gamma)+pow((1-py_a_pos[ii]),gamma)),1/gamma);
> w_y_b_pos[ii]<-pow(py_b_pos[ii],gamma)/pow((pow(py_b_pos[ii],gamma)+pow((1-py_b_pos[ii]),gamma)),1/gamma);

use

functions {
real foo(real p, real gamma) {
return pow(p, gamma) / pow(pow(p, gamma) + pow(1 - p, gamma), inv(gamma));
}
}

and replace all that cut-and-paste with

w_x_a_pos[ii] <- foo(px_a_pos[ii], gamma);
...

Of course, choose a good name for foo. I have no idea what this model
is doing, though --- power functions can be very unstable if the exponents
are too large or too small.

I'd also recommend replacing all of this cut-and-paste

> vector[items_n_neg-items_n_pos] x_b_neg;
> vector[items_n_neg-items_n_pos] y_b_neg;
> vector[items_n_neg-items_n_pos] px_b_neg;
vector[items_n_neg-items_n_pos] py_b_neg;

with further array indexing, such as

vector[items_n_neg - items_n_pos] neg[2,2];

where neg[1,1] = x_a_neg, neg[1,2] = x_b_neg,
neg[2,1] = y_a_neg, etc. It'll let you replace all
the cut-and-paste with loops and should make the code
much easier to read.


Are you sure that this:

> binval[i]<-exp(luce*vf_b[i])/(exp(luce*vf_a[i])+exp(luce*vf_b[i]));

is going to be identified? Usually a binary logistic like this
is coded using only a single predictor, whereas you have three degrees
of freedom here with luce, vf_b and vf_a.

- Bob





> On Apr 28, 2015, at 5:48 AM, Daniel Berman <daniel....@gmail.com> wrote:
>