no matches for normal(real, real)

[Christine C]

I am attempting to build a model describing economic choices.  I would like to model noise in subjects' estimates of probability and value.  It is known that people exhibit scalar variability, i.e. the standard deviation in their estimate of magnitude grows linearly with the magnitude.  So, I would like to model the subjects' perception of reward amount/probability, accounting for the linear noise scaling.  Something like this:

perceived_rprob[trial_i] <- normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset);

However, I get the following error: no matches for normal(real, real).  Function normal not found.  So I must be miscasting the variables? I have also tried the line below instead, but from what I can tell from this group and the documentation, the ~ increments the log posterior, and does not mean "drawn from."  This throws an error that the dimension of the right side of the statement is 0.

perceived_rprob[trial_i] ~ normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset); 

I am using rstan; my code is below.  I apologize if I am doing something super dumb, any help is very much appreciated.  Thank you so much!!!

model_code <- "data {

    int nTrials;

    int choices[nTrials];  //values of 1 (chose left) or 2 (chose right)    

    real<lower=0, upper=1> leftprob[nTrials]; 

    real<lower=0, upper=1> leftvol[nTrials];  //rewards are normalized by max reward

    real<lower=0, upper=1> rightprob[nTrials];

    real<lower=0, upper=1> rightvol[nTrials];  //rewards are normalized by max reward

}

transformed data {

// this executes once vs on each sample in the model block, so better to transform here and pay for the conditional once

int<lower=0,upper=1> choices_01[nTrials];

for (i in 1:nTrials){

choices_01[i] <- (choices[i]==2);

}

}

parameters {

    real rho; //exponential of utility function

    real alpha;  //parameter for probability weighting function

    real<lower=0> lambda; //softmax parameter

    real<lower=0> prob_m; //slope of scalar variability relationship

    real<lower=0> prob_offset; //offset of scalar variability relationship

    real<lower=0> vol_m; 

    real<lower=0> vol_offset;

}

model {

    real Rval[nTrials];

    real Lval[nTrials];

    real perceived_rprob[nTrials];

    real perceived_lprob[nTrials];

    real perceived_rvol[nTrials];

    real perceived_lvol[nTrials];    

    rho ~ cauchy(0, 2.5);

    for (trial_i in 1:nTrials) {    

    perceived_rprob[trial_i] <- normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset);

    perceived_lprob[trial_i] <- normal(leftprob[trial_i], leftprob[trial_i]*prob_m+prob_offset);

    perceived_rvol[trial_i] <- normal(rightvol[trial_i], rightvol[trial_i]*vol_m+vol_offset);

    perceived_lvol[trial_i] <- normal(leftvol[trial_i], leftvol[trial_i]*vol_m+vol_offset);  

    Rval <- exp(-(-log(perceived_rprob[trial_i]))^alpha)*(perceived_rvol[trial_i]^rho);

    Lval <- exp(-(-log(perceived_lprob[trial_i]))^alpha)*(perceived_lvol[trial_i]^rho);

    choices_01[trial_i] ~ bernoulli_logit(lambda*(Rval-Lval));

    }

}"

[Ben Lambert]

Thanks for your message. I think you are right - essentially  'rightprob[trial_i]*prob_m+prob_offset' is of type real because you have:

real<lower=0,upper=1> * real<lower=0> + real<lower=0>, which gets promoted to a real, rather than real<lower=0>, as you might expect. 

Am not sure this'll work - but it is worth trying. Define rightprob as a real, rather than real<lower=0,upper=1>; since you're passing the data, you will need to check that it satisfies your constraints before the data goes into Stan. Based on my own models (which have heteroscedasticity) I think this may work. Let me know.

Alternatively, you could try defining a transformed parameter:

transformed parameter{

    real<lower=0> aSigma;

    aSigma <- rightprob[trial_i]*prob_m+prob_offset;

}

And using aSigma instead for your normal standard deviation? You may suffer the same issue again: Stan may not like trying to coerce the RHS to a type real<lower=0>.

If neither of these work, my apologies! (If you attach a R data file with all the necessary data variables, then this may make debugging quicker.)

Best,

Ben

[Christine C]

Hi Ben,

Thanks so much for your reply!  I tried both of your suggestions, but received the same error message each time: no matches for normal(real, real).  I am attaching an RDS file with simulated data.  Please let me know if you have any other ideas, suggestions, or would prefer the data in a different format.  (Disclaimer: I have just started using R so I can use the diagnostics in rstan, so apologies if the .RDS format is weird?)

