# Discrete choice with correlated random parameters

**URL:** <https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185>\
**Category:** v5\
**Tags:** modeling\
**Created:** [July 10, 2025, 4:12pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185 "2025-07-10T16:12:39Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![dekkert](https://avatars.discourse-cdn.com/v4/letter/d/aeb1de/32.png) [@dekkert](https://discourse.pymc.io/u/dekkert)\
**Post date:** [July 10, 2025, 4:12pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/1 "2025-07-10T16:12:39Z")

</div>

Dear all,

After years of working with and teaching Gibbs Samplers in the context of discrete choice (i.e. multinomial logit/probit), I’ve decided to explore the world of Pymc, and for what I’ve seen I’m very impressed by its versatility.

I have a mixed logit (or hierarchical bayes) running with two random parameters with the following random parameter structure:  
#RP 1:  
mu\_asc = pm.Normal(“mu\_asc”, 0, 10)  
sigma\_asc= pm.InverseGamma(“sigma\_asc”, alpha=3, beta=0.5)  
asc\_alt1 = pm.Normal(“asc\_alt1”, mu=mu\_asc, sigma=np.sqrt(sigma\_asc),dims=(“individuals”))

#RP2:  
mu\_tc = pm.Normal(“mu\_tc”, 0, 10)  
sigma\_tc= pm.InverseGamma(“sigma\_tc”, alpha=3, beta=0.5)  
beta\_tc = pm.Normal(“beta\_tc”, mu=mu\_tc, sigma=np.sqrt(sigma\_tc), dims=(“individuals”))

My intended extension is to estimate the correlation matrix between these. Normally, I’d use the inverseWishart but that is not supported using nuts samplers. The LKJ setup from the examples (corr across alternatives) looks suitable but I’m lost in two aspects regarding implementation:

1. dimensions: I’ve got dims individuals and dims random\_vars, so for each individual 2 random parameters. What is the relevant dims to include in the pm.MvNormal()? and subsequently when referring to it in pm.Model?  
rps = pm.MvNormal(“rp”, mu=mu\_rp, chol=chol\_rp, dims=“random\_vars”)  
u1 = rps[person\_indx][0] + beta\_tt \* database[“tt1”] - pt.exp(rps[person\_indx][1]) \* database[“tc1”]

2. I find the labelling of the pm.LKJCholeskyCov() slightly confusing. How do you properly setup the prior for the 2x2 cholesky matrix chol\_rp referred to above in pm.MvNormal?

chol\_rp, corr, stds = pm.LKJCholeskyCov(  
“chol\_rp”, n=2, eta=2.0, sd\_dist=pm.HalfNormal.dist(10)  
)

Many thanks for your feedback,  
Thijs

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 10, 2025, 5:18pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/2 "2025-07-10T17:18:16Z")

</div>

The dimension of length 2 (random\_vars)? should be on the right of your MvNormal, so something like `dims=("individuals", "random_vars")` to get (n, 2) draws where the correlation structure acts on the last dimension. More details in [Distribution Dimensionality — PyMC 5.23.0 documentation](https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/dimensionality.html)

The LKJ prior looks correct, what are you unsure about?

---

<div class="post-metadata">

