kernel density estimation of a Dirichlet distribution

Question

I have a Dirichlet distribution over the parameters of a multinomial, and I want to estimate its posterior density given some set of evidence. Due to some pecularities of my problem (e.g. presence of "soft evidence", among others), I unfortunately cannot use a simple update counts on the parameters, but I have to resort to sampling techniques.

Once I have collected enough samples for the posterior distribution, I then need to re-estimate its density. Note that I don't necessarily need to estimate a Dirichlet again, the only thing I need is a probability density function. I thought of two possible strategies, one parametric and one non-parametric:

1. Try to re-estimate the Dirichlet parameters (alpha counts) from the samples. It seems theoretically possible (see here for instance), but quite difficult to implement, since there is no closed-form solution.
2. Alternatively, I could give up trying to re-estimate the Dirichlet counts entirely, and resort to kernel density estimations based on my samples.

I'm currently trying to implement the second strategy, based on a kernel density estimation with Gaussian kernels (see e.g. here). I unfortunately don't seem to get the right result, the calculated densities being too high.

As I understand it, the problem might be that the samples of my Dirichlet are not standard multivariate values: namely, they are constrained to be between 0 and 1 and sum up to 1. They are therefore very strongly correlated. For instance, for a Dirichlet with 2 dimensions (aka a Beta distribution), if the first variable has a value 0.6, the second variable is constrained to have a value 0.4. All the other values have a null density.

I am therefore wondering if it is possible to perform a multivariate kernel density estimation for values drawn from a Dirichlet distribution. If yes, how do you think I should do it? Should I use another kernel? If no, do you have any alternative to KDE to suggest?

Thanks!

Answer 1

You probably can do better than both of those. What people often do is average together the full conditional distributions they get while running their Markov chain. If you can do that - and you can, provided your Dirichlet update is conjugate - this is a better strategy. So, for example, suppose I want the posterior density of $\theta$, give the data. I have run a chain sampling $\theta$, $\eta$, using a Gibbs sampler. The full conditional of $\theta$ is $f(\theta | \eta, y)$ where $y$ is my data. Then, use $$\frac 1 M \sum_{m = 1}^M f(\theta | \eta_m, y)$$ to estimate the density, where $m$ indexes the iteration of the chain, with M total iterations. In the big, the ergodic theorem says that point wise this converges to its expectation under $\eta|y$ which is $f(\theta|y)$, i.e. what you want.

Answer 2

You should be able to use KDR to estimate a Dirichlet type distribution, you just need to ensure you capture the constraints correctly.

First off, drop the dimensionality by 1. As you correctly point out, in the 2D case, you have a beta distribution, hence there is only 1 degree of freedom. In this case the correct thing to do would be to estimate the beta shaped pdf of one of the two variables. In higher dimensions, once you have $n-1$ of your probabilities, you know the last one.

Now you only have a boundary problem. For the 2D case, this is just the edges at 0 and 1. If you use the normal Kernel it will "spill over" beyond this which will have 2 effects

1. The support will be incorrect, if you truncate it you won't sum to zero, if you rescale you will over estimate the centre, underestimate the edges.
2. The pdf estimate at the edges will tend to be pulled towards the middle - there will be no values to the left of $0$ or to the right of $1$ to balance the kernel smooth. I.e. imagine that the distribution rises to a point at 1, then the kernel will lower that point by averaging with values at 0.95, biasing the estimate at 1 downwards.

The second of these is difficult (impossible?) to deal with in generality, although there have been many attempts. Given that your reason for choosing this route was for simply implementation, we shall ignore those. The first problem can be mitigated by simply using a kernel that truncates at 0 and 1. One way of doing this is to use a beta distribution as your kernel, with quite large parameters. This approximates a normal away from 0 and 1, but then "bunches up" when at the edges.

Moving on to higher dimensions, in 3D you have support on a right angled triangle from (0,0) - (0,1) - (1,0). By analogy with the above, you could use a dirichlet kernel. This behaves like a normal distribution in the middle of the triangle, but bunches up at the boundaries. The same then clearly applies in higher dimensions.