Thanks!

Christine

[Ben Lambert]

Sorry - it was a silly mistake on my part. The issue is the following: if you use an assignment operator '<-', you need to use normal_rng(). Are you trying to generate fake data, then sample from it? If so, then you want to generate fake data using normal_rng because you want Stan to generate (pseudo)-random samples from a normal. 

If this is the case, may I suggest the following: you generate the data in one Stan file, using the 'generated quantities' block as below (obviously you need to pass it data, etc.) To do this use the algorithm='Fixed_param' option when you call stan in R; you won't have any parameters in your model so you need to tell Stan not to look for them.

generated quantities{

 real perceived_rprob[nTrials];

  ....

  for(trial_i in 1:nTrials){

       perceived_rprob[trial_i] <- normal_rng(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset);

        ...   

   }

}

You then collect the samples for your various objects, and then fit your model using (in a separate Stan file):

model{

    for (trial_i in 1:NTrials){

        ....

        choices_01[trial_i] ~ bernoulli_logit(lambda*(Rval-Lval));

    }

}

There may be a way to fit it all in one Stan file, but think doing the above makes life easier.

Let me know if the above is what you want, and makes sense.

Best,

Ben

[Christine C]

Thanks Ben!  One thing that I am confused about is that prob_m and prob_offset are model parameters that I am trying to fit with stan (these define the slope and offset of the linear relationship between signal and noise).  If I understand your suggestion, you are saying I should generate a bunch of fake data first, but that would require me to specify a particular prob_m and prob_offset for that generated dataset.  That is problematic, no?

I was not thinking of this as a hierarchical model, but I wonder if I should be thinking of these as hyperparameters?  

Please let me know if I am just confused (which is likely).  Thanks for your time and feedback!

Christine

[Ben Lambert]

Hi Christine,

Thanks for your message. I am a bit confused, what is  real perceived_rprob[nTrials]? Is this fake data, or something else?

If it isn't fake data, then where is it coming from? I'm a bit lost as to why you would then fit a binomial to generated data?

Best,

Ben

[Christine C]

On a given trial, a subject is presented with sensory cues that tell them the probability of reward on the right side (real rprob). Because of the particular sensory cues we are using, there is some noise in their estimate of real rprob, and so on a given trial, their estimate of real rprob (perceived_rprob) can be described as a random draw from a normal distribution with mean=real rprob, and with st. dev.=real rprob*m_prob +prob_offset.  

I have made a model in stan that accounts for the trial by trial choice data of subjects based on real probabilities and real reward amounts.  But we know that there is noise in subjects' estimate of probabilities and reward amounts, and that the noise scales in a particular way (i.e., the st. dev. increases linearly with the magnitude), based on previous psychophysical data.  I was hoping to have a model that could account for the trial by trial choice data while simultaneously accounting for the perceptual noise, which we know exhibits a linear relationship.

Does that make sense?  If this is impossible in stan, or completely misguided, please let me know.  I previously tried to do this using gradient descent, but introducing rngs is problematic for such a model.  That was my main motivation for building the model in stan instead.

Thanks again,

Christine

[Ben Lambert]

Sorry, just trying to understand...do you have participants' estimated rprob?

[Christine C]

No, that is one of the things I am trying to estimate with the model.  I just have their choices (chose left or chose right).

[Ben Lambert]

Ok, thanks. But how come you are passing rightprob in the data? What's the difference between 'rightprob' and 'perceivedprob'? Sorry if am being a bit slow here!

[Christine C]

You are not being slow!  I am passing the trial information into the data (what were the real reward probabilities and amounts for the left and right sides).  Perceived_rprob is some filtered version of rightprob that I don't have access to but would like to infer.  More specifically, I would like to model the subjects' percept of the probability on each side as coming from a normal distribution with mean corresponding to the real probability, and with st. dev. as a linear function of the mean (or real rprob).  The subject then compares the perceived relative value of the two options (Rval-Lval) and makes a choice based on that.  

So, rightprob is the actual probability of receiving reward if the subject chooses right.  I know that value because I generated it.  However, I convey it to the subject with some sensory cues that the subject may perceive with some noise.  Perceivedprob is a latent variable that the subject uses to make his/her choice.  Based on previous work, I know how to parameterize the relationship between rightprob and perceivedprob, namely, that the st. dev. in perceivedprob grows linearly with rightprob.