**Author:** ![dekkert](https://avatars.discourse-cdn.com/v4/letter/d/aeb1de/32.png) [@dekkert](https://discourse.pymc.io/u/dekkert)\
**Post date:** [July 10, 2025, 9:43pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/3 "2025-07-10T21:43:01Z")

</div>

Thanks for the quick reply.

Re LKJ, why is there the need to start with “_chol, corr, stds =_” and not just the one term intended for use?

The model now progresses to the compilation stage and replicates my GS results :

```python
N = database.shape[0]
observed = pd.Categorical(database["choice"]).codes
person_indx, uniques = pd.factorize(database["ID"])
coords = {      
    "alts_probs": ["1", "2"],
    "random_vars": ["asc_alt1", "beta_tc"],
    "individuals": uniques,
    "obs": range(N),
}

with pm.Model(coords=coords) as model_1:
    
    # Priors for Fixed parameters
    beta_tt = pm.Normal("beta_tt", 0, 10)    
    beta_hw = pm.Normal("beta_hw", 0, 10)
    beta_ch = pm.Normal("beta_ch", 0, 10)

    # Hierarchical priors on the random parameters
    mu_rp = pm.Normal("mu_rp", 0, 10, dims="random_vars")
    chol_rp, corr, stds = pm.LKJCholeskyCov(
        "chol_rp", n=2, eta=2.0, sd_dist=pm.HalfNormal.dist(10)
    )
    rps = pm.MvNormal("rp", mu=mu_rp, chol=chol_rp, dims=("individuals", "random_vars"))
 
    ## Construct Utility matrix and Pivot
    u1 = rps[person_indx,0] + beta_tt * database["tt1"] - pt.exp(rps[person_indx,1]) * database["tc1"] + beta_hw * database["hw1"] + beta_ch * database["ch1"]
    u2 = beta_tt * database["tt2"] - pt.exp(rps[person_indx,1]) * database["tc2"] + beta_hw * database["hw2"] + beta_ch * database["ch2"]
       
    s = pm.math.stack([u1, u2]).T

    ## Apply Softmax Transform
    p_ = pm.Deterministic("p", pm.math.softmax(s, axis=1), dims=("obs", "alts_probs"))

    ## Likelihood
    choice_obs = pm.Categorical("y_cat", p=p_, observed=observed, dims="obs")

    idata_m1 = pm.sample(tune=4000, draws=4000)    

pm.model_to_graphviz(model_1)
```

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 10, 2025, 10:08pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/4 "2025-07-10T22:08:40Z")

</div>

