8

this is a head-scratcher for me, but a very interesting problem. So I have a stochastic simulation model for a hiring process. Basically different groups get hired into a company with different probabilities. So there is a multinomial process for hiring, where $n$ number of available positions go to members from $m$ number of groups in each year. So in a Bayesian sense the model would look like below--I left out the details for the prior parameters since they are not really important:

$$ n \sim Binomial(\cdot, \cdot) \\ rate1 \sim TruncatedNormal(\cdot, \cdot)\\ rate2 \sim TruncatedNormal(\cdot, \cdot)\\ rate3 \sim TruncatedNormal(\cdot, \cdot)\\ rate4 \sim TruncatedNormal(\cdot, \cdot)\\ \\ \\ hires \sim Multinomial(n, [rate1, rate2, rate3, rate4]) \\ $$

I have hiring data, and I was going to use Bayesian MCMC sampling to estimate the hiring rates for each group from the time-series data.

The challenge is that when I draw samples from the prior for each of the hiring rates, those rates have to sum to 1--otherwise the multinomial draw will not work. So I need to explore the set of hiring rates in such a way that all of the rates sum to 1.

I am not sure how to handle this kind of constraint in an MCMC scheme. Is there a distribution that I would sample the rates from such that they all sum to 1? It seems like a very interesting problem, but I only discovered it while debugging my own code, so I have not had much chance to think about it yet.

Any thoughts would be appreciated.

krishnab
  • 1,008
  • 7
  • 12
  • 1
    Have you considered looking at the [Dirichlet distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution)? It is a probability distribution over discrete probability vectors. This would have the added bonus of [Dirichlet-multinomial conjugacy](https://en.wikipedia.org/wiki/Dirichlet-multinomial_distribution). It might be important to check if that is the multinomial you had in mind, as it can be unclear whether authors are referring to the categorical or multinomial distribution. – microhaus Jun 03 '21 at 23:32
  • I have extremely limited experience with MCMC, however, I have seen a paper before where the authors implement a sum-to-zero constraint trick in their scheme. Even though it's not sum-to-one, I will post it here when I have managed to find it, as it also shows how to implement this in code, and may be of interest. – microhaus Jun 03 '21 at 23:38
  • @microhaus thanks for the suggestion here. I can think about the Dirichlet distribution, as that is a distribution over multinomials. Someone must have implemented something like that before, right. The other thing I can kind of think of is those convolutions in probability distributions, where you have the Probability of P(Z) given X + Z = 1. I think they solved those with characteristic functions or such, but it has been a while since I looked at that stuff :). – krishnab Jun 04 '21 at 00:08
  • @microhaus you know, I did a little checking on the applications of the Dirichlet and it does seem to be used commonly for these types of problems. I have run into the Dirichlet before, but never had to use it myself. So that was a good call. Do you think it makes sense to keep this question open, or should I cancel it? – krishnab Jun 04 '21 at 01:22
  • 2
    What I was referring to in my earlier comment is a blog post [here](http://danielweitzenfeld.github.io/passtheroc/blog/2014/10/28/bayes-premier-league/). The blog uses the strategy from [*Baio, G., Blangiardo, M. (2011). Bayesian hierarchical model for the prediction of football results*](https://discovery.ucl.ac.uk/id/eprint/16040/1/16040.pdf) on a different data set. There are examples of how to code a sum-to-zero constraint (but not sum-to-one constraint for probability vectors) in both. IMO there is no need to close the question, you might receive concrete advice from MCMC experts. – microhaus Jun 04 '21 at 01:31

1 Answers1

8

The problem does not seem to stand with MCMC but with the prior modelling. If the data comes from a Multinomial distribution $$\mathcal D_k(n,p_1,\ldots,p_k)$$ where the probability vector $\mathbf{p}=(p_1,\ldots,p_k)$ belongs to the $k$-dimensional simplex, the prior on $\mathbf{p}$ must put some mass on the simplex (and should logically be allocating all its mass to the simplex. Using truncated Normals is thus inadequate as the prior puts zero mass on the simplex.

As suggested by microhaus, a natural family of distributions in this setting is the family of Dirichlet priors. An alternative (or a generalisation) is to over-parametrise the model and set the constraint later. For instance, take $$p_1=\dfrac{\rho_1}{\sum_{i=1}^k \rho_i},\ldots,p_k=\dfrac{\rho_k}{\sum_{i=1}^k \rho_i}$$ with$$\rho_i\sim\mathcal N^+(a_i,b_i^2)$$ (when using a Gamma prior with a constant scale instead of a truncated Normal, this returns a Dirichlet).

Xi'an
  • 90,397
  • 9
  • 157
  • 575
  • 1
    would the overparametrised version not suffer from convergence issues? – seanv507 Jun 04 '21 at 09:00
  • Since the prior is proper, the posterior is proper as well, hence MCMC has no convergence issue in theory. In the opposite, since the scale of the $\rho_i$'s is unconstrained by the data, the chain should converge better. – Xi'an Jun 04 '21 at 09:03
  • What I am trying to refer to is slow convergence eg https://mc-stan.org/docs/2_20/stan-users-guide/sampling-difficulties-with-problematic-priors.html - example with weak prior. https://mc-stan.org/docs/2_27/stan-users-guide/collinearity-section.html#collinearity.section where they discuss softmax as an example.. – seanv507 Jun 04 '21 at 09:36
  • I see no issue in this case: the dimension does not increase (or just by one), the overparameterised model is less constrained, the prior is proper, and generic MCMC algorithms should thus move more freely in the overparameterised case. In addition, the scale of the $p_i$'s does not need to converge. – Xi'an Jun 04 '21 at 09:42
  • @Xi'an Thank you so much. I really appreciate the clear explanation of the reasoning in the answer. Yes, I can try what you suggest. I was worried that with the rescaling of every $p_1$ that I would not be able to obtain the original unscaled $\rho_i$. So at the end of the day, would I just take the $p_k$ as the parameter estimates for the hiring rates, or would I used the $\rho_k$? – krishnab Jun 04 '21 at 13:57
  • @krishnab: just consider which ones between the $p_k$'s and the $\rho_k$'s can be used in the multinomial distribution. – Xi'an Jun 04 '21 at 13:58
  • 1
    Yep, so I would have to use the normalized $p_k$s, as they are the only ones that fall within the probability simplex, I imagine. In other words, the $\rho_k$ seem to be thrown away, or just a means to get the normalized $p_k$ in the simplex. And that makes sense. – krishnab Jun 04 '21 at 14:58
  • 1
    Haha, also that trick about the Gamma with constant scale is actually a Dirichlet is ingenious. That is great. I just noticed that. – krishnab Jun 04 '21 at 15:16