The more I think about it, the more I think the slope/offset parameters that relate rightprob to perceivedprob should be hyperparameters in a hierarchical model.  Let me know if that makes more sense, and I apologize for all of the confusion.

[Ben Lambert]

Sorry just reread your above answer. Is rightprob the 'true' probability?

If so, then I would think you want to do something like the following:

rightprob[trial_i] ~ normal(perceived_rprob[trial_i] , perceived_rprob[trial_i]*prob_m+prob_offset);

This basically does what you want, but will allow estimates of perceived_rprob to be obtained (in which case you want to define it in the parameters block.) Would that work?

Best,

Ben

[Ben Lambert]

Ok, see my above answer - not sure about hyper-parameters etc. here. But the above will probably work as you want it to...

[Christine C]

Yes, rightprob is the "true" probability.

I think that is very similar to what I initially tried to implement (I assumed that the perceivedprob is drawn from a normal with mean = real rprob).  The model parameters I previously defined were the slope and offset of the line relating the true probability to the perceived probability.  Is it appropriate to also include perceivedprob as a model parameter?  Because that would be a unique value for the left and right side for each trial.  I really just want to find the slope and offset relating real to perceived probability, not the perceivedprob on each trial (that would be an insane number of model parameters!)

Also, I think I still will get the error message that there is no match for normal(real, real).

Thank you for all of your input.

Christine

[Ben Lambert]

Hi Christine,

You won't get the error message if you use the '~' statement rather than '<-', I hope...

How many parameters/trials are we talking about? I often fit models with tens of thousands of parameters... You said earlier that you wanted to infer Perceived_rprob? Do you mean, it's more of a nuisance parameter? (Necessary for the resultant inferences, but not useful in itself.) If that is the case, then you are right to define it at the top of the model block, as opposed to in the parameters block. I do think that hierarchical parameters would work nicely here, yes.

I have just tried the following, and it runs. However, it looks like there are issues with the model still, but could be due to your generated data...(Are the probabilities supposed to be 1s and 0s mostly?) Also I think that since you are dealing with a variable constrained to lie between 0 and 1, you need to either:

- Take the logit of your rightprob variables etc. to make these unconstrained.

- Choose a different sampling model, for example, a beta distribution.

model {

    real Rval[nTrials];

    real Lval[nTrials];

    real perceived_rprob[nTrials];

    real perceived_lprob[nTrials];

    real perceived_rvol[nTrials];

    real perceived_lvol[nTrials];    

    rho ~ cauchy(0, 2.5);

    for (trial_i in 1:nTrials) {    

    rightprob[trial_i] ~ normal(perceived_rprob[trial_i],perceived_rprob[trial_i]*prob_m+prob_offset);

    leftprob[trial_i] ~ normal(perceived_lprob[trial_i], perceived_lprob[trial_i]*prob_m+prob_offset);

    rightvol[trial_i] ~ normal(perceived_rvol[trial_i], perceived_rvol[trial_i]*vol_m+vol_offset);

    leftvol[trial_i] ~ normal(perceived_lvol[trial_i], perceived_lvol[trial_i]*vol_m+vol_offset);  

    Rval[trial_i] <- exp(-(-log(perceived_rprob[trial_i]))^alpha)*(perceived_rvol[trial_i]^rho);

    Lval[trial_i] <- exp(-(-log(perceived_lprob[trial_i]))^alpha)*(perceived_lvol[trial_i]^rho);

    choices_01[trial_i] ~ bernoulli_logit(lambda*(Rval[trial_i]-Lval[trial_i]));

    }

}

[Michael Betancourt]

Christine,

I believe this would be simplified if you think more generatively.  The idea is to

build your model to mimic the data generating process itself.  For example, in

a simple Binomial experiment we have some latent “true” probability of success

that we measure in a finite number of trials,

k ~ Binomial(N, p);

and we infer p from k and N.  This can be made much more complicated, with

different p for many different individuals or trials, or having p be a function of 

covariates (leading to logistic regression), etc.

I don’t fully understand the model that you are trying to build, but it sounds like

there are some data censoring issues.  For example, having just (k / N) instead

of k and N separately makes this drastically more complicated because you

can’t exactly recover the likelihood.  

Having some empirical p = k / N, call it p_hat, along with some estimate of the 

standard deviation also runs into problems because the Gaussian approximation 

