Categorical Logit, with exponent in utility function

Modeling
1

Hello Stan users,

I’m trying to write Stan code that is both a) a categorical logit and b) is a ‘structural model’ (economics terminology). I can find useful tutorials/code on either the former (e.g. https://github.com/ksvanhorn/ART-Forum-2017-Stan-Tutorial) or the latter (e.g. https://github.com/CCS-Lab/hBayesDM/blob/master/exec/ra_prospect.stan) but not the combination, and have got stuck in adapting the code for my purpose.

My data comes from 300 people choosing one of 11 options, in two rounds. It is a risky choice, with a simplified utility function below. (I’m aware ‘normal’ categorical logits need K-1 for identification, I don’t think that is true here.)

Previously, I calculated the utility of each option and the code took ~3 hours to run. It appears vectorising is a sensible option to improve speed, but I’ve come up against the problem that you can’t use exponents on vectors. [error: No matches for: pow(vector, real)]

The answer appears to be looping over choices in the model block, but I’m getting confused between the different indices, as R respondents make a choice between C alternatives in S rounds. Furthermore, this seems to defeat the point of vectorising. If I should loop over choices to use the exponential, how do I do this? If I shouldn’t, or a better set up would be optimal, what combination of vector/array/indexing options seems best? Parameters will either be the same for everyone, or have some distribution over respondents (but constant across scenarios and alternatives).

Current code below

data {
	int<lower=2> C; // choices/alternatives per scenario=11
	int<lower=2> R; // respondents=300ish
	int<lower=2> S; // scenarios=2
	int<lower=1> N;	// observations=R*S
	int<lower=1,upper=C> Y[N]; // number of chosen alternative
	vector[C] gain[N]; 
	vector[C] loss[N]; 
	vector[C] EV[N]; 
}

parameters {
	real<lower=0> lambda[N];
	real mean_lambda;  
	real<lower=0> sigma_lambda;
	real<lower=1> r;
}

model {
  mean_lambda ~ cauchy(5, 1);  
  sigma_lambda ~ cauchy(0, 1);  
  lambda ~ normal(mean_lambda, sigma_lambda);
  r ~ cauchy(1,1);
  
	for (i in 1:N){
		Y[i] ~ categorical_logit(EV[i] + gain[i] - lambda[i]*pow(loss[i],r));
}}
1 Like
2

This should work (make a temporary to hold the output of pow):

for (i in 1:N){
  vector[C] ptmp;
  for(c in 1:C) {
    ptmp[c] = pow(loss[i, c], r);
  }

  Y[i] ~ categorical_logit(EV[i] + gain[i] - lambda[i] * ptmp);
}

A little messy for now, but pow will get vectorized eventually :D

3

Oh, messy is definitely fine for now!

Many thanks: the model is now running smoothly.

4

We also added the syntax loss[i, c]^r if you prefer.

And yes, we shold vectorize this sooner rather than later. Rayleigh was working on a general case, but we should probably just vectorize pow itself sooner.