3
$\begingroup$

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).

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Possible duplicate of Dirichlet posterior $\endgroup$
    – Xi'an
    Jun 8, 2016 at 4:55
  • $\begingroup$ @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. $\endgroup$
    – TuanDT
    Jun 8, 2016 at 12:32
  • 5
    $\begingroup$ 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. andrewgelman.com/2009/04/29/conjugate_prior $\endgroup$
    – alberto
    Jun 8, 2016 at 13:07

1 Answer 1

Reset to default
3
$\begingroup$

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.