is valid only for true p away from the 0 and 1 boundaries.  To avoid issues we

recommend treating such convolutions on the logit scale, for example

logit_obs_prop ~ normal(logit_latent_prob, m * log_latent_prob + offset).

That said, I would be very hesitant in this particular model as it is very hard to

derive from any reasonable generative model (I understand the psycophysical

literature might support it, but statistical misunderstandings are rampant in

many, many fields).  

So long story short, the kind of model you want (latent parameters that aren’t

directly observed with hierarchical structure) is very much appropriate for Stan

but at the moment you’re likely being limited by some of the assumptions in

your model.

It would help if you could simplify your model to just the essence of what

you’re trying to fit and then it would be much easier for us to comment and

make more explicit recommendations.

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

[Christine C]

Thank you Ben and Michael.  I apologize for the confusion, I have not done a good job of describing the model.

According to some economic theories, humans choose between two prospects, or gambles, based on the relative value of those prospects (in my case, Rval and Lval), where the value is the product of a nonlinear subjective function relating true probability and perceived probability (in my case, exp(-(-log(p))^alpha)), and a nonlinear function relating true reward amount and perceived reward amount (in my case, vol^rho).  My primary interest is in inferring the parameters alpha and rho from many thousands of choices between different gambles.  

However, in my experiment, I am using certain sensory cues to convey reward probability and reward amount to the subject (e.g., # light flashes conveys the true probability).  It is known that subjects' estimates of such stimuli can be noisy.  Specifically, the standard deviation in the subjects' estimate of # flashes grows linearly with the # flashes presented.

My concern is that if I don't account for this sensory noise, it could mess up my estimates of alpha and rho.  So I was hoping to compare the alpha and rho parameters before accounting for noise (that model already works in stan) and afterwards, to see if it is something I should worry about.

Does that make sense?  Please let me know if I should clarify anything, and thanks so much for your feedback!

[Michael Betancourt]

Thank you Ben and Michael.  I apologize for the confusion, I have not done a good job of describing the model.

No worries, just for future reference we can help much better

with short, simple examples that capture your problem.

According to some economic theories, humans choose between two prospects, or gambles, based on the relative value of those prospects (in my case, Rval and Lval), where the value is the product of a nonlinear subjective function relating true probability and perceived probability (in my case, exp(-(-log(p))^alpha)), and a nonlinear function relating true reward amount and perceived reward amount (in my case, vol^rho).  My primary interest is in inferring the parameters alpha and rho from many thousands of choices between different gambles.  

However, in my experiment, I am using certain sensory cues to convey reward probability and reward amount to the subject (e.g., # light flashes conveys the true probability).  It is known that subjects' estimates of such stimuli can be noisy.  Specifically, the standard deviation in the subjects' estimate of # flashes grows linearly with the # flashes presented.

My concern is that if I don't account for this sensory noise, it could mess up my estimates of alpha and rho.  So I was hoping to compare the alpha and rho parameters before accounting for noise (that model already works in stan) and afterwards, to see if it is something I should worry about.

Does that make sense?  Please let me know if I should clarify anything, and thanks so much for your feedback!

Yeah, that’s conceptually no problem for Stan, you just have

to be careful in how you build the model.  You start with

some latent p_true and then some observer-specific p_true_obs 

that gets corrected for the observer-dependent environment,

and then some data generated from each p_true_obs.  

p_true and all of the p_true_obs would be parameters that

you fit jointly.

The problem with working with probabilities is that they

have to stay within 0 and 1, so you can’t apply a correction

that gives p_true_obs outside of these ranges.  I _highly_

recommend reformulating your model in terms of logit(p),

which is unconstrained and then constraining back to

a probability with a logistic only when constructing the

measurement part of the likelihood.

So, as a very informal example, you might have something

like

logit_p_true ~ normal(0, 5); // covers the bulk of the probabilities

// alpha quantifies how easily logit_p_true is conveyed to each observer

// beta1 and beta1 quantifies how much the rewards bias the observer-specific probability

mu_logit_p_obs[i] = alpha * logit_p_true + beta1 * positive_reward_strength - beta2 * negative_reward_strength;

logit_p_obs[i] ~ normal(mu_logit_p_obs[i], sigma);

k[i] ~ binomial_logit(N[i], logit_p_obs[i]);

And you could make sigma depend on the latent probability,

too.  I’m sure the model would become very complicated

very quickly, but by starting with something simple and adding

complexity bit by bit you should be able to develop a robust,

principled model.

[Christine C]

Thank you so much for your help!  I apologize again for not explaining the model clearly enough.  

[Christine C]

Okay, I have been thinking about this and I have a few questions, if that's okay.

1. In terms of reformulating the model to be in terms of logit(p), I have a large number of trials in which p=1, and logit(p) = inf.  Do you have any suggestions about how I should handle those trials?  There are many of them, and I would prefer to not exclude them.

2. For k[i] ~ binomial_logit(N[i], logit_p_obs[i]), what is N[i]?  I couldn't find this in the manual/documentation.  Does it correspond to that particular trial number?  Or the total number of trials?

Thank you so much!

Christine

[Michael Betancourt]

>
> 1. In terms of reformulating the model to be in terms of logit(p), I have a large number of trials in which p=1, and logit(p) = inf. Do you have any suggestions about how I should handle those trials? There are many of them, and I would prefer to not exclude them.

But aren’t the p unobserved? If they’re parameters then unless the data are all 0’s and 1’s
then the latent p would never quite reach 0 or 1. logit(p), which can get very close to 0 and 1
would the be okay.

If there are some latent p that are known then you’ll have to use the nominal scale between 0
and 1. But then if you convolve those known p to get some other latent, unobserved p (as in
those individual belief p’s) then you’ll have to be careful to ensure that you stay within those
boundaries.

> 2. For k[i] ~ binomial_logit(N[i], logit_p_obs[i]), what is N[i]? I couldn't find this in the manual/documentation. Does it correspond to that particular trial number? Or the total number of trials?

Total number of trials. If you have N trials with probability of success p, then
the number of successes is distributed as

k ~ binomial(N, p);

Or

k ~ binomial_logit(N, logit(p));

[Christine C]

Thanks for your reply!  I apologize for the delay, I got derailed this week.  

The true p are observed, and are often 1. I was having trouble thinking about things in terms of logit(p), so instead I am treating latent_p = tanh(flash_estimate), to ensure that the latent p never exceeds one.   I am attempting to think more generatively about the model; my code is below.  The model compiles, but I get the following error message:

"Rejecting initial value:"                                                                                                              

"Error evaluating the log probability at the initial value."   

"Exception thrown at line 50: stan::math::normal_log: Random variable is nan, but must not be nan!"                                           

"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

Stan model '2dd970cfdc0651a12bf6907215f11fd1' does not contain samples.

Error in dimnames(x)[[3]] : subscript out of bounds

I don't understand why the random variable/which random variability is repeatedly nan.  If you have any insight, it would be very much appreciated!  Thanks so much!!!  

model_code <- "data {

    int nTrials;

    int choices[nTrials];  //values of 1 (chose left) or 2 (chose right)   

    real<lower=0, upper=1> leftprob[nTrials]; 

    real<lower=0, upper=1> leftvol[nTrials];  //rewards are normalized by max reward

    real<lower=0, upper=1> rightprob[nTrials];

    real<lower=0, upper=1> rightvol[nTrials];  //rewards are normalized by max reward

}

transformed data {

// this executes once vs on each sample in the model block, so better to transform here and pay for the conditional once

int<lower=0,upper=1> choices_01[nTrials];

for (i in 1:nTrials){

choices_01[i] <- (choices[i]==2);

}

}

parameters {

    real rho; //exponential of utility function

    real alpha;  //parameter for probability weighting function

    real<lower=0> lambda; //softmax parameter

    real<lower=0> prob_m; //slope of scalar variability relationship

    real<lower=0> prob_offset; //offset of scalar variability relationship

    real<lower=0> vol_m; 

    real<lower=0> vol_offset;

}

model {

    real Rval[nTrials];

    real Lval[nTrials];

    real rflashes[nTrials]; //sensory cues conveying probability

    real lflashes[nTrials];

    real rclicks[nTrials]; //sensory cues conveying reward amount

    real lclicks[nTrials];

    int latent_rprob[nTrials]; //latent estimate of probability

    int latent_lprob[nTrials];

    int latent_rvol[nTrials]; //latent estimate of reward amount

    int latent_lvol[nTrials];

    

    rho ~ cauchy(0, 2.5);

    for (trial_i in 1:nTrials) {

    rflashes[trial_i] ~ normal(rightprob[trial_i], rightprob[trial_i]*prob_m + prob_offset); //latent variable: number of flashes are estimated with gaussian noise

    lflashes[trial_i] ~ normal(leftprob[trial_i], leftprob[trial_i]*prob_m + prob_offset);    

    latent_rprob[trial_i] ~ binomial(nTrials, tanh(rflashes[trial_i])); //noisy estimate of flashes maps onto latent estimate of probability

    latent_lprob[trial_i] ~ binomial(nTrials, tanh(lflashes[trial_i]));    

    lclicks[trial_i] ~ normal(leftvol[trial_i], leftvol[trial_i]*vol_m + vol_offset);

    rclicks[trial_i] ~ normal(rightvol[trial_i], rightvol[trial_i]*vol_m + vol_offset);    

    latent_rvol[trial_i] ~ binomial(nTrials, tanh(rclicks[trial_i]));

    latent_lvol[trial_i] ~ binomial(nTrials, tanh(lclicks[trial_i]));       

    Rval[trial_i] <- exp(-(-log(latent_rprob[trial_i]))^alpha)*(latent_rvol[trial_i]^rho);

    Lval[trial_i] <- exp(-(-log(latent_lprob[trial_i]))^alpha)*(latent_lvol[trial_i]^rho);

        choices_01[trial_i] ~ bernoulli_logit(lambda*(Rval[trial_i]-Lval[trial_i]));

    }

}"

[Christine C]

I think at least one problem is that binomial and binomial_logit functions take integers as inputs, and I am giving them reals (probabilities)...

[Bob Carpenter]

That is definitely a problem, but that would mean the
program wouldn't compile.

Your problem is that you're evaluating normal(y|...)
with an undefined y at line 50 of your code. You need
to define variables before using them.

- Bob

[Christine C]

Argh, I am sorry for being so dense.  

You are referring to this line?

rflashes[trial_i] ~ normal(rightprob[trial_i], rightprob[trial_i]*prob_m + prob_offset)

But I initialized real rflashes[nTrials] at the beginning of the model block.  How else am I supposed to define it?

I apologize if I this is a stupid question.  Thank you for your quick and helpful responses.

Christine

[Jonah Gabry]

At first glance it looks like you've declared it but not given it a value. Then when you have the "~" statement it tries to compute the log density at the value on the left of the "~", which isn't going to work if it doesn't yet have a value. You'll either need to declare that variable as an unknown parameter in the parameters block, or keep it in the model block but give it a value. If it's something you're trying to estimate then the former and if it's observed (or a function of observed variables) then the latter.

Hope that helps,

Jonah

[Michael Betancourt]

Christine,

Would it be possible to reduce this one more step to just have one sequence

of probabilities?  It would be easier to then see the sequence that you’re

trying to model and facilitate discussion about the modeling.  In addition to

the implementation problem that you’re having I think there’s going to be

plenty of discussion about how to convolve probabilities forward like this.

[Christine C]

Hi Jonah,

Thank you for your reply!  That was super helpful!

Michael,

What do you mean by having one sequence of probabilities?  I am happy to reduce things one more step, if you don't mind clarifying what would be helpful?  Would you like me to clearly describe how probabilities are inferred on a single trial?

Thanks so much for all the feedback!

Christine

On Monday, June 27, 2016 at 1:07:29 PM UTC-4, Christine C wrote:

I am attempting to build a model describing economic choices.  I would like to model noise in subjects' estimates of probability and value.  It is known that people exhibit scalar variability, i.e. the standard deviation in their estimate of magnitude grows linearly with the magnitude.  So, I would like to model the subjects' perception of reward amount/probability, accounting for the linear noise scaling.  Something like this:

perceived_rprob[trial_i] <- normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset);

