Advice on optimizing Stan code: Nested mixed logit

[al hanz seiji lidasan]

Hi, we are currently testing the application of Stan for discrete choice models for mode choice in transportation using simple simulated data.
Specifically, the model we are testing is nested mixed logit with 3 choices namely, "car", "rail", "bus" with rail and bus correlated in one nest due to their public transit nature and car in another nest.
The correlation is modeled as error terms (eps) that comes from a single distribution N(0,sigma).
We were able to recover the results of the toy model by comparing the Stan model to a maximum-likelihood estimate.
However, we are concerned regarding the very low effective number of samples.
Is there a way to improve the estimate of the effective number of samples?
I have attached our Stan code and the simulated data so it can be replicated.
The Stan code, result of the Stan code and the maximum likelihood estimate are shown below.


data <- read.csv(file="35.csv",header=T)
ns <- nrow(data); n.alt <- 3; n.var <- 1

mnl_code <-  "
data {
    int<lower=2> K;  //number of alternatives
    int<lower=0> N;  //number of samples
    int<lower=1> D;  //number of variables
    int<lower=0,upper=1> y[N,K];  //outcome for each sample
    matrix [K,D]  x[N];  //matrix of variables
}

parameters {
    vector[D] beta;
    real<lower=0> sigma;  //hypersigma
    vector[N] eps; //error terms
}

//transformed parameters {
//}

model {
    beta ~normal(0,10);
    eps ~ normal(0,sigma);
    sigma ~ normal(0,1);
    for (n in 1:N) {
           vector[K] v; vector[K] lp;
               v = x[n]*beta;
                  lp[1] = log( exp(v[1]        )/( exp(v[1])+exp(v[2]+eps[n])+exp(v[3]+eps[n])) );
                  lp[2] = log( exp(v[2]+eps[n] )/( exp(v[1])+exp(v[2]+eps[n])+exp(v[3]+eps[n])) );
                  lp[3] = log( exp(v[3]+eps[n] )/( exp(v[1])+exp(v[2]+eps[n])+exp(v[3]+eps[n])) );
                  target += lp[1]*y[n,1]+lp[2]*y[n,2]+lp[3]*y[n,3];
    }
}
"


Stan

1 chains, each with iter=4000; warmup=2000; thin=1;

post-warmup draws per chain=2000, total post-warmup draws=2000.

            mean se_mean     sd     2.5%      25%      50%      75%    97.5% n_eff Rhat

beta[1]     5.36    0.04   0.36     4.71     5.11     5.34     5.59     6.07    99 1.00

sigma       1.85    0.04   0.26     1.38     1.67     1.82     2.01     2.41    43 1.01

lp__    -1734.57   18.15 115.77 -1962.33 -1814.64 -1726.38 -1656.17 -1504.12    41 1.01

Maximum Likelihood

Log-Likelihood: -783.4759

2  free parameters

Estimates:

     Estimate Std. error t value  Pr(> t)   

[1,]   5.3705     0.4036  13.308  < 2e-16 ***

[2,]   1.8679     0.3187   5.861 4.61e-09 ***

Would really appreciate some advice and inputs.

[Sebastian Weber]

Hi!

I really don't get your model, but looking at the code you may want to revise the inner loop like this:

    for (n in 1:N) {
           vector[K] v; vector[K] lp;
               v = x[n]*beta;

               real lscale;
  // pre-calcluate common factor on the log-scale and switch to non-centered parametrization
               lscale <- log_sum_exp( v[1], log_sum_exp(v[2] + sigma * xi[n], v[3] + sigma * xi[n])  );
                  lp[1] =v[1] - lscale;

                  lp[2] = v[2]+sigma * xi[n] - lscale;
                  lp[3] = v[3] + sigma * xi[n] - lscale;

                  target += lp[1]*y[n,1]+lp[2]*y[n,2]+lp[3]*y[n,3];
    }

then define instead of eps the variable xi and assign it a normal(0, 1) prior.

Checkout the manual on the centered vs non-centered parametrization and try to avoid switching natural / log space;prefer to stay on the log scale whenever you can and use specialized functions like log_sum_exp as needed.

Sebastian

[Bob Carpenter]

First off, 2.10 introduced a bug which we have already
patched. So hopefully you're using 2.9 or earlier or already
have 2.11.

You can't diagnose convergence with one chain---you want to use at
least two and preferably four.

In terms of n_eff, the non-centered parameterization is the
best thing to try.

But you also have an extra degree of freedom, so the model's
not strictly identifiable in the likelihood, so the priors are
driving convergence. There's a discussion in the manual regression
chapter under heading "identifiability" (p. 60 for v2.10).

You also want to make the arithmetic more robust. Basically,
avoid use of exp() which can overflow and underflow easily. See below.

I take it the eps is the effect of public transportation?

You can use matrix multiply for x * beta:

vector[K] v;
v = x * beta;

You should gather and share all subterms rather than reapplying exp(), etc.

We have a built-in function to do all of this:

vector[3] zz;

zz[1] = v[n, 1];
zz[2] = v[n, 2] + eps[n];
zz[3] = v[n, 3] + eps[n];
target += y[n] * log(softmax(zz));

That should also be a bit more stable arithmetically.

This is where the extra degree of freedom is. You can add
any constant to v and get the same answer. Then if
you pin v[n,1] = 0, you get the same problem in that
you can't identify v and eps because you can add to one
and subtract from the other.

- Bob

> }
> }
> "
>
>
> Stan
> 1 chains, each with iter=4000; warmup=2000; thin=1;
> post-warmup draws per chain=2000, total post-warmup draws=2000.
> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
> beta[1] 5.36 0.04 0.36 4.71 5.11 5.34 5.59 6.07 99 1.00
> sigma 1.85 0.04 0.26 1.38 1.67 1.82 2.01 2.41 43 1.01
> lp__ -1734.57 18.15 115.77 -1962.33 -1814.64 -1726.38 -1656.17 -1504.12 41 1.01
>
>
> Maximum Likelihood
> Log-Likelihood: -783.4759
> 2 free parameters
> Estimates:
> Estimate Std. error t value Pr(> t)
> [1,] 5.3705 0.4036 13.308 < 2e-16 ***
> [2,] 1.8679 0.3187 5.861 4.61e-09 ***
>
>
> Would really appreciate some advice and inputs.
>

> --
> 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.
> <35.csv><MXL_stan_rnd3_2r.R>