> [@dekkert](#):
>
> Re LKJ, why is there the need to start with “_chol, corr, stds =_” and not just the one term intended for use?

People like to see them I guess. There’s some keyword arguments to disable them but then you need to call expand\_triangular manually to get the cholesky.

---

<div class="post-metadata">

**Author:** ![bob-carpenter](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bob-carpenter/32/8124_2.png) [@bob-carpenter](https://discourse.pymc.io/u/bob-carpenter)\
**Post date:** [July 15, 2025, 2:56pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/5 "2025-07-15T14:56:18Z")

</div>

Now that you’ve escaped Gibbs, you can have more flexible priors—you’re not required to use conjugate priors on variance like the inverse gamma. I find it easier to put priors directly on the scales (standard deviation for normals) because they have the same units as the mean and the observations.

You can also put your own covariance matrix together with an LKJ prior on the correlations (as parameter increases, the prior concentrates around unit correlation matrices) and whatever prior you want on the scales by multiplying after you’re done—I don’t know how flexible the `sd_dist` argument is.

In 2D, you have a lot of options. You can model two scales \sigma\_1, \sigma\_2 \> 0 and a correlation \rho \in (-1, 1) and create the covariance matrix [[\sigma\_1^2 \quad \sigma\_1 \cdot \sigma\_2 \cdot \rho], [\sigma\_1 \cdot \sigma\_2 \cdot \rho \quad\sigma\_2^2]]. This lets you put whatever priors you want on \sigma and \rho as long as they respect the constraints. The LKJ corresponds roughly to taking \phi \sim \textrm{beta}(a, a) for a \>\> 1 and then setting \rho = 2 \cdot \phi - 1 \in (-1, 1). But in the general case you can use a \beta(a, b) prior if you want asymmetry, you can impose additional constraints (like the correlation must be positive), etc.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 15, 2025, 3:31pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/6 "2025-07-15T15:31:49Z")

</div>

> [@bob-carpenter](#):
>
> Now that you’ve escaped Gibbs, you can have more flexible priors

That’s gonna go on my quote wall

---

<div class="post-metadata">

**Author:** ![dekkert](https://avatars.discourse-cdn.com/v4/letter/d/aeb1de/32.png) [@dekkert](https://discourse.pymc.io/u/dekkert)\
**Post date:** [July 15, 2025, 3:54pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/7 "2025-07-15T15:54:46Z")

</div>

Ha! Fully appreciate that Bob, and I will dig deeper over the summer. Here, I stayed in the conjugate world for replicability of my ‘_trusted_’ Gibbs results.

On a related note, I’m doing some portfolio choice models where Gibbs makes the application scalable to larger choice sets by augmenting the utilities of the components of the portfolios instead of working with the 2^J potential combinations! There’s still life in the old beast 🙂

---

<div class="post-metadata">

**Author:** ![bob-carpenter](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bob-carpenter/32/8124_2.png) [@bob-carpenter](https://discourse.pymc.io/u/bob-carpenter)\
**Post date:** [July 15, 2025, 4:14pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/8 "2025-07-15T16:14:07Z")

</div>

What led you to trust the Gibbs results?

There are some highly specialized problems in which Gibbs can be custom tuned with blocking to be faster than HMC. But the fundamental problem is that it scales very poorly in dimension and with correlated parameters.

> [@dekkert](#):
>
> Gibbs makes the application scalable to larger choice sets by augmenting the utilities of the components of the portfolios instead of working with the 2^J potential combinations!

I’m not exactly sure what you mean, but this sounds like a standard random effects approach and is something you can do more efficiently with HMC than with Gibbs. The problem with Gibbs is that it scales very poorly in dimension and very poorly in the face of correlated parameters that cannot be blocked and sampled jointly. For example, one cannot sample from a 100 dimensional highly correlated normal with Gibbs—it’ll take forever to make progress (see the plot in the original Hoffman and Gelman NUTS paper, for example, or try it yourself).

When we started building Stan, we first pulled the lid off JAGS and sped that up considerably (it’s very inefficiently implemented on top of R’s C++ library) until we realized that no matter how fast we made Gibbs, it just wouldn’t sample even moderately hard problems in higher dimensions if the correlations were too high.

---

<div class="post-metadata">

**Author:** ![dekkert](https://avatars.discourse-cdn.com/v4/letter/d/aeb1de/32.png) [@dekkert](https://discourse.pymc.io/u/dekkert)\
**Post date:** [July 15, 2025, 9:50pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/9 "2025-07-15T21:50:55Z")

</div>

> [@bob-carpenter](#):
>
> What led you to trust the Gibbs results?

Rephrase…the gibbs sampler aligns with results from MLE and EM; combined with enough experience to see the sampler go off in cases where there is strong correlation across the blocks.

In the portfolio models we ask people to allocate gov budget B to a range of projects varying in cost. This leaves three sets of projects, M - those in the portfolio, K - those not chosen but still affordable (true zero), and Z - those not chosen but not affordable given the remaining budget. Each project has a random utility, and using KKT conditions we get a ranking of the utility per unit of budget determining the likelihood - see [Dekker et al. (2024)](https://doi.org/10.1016/j.jocm.2024.100507) case 1 for details.

V\_{m}+\epsilon\_{m}-ln(p\_m)\>V\_0+\epsilon\_{0}\>V\_{k}+\epsilon\_{k}-ln(p\_{k}), \forall m,k and  
V\_{m}+\epsilon\_{m}-ln(p\_m)\>V\_{z}+\epsilon\_{z}-ln(p\_{z}), \forall m,z

In the above good zero is the numeraire of the remaining budget with price of 1. We get a closed form solution using extreme value distributions, but this involves set theory and enumerating the full factorial of the sets Z and K (eq 18 and appendix A). This is where the issues start with increasing number of alternatives that go into the decision problem. This works but slow with MLE and results verified with MH random walk.

Gibbs is a natural candidate (especially when switching to normals) as by augmenting V\_{j}+\epsilon\_{j}-ln(p\_{j}) the model reduces to a SUR structure to which we ultimately could introduce correlation. Indeed, normalisation is needed an probably best to just assume V\_{0}+\epsilon\_{0}=0. Even just augmenting using extreme value + MH random walk speeds up and enables replication of above results.

I’d be happy to hear how you would approach this estimation challenge using HMC without having to work with the inconvenient likelihood function presented in the referred paper.

---

<div class="post-metadata">

**Author:** ![bob-carpenter](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bob-carpenter/32/8124_2.png) [@bob-carpenter](https://discourse.pymc.io/u/bob-carpenter)\
**Post date:** [July 16, 2025, 2:31pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/10 "2025-07-16T14:31:03Z")

</div>

> [@dekkert](#):
>
> the gibbs sampler aligns with results from MLE and EM

Aligns in the sense that posterior mean estimates from the sampler correspond to the MLE? Or that the Bayesian posterior intervals correspond to the frequentist confidence intervals? Either way, if that’s your measure of success, why not just stick with frequentist methods? That’s not a rhetorical question—if you have a model you can fit with EM, it’s going to be way faster and more robust in most cases than sampling.

I’m afraid that paper’s way too detailed for me to have time to figure it out—it’s really hard jumping into an applied field to try to understand the language the use to talk about models and the underlying conventions they follow.

If it’s just an augmented variable approach, you can do that with HMC, too. As a simple example, suppose we have something like y\_n \sim \textrm{negBinomial}(\alpha, \beta). We can replace that with \lambda\_n \sim \textrm{gamma}(\alpha, \beta) and y\_n \sim \textrm{Poisson}(\lambda\_n). In many cases, the geometry of the posterior (in terms of conditioning of Hessian and stability of its eigenstructure around the posterior) is much better, despite the huge increase in dimensionality. For example, I’m able to do this for 20 replicates in a genomics experiment with 20K splice variants, or about 400K latent \lambda and it fits OK (not super fast, of course, but we’re talking an hour, not weeks).

Is “MH random walk” a Metropolis-Hastings Markov chain with accept/reject? If MH is enough for your problem, why not stick with that? Is it too slow? Will it not scale in dimension (MH or Gibbs rarely scale well in dimension unless there’s a lot of conditional independence in the model)? Is it too hard to generalize?

---

<div class="post-metadata">

**Author:** ![Nathaniel\_Forde](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/nathaniel_forde/32/4041_2.png) [@Nathaniel\_Forde](https://discourse.pymc.io/u/Nathaniel_Forde)\
**Post date:** [July 17, 2025, 12:47pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/11 "2025-07-17T12:47:47Z")

</div>

> [@dekkert](#):
>
> Dekker et al. (2024)

Love to see it! Glad to see people using PyMC for discrete choice

---

<div class="post-metadata">

**Author:** ![dekkert](https://avatars.discourse-cdn.com/v4/letter/d/aeb1de/32.png) [@dekkert](https://discourse.pymc.io/u/dekkert)\
**Post date:** [July 17, 2025, 9:02pm UTC](https://discourse.pymc.io/t/discrete-choice-with-correlated-random-parameters/17185/12 "2025-07-17T21:02:32Z")

</div>

Lots to digest here @bob-carpenter and no expectation at all for you or anyone else to work their way through the paper (more than welcome of course if interested ;-)).

Ultimately, the problem is similar to the multinomial probit model which has received some attention on Stan forums. In the MNP case the ordering of the utilities for different alternatives determines choice (select the one with the max utility). In our case being in or out of the portfolio is determined by relative attractiveness (i.e. add those with high util and drop those with low utils and stay within the budget constraint with the overall portfolio cost). From the Stan discussions I now understand why HMC struggles with this.

Nevertheless, this has been helpful and I have some ideas for simplifying the model structure and then combine it with the hierarchical LKJ correlation structure which should produce (close to) identical results. Will open a new thread when I get stuck implementing this with CustomDist().