However, I get the following error: no matches for normal(real, real).  Function normal not found.  So I must be miscasting the variables? I have also tried the line below instead, but from what I can tell from this group and the documentation, the ~ increments the log posterior, and does not mean "drawn from."  This throws an error that the dimension of the right side of the statement is 0.

perceived_rprob[trial_i] ~ normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset); 

I am using rstan; my code is below.  I apologize if I am doing something super dumb, any help is very much appreciated.  Thank you so much!!!

model_code <- "data {

    int nTrials;

    int choices[nTrials];  //values of 1 (chose left) or 2 (chose right)    

    real<lower=0, upper=1> leftprob[nTrials]; 

    real<lower=0, upper=1> leftvol[nTrials];  //rewards are normalized by max reward

    real<lower=0, upper=1> rightprob[nTrials];

    real<lower=0, upper=1> rightvol[nTrials];  //rewards are normalized by max reward

}

transformed data {

// this executes once vs on each sample in the model block, so better to transform here and pay for the conditional once

int<lower=0,upper=1> choices_01[nTrials];

for (i in 1:nTrials){

choices_01[i] <- (choices[i]==2);

}

}

parameters {

    real rho; //exponential of utility function

    real alpha;  //parameter for probability weighting function

    real<lower=0> lambda; //softmax parameter

    real<lower=0> prob_m; //slope of scalar variability relationship

    real<lower=0> prob_offset; //offset of scalar variability relationship

    real<lower=0> vol_m; 

    real<lower=0> vol_offset;

}

