Expectation of scaled dirichlet w.r.t. the aitchinson base

scaled_dirichlet.R

p <- MCMCpack::rdirichlet(1, rep(1, 4))[1,]

N <- 10000
g2 <- matrix(rnorm(length(p) * N), nrow = length(p))^2

pg2 = sweep(g2, MARGIN = 1, STATS = p, FUN = "*")
pg2_norm <- sweep(pg2, MARGIN = 2, STATS = colSums(pg2), FUN = "/")


# ILR Transform to coordinates w.r.t an orthonormal basis of the simplex
pg2_ilr_t <- compositions::ilr(t(pg2_norm))

# Expectation in the transfromed space
E_a_ilr <- colMeans(pg2_ilr_t)
# Backtransform the expectation
E_a <-  compositions::ilrInv(E_a_ilr)
as.numeric(E_a) - p

E <- rowMeans(pg2_norm)
if(length(p)==2) {
  z <- c(p[2]/p[1], p[1]/p[2])
  E_est <- 1/(1 + sqrt(z))
  print(E_est - E)
}



integrand <- Vectorize(function(u) {
  n <- length(p)
  u^(n / 2 - 1) / prod(sqrt(1 + u * (1/p[2:n] - 1))) / (1 + u*(1/p[1] - 1))^1.5
})
int_res <- integrate(integrand, lower = 0, upper = 1)


fixed <- 0.5 /prod(sqrt(p))

int_res$value * fixed - E[1]