Computation of ratio with Dirichlet distribution

For Dirichlet random variates $x_1,...x_K$ with concentration parameters $\alpha_1...\alpha_K$,

$$y_{i,j}=\frac{x_i}{x_i+x_j}\sim\text{Beta}(\alpha_i,\alpha_j)$$ and $$\frac{x_i}{x_j}=\frac{y_{i,j}}{1-y_{i,j}}\sim\beta'(\alpha_i,\alpha_j)$$ which is the beta prime distibution. The mean exists for $\alpha_j>1$, and the variance exists for $\alpha_j>2$.

$$P(x_i<x_j)=P\bigg(\frac{x_i}{x_j}<1\bigg)=I_{1/2}(\alpha_i,\alpha_j)$$

where $I$ is the regularized incomplete beta function.

pbeta(0.5, 1/sum(1:6), 2/sum(1:6)) #> [1] 0.6677137 

---

Estimates of $\mathbf{\alpha}$

I didn't look into the specifics of DirichReg, but for estimates of $\mathbf{\alpha}$, the moments estimates do well here, from which we can use optim for the MLE:

# load R package library(DirichletReg) # simulate data (6 proportions, 10000 individuals) pct0 = (1:6)/sum(1:6) set.seed(12345) n <- 1e4L data = rdirichlet(n, pct0) data = as.data.frame(data) colnames(data) = paste0('pct',1:6) alpha0 <- colMeans(data) alpha0 <- alpha0*(alpha0[3]*(1-alpha0[3])/var(data[,3]) - 1) tLogData <- t(log(data)) alpha <- exp( optim( log(alpha0), function(alpha) { alpha <- exp(alpha) sum((1 - alpha)*tLogData) + n*(sum(lgamma(alpha)) - lgamma(sum(alpha))) }, method = "L-BFGS-B" )$par ) matrix(c(pct0, alpha0, alpha), 3, 6, 1, list(c("pct0", "MoM", "MLE"), colnames(data))) #> pct1 pct2 pct3 pct4 pct5 pct6 #> pct0 0.04761905 0.09523810 0.1428571 0.1904762 0.2380952 0.2857143 #> MoM 0.04847265 0.09570196 0.1431152 0.1910795 0.2377477 0.2799800 #> MLE 0.04797821 0.09448738 0.1427506 0.1881108 0.2355976 0.2822832 

---

Computing $P(\alpha_i>\alpha_j)$

From the comments:

...ultimately I would like to be able to test for differences between ratios. Say the % correspond to budget shares spent on different types of products (20% on A; 30% on B; 10% on C and 40% on D), we could say that people spent 1.5 more on B than A and 2 times more on D than A - But is 1.5 significantly different from 2?

Comparing the ratio of the proportion spent on B vs A to the ratio of the proportion spent on D vs A is equivalent to comparing the proportion spent on B to the proportion spent on D:

$$P\bigg(\frac{\alpha_D}{\alpha_A}>\frac{\alpha_B}{\alpha_A}\bigg)=P(\alpha_D>\alpha_B)=P(E[x_D]>E[x_B])$$

In what follows, I will use the notation $x_{(i)k}$ to refer to the proportion of class $i$ in the $k^{\text{th}}$ observation out of $n$ total observations.

Recall:

$$y_{(i,j)k}=\frac{x_{(i)k}}{x_{(i)k}+x_{(j)k}}\sim\text{Beta}(\alpha_i,\alpha_j)$$

We can obtain $P(\alpha_i>\alpha_j)$ by numerically integrating the joint posterior:

$$P(\alpha_i>\alpha_j)=\frac{\int_0^{\infty}\int_0^{\alpha_i}\pi(\alpha_i,\alpha_j)\prod_{k=1}^n\frac{y_{(i,j)k}^{\alpha_i-1}(1-y_{(i,j)k})^{\alpha_j-1}}{B(\alpha_i,\alpha_j)}d\alpha_jd\alpha_i}{\int_0^{\infty}\int_0^{\infty}\pi(\alpha_i,\alpha_j)\prod_{k=1}^n\frac{y_{(i,j)k}^{\alpha_i-1}(1-y_{(i,j)k})^{\alpha_j-1}}{B(\alpha_i,\alpha_j)}d\alpha_jd\alpha_i}$$

where $\pi(\alpha_i,\alpha_j)$ is the prior distribution of $\alpha_i$ and $\alpha_j$.

The following R function performs this calculation:

library(cubature) alpha_compare <- function(x, plogab) { # Performs a pairwise comparison of Dirichlet parameters # INPUT: # x: Dirichlet-distributed data, with observations along the rows and # categories along the columns. # plogab: function returning the natural log of the prior distribution # OUTPUT: a ncol(x)-by-ncol(x) matrix, p, where p[i, j] = P(alpha_i > alpha_j) n <- nrow(x) m <- ncol(x) alpha0 <- colMeans(x) # MoM parameter estimation alpha0 <- alpha0*(alpha0[m%/%2]*(1-alpha0[m%/%2])/var(x[,m%/%2]) - 1) # precalculate sum(log(y)) and sum(log(1 - y)) sum_ln_y <- outer(1:m, 1:m, Vectorize(function(i, j) sum(log(x[,i]/(x[,i] + x[,j]))))) p <- diag(NA, m) # posterior likelihood f <- function(a) function(b) exp(n*(lgamma(a + b) - lgamma(a) - lgamma(b)) + (a - 1)*sum_ln_y[i, j] + (b - 1)*sum_ln_y[j, i] + plogab(a, b) - C) g <- function(u = NA) Vectorize(function(a) cubintegrate(f(a), 0, if (is.na(u)) a else u)$integral) for (i in 1:(m - 1)){ for (j in (i + 1):m) { # use a scaling constant for numeric stability C <- n*(lgamma(alpha0[i] + alpha0[j]) - lgamma(alpha0[i]) - lgamma(alpha0[j])) + (alpha0[i] - 1)*sum_ln_y[i, j] + (alpha0[j] - 1)*sum_ln_y[j, i] + plogab(alpha0[i], alpha0[j]) p[i, j] <- cubintegrate(g(), 0, Inf)$integral if (p[i, j] != 0) p[i, j] <- p[i, j]/cubintegrate(g(Inf), 0, Inf)$integral p[j, i] <- 1 - p[i, j] } } p } 

Demonstrated using $\pi(\alpha_i,\alpha_j)=e^{-\alpha_i-\alpha_j}$, first with the values of $\mathbf{\alpha}$ from the OP's simulated data:

set.seed(69916044) n <- 1e4L m <- 6L x <- matrix(rgamma(m*n, rep((1:m)/sum(1:m), each = n)), n, m) x <- x/rowSums(x) colMeans(x) #> [1] 0.05181357 0.09296486 0.14159702 0.18714575 0.23150131 0.29497748 logprior <- function(a, b) -a - b alpha_compare(x, logprior) #> [,1] [,2] [,3] [,4] [,5] [,6] #> [1,] NA 0.0000000 0.000000e+00 0.000000e+00 0.000000000 0.000000e+00 #> [2,] 1 NA 3.088801e-06 0.000000e+00 0.000000000 0.000000e+00 #> [3,] 1 0.9999969 NA 9.619879e-104 0.000000000 0.000000e+00 #> [4,] 1 1.0000000 1.000000e+00 NA 0.000132252 6.671738e-189 #> [5,] 1 1.0000000 1.000000e+00 9.998677e-01 NA 1.758110e-58 #> [6,] 1 1.0000000 1.000000e+00 1.000000e+00 1.000000000 NA 

With a large number of samples and a large separation between the values of $\mathbf{\alpha}$, the likelihoods are unsurprisingly close to 0 and 1.

A second example with a smaller number of samples and less separation between the values of $\mathbf{\alpha}$:

n <- 100L x <- matrix(rgamma(m*n, rep(m + (1:m)/sum(1:m), each = n)), n, m) x <- x/rowSums(x) colMeans(x) #> [1] 0.1557781 0.1542671 0.1676780 0.1816739 0.1702702 0.1703327 alpha_compare(x, logprior) #> [,1] [,2] [,3] [,4] [,5] [,6] #> [1,] NA 0.6131259 0.1542035 0.010199254 0.03471897 0.06467440 #> [2,] 0.3868741 NA 0.1090327 0.005366352 0.01768526 0.04666039 #> [3,] 0.8457965 0.8909673 NA 0.105368051 0.22317977 0.32570371 #> [4,] 0.9898007 0.9946336 0.8946319 NA 0.77107643 0.82310036 #> [5,] 0.9652810 0.9823147 0.7768202 0.228923568 NA 0.60937322 #> [6,] 0.9353256 0.9533396 0.6742963 0.176899641 0.39062678 NA 

You can obtain the exact mean and variance of the ratio of any two elements of a Dirichlet distributions (assuming that the mean and variance exists which depends on the values of the parameters).

I'm going to take the lazy way out and use Mathematica.

First define the distribution of the ratio of the first two elements:

dist = DirichletDistribution[{a1, a2, a3, a4, a5, a6}]; dist12 = TransformedDistribution[x1/x2, {x1, x2, x3, x4, x5} \[Distributed] dist]; 

Now determine the pdf of the ratio:

pdf = PDF[dist12, r] 

$$\frac{r^{a_1-1} (r+1)^{-a_1-a_2} \Gamma (a_1+a_2)}{\Gamma (a_1) \Gamma (a_2)}$$

for $r>0$ and 0 elsewhere.

The mean and variance when those exist are

Mean[dist12] 

$$a_1/(a_2-1)$$

Variance[dist12] 

$$\frac{a_1 (a_1+a_2-1)}{(a_2-2) (a_2-1)^2}$$

The mean will only exist if $a_2>1$ and the variance will only exist if $a_2>2$. Your example has both $a_1$ and $a_2$ less than 1 so neither the mean nor variance exist in that case.