model {

    real Rval[nTrials];

    real Lval[nTrials];

    real perceived_rprob[nTrials];

    real perceived_lprob[nTrials];

    real perceived_rvol[nTrials];

    real perceived_lvol[nTrials];    

    rho ~ cauchy(0, 2.5);

    for (trial_i in 1:nTrials) {    

    perceived_rprob[trial_i] <- normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset);

    perceived_lprob[trial_i] <- normal(leftprob[trial_i], leftprob[trial_i]*prob_m+prob_offset);

    perceived_rvol[trial_i] <- normal(rightvol[trial_i], rightvol[trial_i]*vol_m+vol_offset);

    perceived_lvol[trial_i] <- normal(leftvol[trial_i], leftvol[trial_i]*vol_m+vol_offset);  

    Rval <- exp(-(-log(perceived_rprob[trial_i]))^alpha)*(perceived_rvol[trial_i]^rho);

    Lval <- exp(-(-log(perceived_lprob[trial_i]))^alpha)*(perceived_lvol[trial_i]^rho);

    choices_01[trial_i] ~ bernoulli_logit(lambda*(Rval-Lval));

    }

}"

[Michael Betancourt]

I was referring to dropping left and right and just having one set of latent

probabilities, although, upon a more careful look, that would be tricky with

