3

I want to obtain posterior distribution for parameters of a Dirichlet distribution $x = (p_1,p_2,p_3) \sim Dir(p_1,p_2,p_3; a_1,a_2,a_3)$ with uniform $P(a_1,a_2,a_3)$ and observed data $X=\{x_1,x_2,...,x_n\}$. What I have is:

$$Pr(a_1,a_2,a_3 | X) \propto P(X|a_1,a_2,a_3)= \Pi_{i}^nP(x_i|a_1,a_2,a_3)$$ $$=\Bigg[\frac{\Gamma(a_1 + a_2 + a_3)}{\Gamma(a_1)\Gamma(a_2)\Gamma(a_3)}\Bigg]^n (\Pi_i^nx_{i1})^{a_1}(\Pi_i^nx_{i2})^{a_2}(\Pi_i^nx_{i3})^{a_3} (*)$$ How can we sample $(a_1,a_2,a_3)$ from this exotic distribution? Is a more intuitive distribution for $$Pr(a_1,a_2,a_3 | X) \propto \Bigg[\frac{\Gamma(a_1 + a_2 + a_3)}{\Gamma(a_1)\Gamma(a_2)\Gamma(a_3)}\Bigg](\frac{\sum_i x_{i1}}{\sum_i x_{i1}+\sum_i x_{i2}+\sum_i x_{i3}})^{a_1}(\frac{\sum_i x_{i2}}{\sum_i x_{i1}+\sum_i x_{i2}+\sum_i x_{i3}})^{a_2}(\frac{\sum_i x_{i3}}{\sum_i x_{i1}+\sum_i x_{i2}+\sum_i x_{i3}})^{a_3}$$ ? Is there a mathematical way to derive this? Empirically, it seems to give a good posterior for the distribution of $(a_1,a_2,a_3)$, but I cannot see how to derive it. Any help is appreciated.

Note that the data $X$ is generated by the Dirichlet distribution, there is no multinomial distribution here. $X$ is a collection of tupples of ratios (the elements of each data points in $X$ is less than 1).

kjetil b halvorsen
  • 63,378
  • 26
  • 142
  • 467
TuanDT
  • 195
  • 6
  • 1
    Possible duplicate of [Dirichlet posterior](http://stats.stackexchange.com/questions/41882/dirichlet-posterior) – Xi'an Jun 08 '16 at 04:55
  • @Xi'an Thank you for the link, very informative. My question is slightly different though, it's edited for clarity. I also misunderstood the problem before, understanding it better now, hence the new phrasing of the question. – TuanDT Jun 08 '16 at 12:32
  • 5
    You have a Dirichlet distribution with a uniform prior over its parameters? Are you sure you want to torture yourself with a uniform prior? You might want to take a look at this post and its comments. http://andrewgelman.com/2009/04/29/conjugate_prior/ – alberto Jun 08 '16 at 13:07

1 Answers1

3

How about just running some MCMC?

library(MCMCpack)
a <- c(2,3,5)
x <- rdirichlet(100,a)
logposterior <- function(a,x) 
  sum(log(ddirichlet(x,a)))
chain <- MCMCmetrop1R(logposterior,c(1,1,1),x=x)
library(coda)
plot(chain)

enter image description here

Jarle Tufto
  • 7,989
  • 1
  • 20
  • 36