2

I have a draw from a Dirichlet distribution that I would like to use as input into a MultinomialDistribution. For instance:

rand = RandomVariate[
  DirichletDistribution[{1, 1, 1, 1, 1, 239, 302, 1}]]

Then I want to take this list and put it into a MultinomialDistribution and draw a sample from that. For instance:

RandomVariate[MultinomialDistribution[1, rand]]

It throws an error.

"The value \
{2.02597*10^-6,0.000659251,0.00210503,0.000496145,0.00063576,0.426032,\
0.564779} at position 2 in \
MultinomialDistribution[1,{2.02597*10^-6,0.000659251,0.00210503,0.\
000496145,0.00063576,0.426032,0.564779}] is expected to be a list of \
numbers greater than or equal to 0 and summing to 1. "

Probably because of:

Total[rand]
0.994709

My understanding is that the sum of the values drawn from a Dirichlet should equal one. However, in this case, they don't and I'm not sure why. Is this a mathematica issue or an issue with my understanding of the Dirichlet. I suspect the prior because I can do this in python (with numpy) okay and they always add to one. 0.99 is pretty close to 1. Is there perhaps a workaround for this?

I can reproduce this on an even smaller set:

bins = {1, 2, 3, 1, 10}
Total[RandomVariate[DirichletDistribution[bins]]]
Length[RandomVariate[DirichletDistribution[bins]]]

Returns:

0.443586
4
Tom Hayden
  • 411
  • 3
  • 9

2 Answers2

5

If a random vector $x= \{x_1, x_2,...,x_k\}$ has Dirichlet distribution of order k>2 with parameters $\alpha_1,\alpha_2,...,\alpha_k$ its kth component is determined by the condition $x_k= 1- x_1 - x_2 - \cdots- x_{k-1}$ (see Wikipedia > Dirichlet Distribution). So you need to append 1-Total[rand] to rand to get the required parameter for the MultinomialDistribution:

rand = RandomVariate[DirichletDistribution[{1, 1, 1, 1, 1, 239, 302, 1}]];
RandomVariate[MultinomialDistribution[1, Join[rand, {1 - Total@rand}]], 5]
(* {{0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 0}, 
    {0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0}} *)

Alternatively, you can use ParameterMixtureDistribution as follows:

p = {p1, p2, p3, p4, p5, p6, p7, 1 - p1 - p2 - p3 - p4 - p5 - p6 - p7};
dist = ParameterMixtureDistribution[MultinomialDistribution[1, p], 
        Distributed[Most@p, DirichletDistribution[{1, 1, 1, 1, 1, 239, 302, 1}]]];

RandomVariate[dist, 5]
(* {{0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 0}, 
    {0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0}} *)
kglr
  • 394,356
  • 18
  • 477
  • 896
  • Thank you! Kudos to mathematica for being so close to the theory even though it's weird from a computational sense :) – Tom Hayden Jan 07 '15 at 14:41
1

Use

Normalize[rand, Total]

I believe I have asked this before in here.

In my case, MMA is very misleading as my probabilities do sum to 1.

Community
  • 1
Chen Stats Yu
  • 4,986
  • 2
  • 24
  • 50
  • Thank you! I apologize for double asking this question! – Tom Hayden Jan 06 '15 at 23:09
  • @TomHayden No problem. Mathematica is very confusion most of the time. – Chen Stats Yu Jan 06 '15 at 23:10
  • 1
    I believe these two questions to be distinct. Chen's question was about why the error message occurred (the numbers added almost infinitesimally close to 1) and didn't use RandomVariates, but in this case the RandomVariates' total clearly differs from 1, and one could wonder whether this isn't a bug. Normalizing restores the correct total to 1, but is that legal? That is to say, do you get a real Dirichlet distribution that way? – Sjoerd C. de Vries Jan 06 '15 at 23:30
  • Along similar lines, is it possible that RandomVariate on DirichletDistribution is not returning all the values. For instance, when I get 5 inputs, I get 4 inputs back: bins = {1, 2, 3, 1, 10} indexes = {1, 2, 3, 4, 5} RandomVariate[DirichletDistribution[bins]] {0.0206046, 0.0435381, 0.112551, 0.145855} – Tom Hayden Jan 06 '15 at 23:52