the final likelihood depending on both.

While I don’t full grok your model, there are definitely some issues.  In

particular, you definitely don’t want to do this:

perceived_rprob[trial_i] <- normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset);

which assigns perceived_rprob[trial_i] the value of the normal density and

doesn’t make any statistical sense.

If I understand this correctly, then the model that you want to build is that each user 

in each trial makes a choice based on the difference between Lval and Rval, 

which themselves are based on the user’s perceived left and right probabilities.  

Those perceptions are assumed to be based on the true left and right probabilities, 

which are assumed to be known exactly by the model.

In this case the first big change that you need is to make the perceived

probabilities (and volumes) latent variables, in other words they should all

be parameters sampled by the model.  In the end you marginalize out over

all of these latent perceptions to get at the parameters you care about

while incorporating the uncertainty in not knowing the users’ perceptions

exactly.

Another big modeling issue (the one I had previously referenced) is that the

normal really isn’t compatible with probabilities.  So something like

perceived_rprob[trial_i] ~ normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset)

will lead to an awkward posterior.  What you want for perceived_rprob[trial_i] is a

distribution over probabilities, such as a beta distribution.  We don’t have a mean/variance

parameterization of the beta so you’d have to do this yourself, but it would conceptually

be similar to the above only more appropriate for probabilities.  

Alternatively you could unconstrained everything and then use the normal without

issue, for example

logit_perceived_rprob[trial_i] <- normal(logit_rightprob[trial_i], exp(logit_rightprob[trial_i] * prob_m + prob_offset))

but then you run into problems where some of the rightprob[trial_i] are exactly 0 and 1.

I know you said that the perception variances are supposed to vary in this way and

I am in no way claiming that I am experienced in your field, but these kinds of

conceptual conflicts are _huge_ red flags that almost always indicate that there’s

tension between claims and the underlying statistics.

Finally, the calculation of Lval and Rval introduces some numerical awkwardness that

is going to induce posteriors that are hard to fit.  There are hints of the Folk Theorem

again, now here suggesting another form, but I’ll just leave a warning here and we

can consider it later.

[Christine C]

Although I had to google the term grok, everything else in your response was very clear :) Thanks!  I think you understand the essence of the model.  

In this case the first big change that you need is to make the perceived probabilities (and volumes) latent variables, in other words they should all be parameters sampled by the model.  In the end you marginalize out over all of these latent perceptions to get at the parameters you care about while incorporating the uncertainty in not knowing the users’ perceptions exactly.

Are you suggesting that the perceived left/right probabilities (and volumes) on each trial should be model parameters?  That would mean I would have to have 4x as many model parameters as I have data: for each trial, there would be 4 model parameters (left/right probabilities/volumes).  Is that what you are suggesting?  That seems insane to me. I don't grok it.

I know you said that the perception variances are supposed to vary in this way and I am in no way claiming that I am experienced in your field, but these kinds of conceptual conflicts are _huge_ red flags that almost always indicate that there’s tension between claims and the underlying statistics.

Just to clarify, perception variances have been reported to vary in this way for potentially unbounded stimuli (e.g., estimation of stimulus intensity, time, magnitudes, etc.)  These variances are not thought to apply to perception of probability.  What is unique to my experiment is that I am using stimuli that cause this particular perception variance to convey probability to the subject.  Anyway, not sure if that is helpful/useful at all, but it's a weird artifact of my experiment (I have my reasons, but we'll leave it there).

I think the first thing I will try is to unconstrain everything, as you suggest.  On trials where p=1, I suppose I can set it to p=.999?

Thanks so much!

Christine

On Monday, June 27, 2016 at 1:07:29 PM UTC-4, Christine C wrote:

I am attempting to build a model describing economic choices.  I would like to model noise in subjects' estimates of probability and value.  It is known that people exhibit scalar variability, i.e. the standard deviation in their estimate of magnitude grows linearly with the magnitude.  So, I would like to model the subjects' perception of reward amount/probability, accounting for the linear noise scaling.  Something like this:

perceived_rprob[trial_i] <- normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset);

    perceived_rprob[trial_i] <- normal(rightprob[trial_i], rightprob[trial_i]*prob_m+prob_offset);

    perceived_lprob[trial_i] <- normal(leftprob[trial_i], leftprob[trial_i]*prob_m+prob_offset);

[Michael Betancourt]

Are you suggesting that the perceived left/right probabilities (and volumes) on each trial should be model parameters?  That would mean I would have to have 4x as many model parameters as I have data: for each trial, there would be 4 model parameters (left/right probabilities/volumes).  Is that what you are suggesting?  That seems insane to me. I don't grok it.

Yes!  The problem is that you can’t quantify the uncertainty about parameters that you care about

without first quantifying the uncertainty in those unobserved perceptions.  In other words, in order

for the model to learn the interesting parameters it must first learn the perceptions.  

A cool side effect of this is that you can start analyzing those learned perceptions to potentially

say some interesting things about the experiential population.

I know you said that the perception variances are supposed to vary in this way and I am in no way claiming that I am experienced in your field, but these kinds of conceptual conflicts are _huge_ red flags that almost always indicate that there’s tension between claims and the underlying statistics.

Just to clarify, perception variances have been reported to vary in this way for potentially unbounded stimuli (e.g., estimation of stimulus intensity, time, magnitudes, etc.)  These variances are not thought to apply to perception of probability.  What is unique to my experiment is that I am using stimuli that cause this particular perception variance to convey probability to the subject.  Anyway, not sure if that is helpful/useful at all, but it's a weird artifact of my experiment (I have my reasons, but we'll leave it there).

You just have to get those reasons to align with the statistics!  ;-)

In a very real sense the logit probabilities are equivalent to unbounded stimuli

and hence it would be appropriate to model the logit perception probabilities

that way.  This is one of the big motivations of general linear models.

I think the first thing I will try is to unconstrain everything, as you suggest.  On trials where p=1, I suppose I can set it to p=.999?

This is the open question to which I do not have a great answer.  Setting

those extremes to 0.999 might work okay but I couldn’t guarantee it.  You

might have to resort to a probabilistic model, for example with beta distributions.

Anyone else have any ideas?

[Jonah Gabry]

I think the first thing I will try is to unconstrain everything, as you suggest.  On trials where p=1, I suppose I can set it to p=.999?

This is the open question to which I do not have a great answer.  Setting

those extremes to 0.999 might work okay but I couldn’t guarantee it.  You

might have to resort to a probabilistic model, for example with beta distributions.

Anyone else have any ideas?

Hi Christine, I don't have a great suggestion for this either, but I think you could start a new thread just for this particular question and probably get some more eyes on it that way (at the moment it's kind of buried in this long thread).

Jonah