Distribution of *conditional* frequencies when frequencies follow a Dirichlet distribution

Question

Context: we have a large number of individuals characterized by two binary traits; call these $T$ with values $\{0,1\}$, and $T'$ with values $\{0',1'\}$. So there are four types of individuals: $00'$, $01'$, $10'$, $11'$, which appear in the population with unknown relative frequencies $f_{00'}$, $f_{01'}$, $f_{10'}$, $f_{11'}$ summing up to one.

Suppose that our degree of belief about these frequencies (assumed continuous) is expressed by a Dirichlet distribution with parameters $(Aa_{00'}, Aa_{01'}, Aa_{10'}, Aa_{11'})$, the $a$s summing up to one: $$\mathrm{p}[f_{00'}, f_{01'}, f_{10'}, f_{11'} \mid A, (a_{00'}, a_{01'}, a_{10'}, a_{11'})] \propto \prod_{i=0}^1\prod_{j'=0'}^{1'} f_{ij'}^{A a_{ij'}-1}\;\delta\bigl({\textstyle\sum_{ij'}}f_{ij'}-1\bigr).$$

We can also consider the marginal frequencies of individuals having trait $T'$ only, for example: $f_{0'} \equiv f_{00'} + f_{10'}$ and $f_{1'} \equiv f_{01'} + f_{11'}$. Owing to the "aggregation" property of the Dirichlet distribution (Kotz & al 2000, also Basu & al 1982), these marginal frequencies also have a Dirichlet distribution with parameters $\bigl(A(a_{00'} + a_{10'}), A(a_{01'} + a_{11'})\bigr)$ (a Beta distribution).

Question: Consider now the conditional frequencies of trait $T$ given $T'$, for example $$f_{1\mid 0'} \equiv \frac{f_{10'}}{f_{00'}+f_{10'}}.$$ What distribution expresses our degree of belief about such a conditional frequency, given the context above?

While I sit down and calculate (or sample), I'd be grateful for any literature or calculation hints on this. Thank you!

Additional motivation: For inference about sequential data, like for example text, speech, genes, some literature express the degree of belief about conditional frequencies $f_{i\mid j}$ (of, say, one word given the previous one) with a Dirichlet distribution (e.g. MacKay & al 1995): $$\mathrm{p}[f_{i \mid j} \mid A, (a_{i\mid j})] \propto \prod_{i} f_{i\mid j}^{A a_{i\mid j}-1}\;\delta\bigl({\textstyle\sum_{i}}f_{i\mid j}-1\bigr), \qquad\text{for every }j.$$ This approach is different from using a Dirichlet distribution for the joint frequencies $f_{ij}$, and I wonder how different is the distribution for the conditional frequencies that we obtain by assuming Dirichlet for the joint frequencies instead, as in my question above.

References:

– Basu, de Bragança Pereira: On the Bayesian analysis of categorical data: the problem of nonresponse (1982) https://doi.org/10.1016/0378-3758(82)90004-0, §§ 3–4.

– Kotz, Balakrishnan, Johnson: Continuous Multivariate Distributions. Vol. 1 (2nd ed. Wiley 2000), §49.1.

– MacKay, Peto: A hierarchical Dirichlet language model (1995) https://doi.org/10.1017/S1351324900000218, https://pdfs.semanticscholar.org/01fa/57bd91f731522c861404d29e4604ba6ac6d3.pdf.

Answer 1

Let $\gamma_{ij} \sim \texttt{gamma}(A a_{ij})$ independently. Recall that construction of the Dirichlet distribution as a normalization of gamma random variables. Then the Dirichlet distribution we start with is equal in distribution to $$ \left(\frac{\gamma_{00}}{\sum \gamma_{ij}}, \ldots, \frac{\gamma_{11}}{\sum \gamma_{ij}}\right) \sim \texttt{dirichlet}(A a_{00}, \ldots, Aa_{11}). $$ But now we have, for example, $$ f_{1 | 0} = \frac{f_{10}}{f_{00} + f_{10}} \stackrel{d}{=} \frac{\gamma_{10}}{\gamma_{00} + \gamma_{10}}. $$ And we can get a similar expersion for $f_{0|0} = \gamma_{00}/(\gamma_{00} + \gamma_{10})$, whence it follows that $$ (f_{0|0}, f_{1|0}) \sim \texttt{dirichlet}(A a_{00}, A a_{10}). $$ Easy.

To answer your question about how the conditional model varies from the joint model: if you specify Dirichlet priors on the conditionals of the first component given the second, and a very particular Dirichlet prior on the marginals of the second, then you get a Dirichlet prior for the joint. So you get a little bit more flexibility by using Dirichlet priors on the marginals and conditionals separately; the most general sort of prior you get this way is a special case of what has been called a hyper-Dirichlet, or Dirichlet-tree prior. The hyper-Dirichlet is also conjugate to multinomial sampling; in addition to the Dirichlet, it also contains all stick-breaking constructions you can get from the Beta distribution, and many other possibilities as well.

Answer 2

Edit: (I decided not to delete this answer, as it contains a proof of the distributive property of the Dirichlet Distribution. I have however now managed to answer the original post, which I put in a separate answer)

I generally think about these problems using the fundamental theorem of calculus. You make a new variable (let's call it a), write down an integral for the probability that $a<A$, $P(a<A)$, differentiate wrt A (this will usually use the fundamental theorem of calculus, or in more complex cases, Leibniz rule for differentiating under the integral sign, a generalisation of the former), and the result gives you the pdf for a evaluated at A, or $p(A)$. Let's see how this works, to derive an expression for $a=q_{1}+q_{2}$, or what you refer to as the associative property of the Dirichlet Distribution.

$P(a<A)=\int_{0}^{A}dq_{1}\int_{0}^{A-q_{1}}dq_{2}\int_{0}^{1-q_{1}-q_{2}}dq_{3}D(q_{1}, q_{2}, q_{3}, 1 - q_{1} - q_{2} - q_{3};\alpha)$

Writing this as $\int_{0}^{A}f(A,q_{1})dq_{1}$ and differentiating, one obtains:

$f(A,A) + \int_{0}^{A}\frac{\partial}{\partial A}f(A,q_{1})dq_{1}$

Because

$f(A,q_{1})=\int_{0}^{A-q_{1}}dq_{2}\int_{0}^{1-q_{1}-q_{2}}dq_{3}D(q_{1},q_{2},q_{3}, 1-q_{1}-q_{2}-q_{3})$

we can see that $f(A,A)=0$, so we only have to worry about

$\int_{0}^{A}dq_{1}\frac{\partial}{\partial A}\int_{0}^{A-q_{1}}dq_{2}\int_{0}^{1-q_{1}-q_{2}}dq_{3}D(q_{1},q_{2},q_{3}, 1-q_{1}-q_{2}-q_{3})$

which is a relatively simple case of the fundamental theorem of calculus, basically replace all instances of $q_{2}$ with $(A-q_{1})$

$\int_{0}^{A}dq_{1}\int_{0}^{1-q_{1}-(A-q_{1})}dq_{3}D(q_{1}, a-q_{1}, q_{3}, 1 - q_{1} - (A-q_{1})-q_{3};\alpha)$

which, explicitly now subbing in the function form of the Dirichlet Distirbution, is given by

$\frac{1}{B(\alpha)}\int_{0}^{A}dq_{1}q_{1}^{\alpha_{1}-1}(A-q_{1})^{\alpha_{2}-1}\int_{0}^{1-A}dq_{3}q_{3}^{\alpha_{3}-1}(1-A -q_{3})^{\alpha_{4}-1}$

This is now a product of integrals rather than a double integral. The first is solved with the substitution $v=\frac{q_{1}}{A}$ and the second with $u=\frac{q_{3}}{1-A}$. The first gives $A^{\alpha_{1}+\alpha_{2}-1}B(\alpha_{1}, \alpha_{2})$ and the second gives $(1-A)^{\alpha_{3}+\alpha_{4}-1}B(\alpha_{3},\alpha_{4})$

Putting these together with the original normalisation constant, you get the desired result.

However for $a=\frac{q_{1}}{q_{1}+q_{2}}$, I've hit a stumbling block which might turn out to be very hard to resolve and mean this doesn't have a closed form solution, or perhaps just that this method doesn't work. The sticking point, is that in order to define the region in 2-d space such that $0< \frac{q_{1}}{q_{1}+q_{2}}<1$, is actually quite hard. For example, if $q_{1}=1$ and $A=\frac{1}{4}$, no (permitted) $q_{2}$ can satisfy $\frac{1}{1+q_{2}}<\frac{1}{4}$.

In fact, if $A>\frac{1}{2}$, this is always possible but if $A<\frac{1}{2}$, then $q_{1}<\frac{A}{1-A}$ is needed. But the condition $q_{1}<\frac{A}{1-A}$ isn't enough, because when $A>\frac{1}{2}$, this constraint is too weak, $q_{1}$ needs to be smaller than a number which is greater than 1, so actually the constraint is $q_{1}<min(1, \frac{A}{1-A})$, and this cannot be differentiated, so I don't see how to differentiate under the integral sign.

Answer 3

Right, I've finally solved this, and decided to add another answer rather than edit my original one as I've edited it a lot of times and it's getting messy.

From the joint distribution $P(\underline{q})=\frac{q_{1}^{\alpha_{1}-1}q_{2}^{\alpha_{2}-1}q_{3}^{\alpha_{3}-1}q_{4}^{\alpha_{4}-1}}{B(\underline{\alpha})}$

defined for $0 \leq q_{i} \leq 1 \forall i$ and $q_{1}+q_{2}+q_{3}+q_{4}=1$

and we define a new variable $a=\frac{q_{1}}{q_{1}+q_{2}}$ and wish to derive its distribution. We do this by finding an expression for $p(a<A)$ and differentiating w.r.t A. This will give us the pdf for a evaluated at A. Constructing this expression for the joint distribution is the hardest part (and the differentiation process is a little messy and requires care).

In particular, to construct $p(a<A)$, we really have to think about all of the regions in $q$-space for which $a<A$, as we have to integrate over all of this space to find $p(a<A)$. This was the part I was struggling with in my previous answers. We can solve this two ways, they equate to the same:

First

$0 \leq \frac{q_{1}}{q_{1}+q_{2}} \leq A$

We rearrange in terms of $q_{1}$ and find that $q_{1}<q_{2}\left(\frac{A}{1-A}\right)$

We also know that $q_{1}+q_{2} \leq 1$ and thus $q_{1} \leq 1 -q_{2}$. As both of these need to be true, we can write $q_{1} \leq min(1-q_{2}, q_{2}\left(\frac{A}{1-A}\right))$

Looking at this, we see that in general, when $q_{2}$ is large, $1-q_{2}$ is going to be the lesser of the two terms and vice versa when $q_{2}$ is small. The crossover occurs when $1-q_{2}=q_{2}\left(\frac{A}{1-A}\right)$ or $q_{2}=1-A$. Because $0\leq A \leq 1$, we know that this cross-over will occur for $0\leq q_{2} \leq 1$, i.e. in the allowed range, so will be important.

For $q_{2} < 1-A$, $q_{1}<q_{2}\left(\frac{A}{1-A}\right)$ and

for $q_{2}>1-A$, $q_{1}<1-q_{2}$

Thus our $(q_{1},q_{2})$ integration limits have to look like: $\int_{0}^{1-A}dq_{2}\int _{0}^{q_{2}\left(\frac{A}{1-A}\right)}dq_{1} + \int_{1-A}^{1}dq_{2}\int _{0}^{1-q_{2}}dq_{1}$

Second (note, these have to end up being equivalent)

We start from the same place, namely

$0 \leq \frac{q_{1}}{q_{1}+q_{2}} \leq A$ but re-arrange in terms of $q_{2}$ such that $q_{2}>\left(\frac{1-A}{A}\right)q_{1}$. Also $q_{1}+q_{2}<1$ and thus $q_{2}<1-q_{1}$, or, put together

$\left(\frac{1-A}{A}\right)q_{1}< q_{2} < 1-q_{1}$

Looking at this, we see that if $q_{1}$ becomes large enough, the upper limit for $q_{2}$ will be smaller than the lower limit, so let's check when that crossover occurs.

$1-q_{1}=\left(\frac{1-A}{A}\right)q_{1}$, which happens at $q_{1}=A$, so yes, it happens in the region we care about. Thus as well as having the above constraint on $q_{2}$, we know $q_{1}<A$. So we can write the integral limits as

$\int_{0}^{A}dq_{1}\int_{\left(\frac{1-A}{A}\right)q_{1}}^{1-q_{1}}dq_{2}$

It's not obvious to me, even when having these two forms in front of me, why they are the same. The latter one is easier to work with however, so I'm going to use the latter.

We then integrate over all $q_{3}$ values for which $q_{1}+q_{2}+q_{3}<1$ and then finally replace $q_{4}$ with $1-q_{1}-q_{2}-q_{3}$ as it's uniquely determined. Thus the term we need to differentiate is given by:

$\int_{0}^{A}dq_{1}\int_{\left(\frac{1-A}{A}\right)q_{1}}^{1-q_{1}}dq_{2}\int_{0}^{1-q_{1}-q_{2}}\frac{dq_{3}}{B(\underline{\alpha})}q_{1}^{\alpha_{1}-1}q_{2}^{\alpha_{2}-1}q_{3}^{\alpha_{3}-1}(1-q_{1}-q_{2}-q_{3})^{\alpha_{4}-1}$

We first write this as $\int_{0}^{A}f(q_{1},A)dq_{1}$, so when we differentiate w.r.t A, we get (using Leibniz rule for differentiating under the integral sign)

$\int_{0}^{A}\frac{d}{dA}f(A,q_{1})dq_{1} + f(A,A)$

because we can write

$f(A,q_{1})=\int_{\left(\frac{1-A}{A}\right)q_{1}}^{1-q_{1}}g(q_{2})dq_{2}$

it turns that that $f(A,A)$ is given by

$\int_{1-A}^{1-A}g(q_{2})dq_{2}=0$

so we get a nice simplification and only have to calculate

$\int_{0}^{A}\frac{d}{dA}f(A,q_{1})dq_{1}$

given by (essentially replace all $q_{2}$ with $q_{1}\left(\frac{1-A}{A}\right)$ and then also take the derivative of $q_{1}\left( \frac{1-A}{A}\right)$ wrt A)

$\int_{0}^{A}q_{1}\cdot -1 \cdot \frac{d}{dA}\left(\frac{1-A}{A}\right)\int_{0}^{1-q_{1}-q_{1}\left(\frac{1-A}{A}\right)}\frac{dq_{3}}{B(\underline{\alpha})}q_{1}^{\alpha_{1}-1}\left(q_{1}\left(\frac{1-A}{A}\right)\right)^{\alpha_{2}-1}q_{3}^{\alpha_{3}-1}(1-q_{1}-q_{1}\left(\frac{1-A}{A}\right)-q_{3})^{\alpha_{4}-1}$

which tidies up to

$\left(\frac{1-A}{A}\right)^{\alpha_{2}-1}\frac{1}{A^{2}}\int_{0}^{A}dq_{1}q_{1}^{\alpha_{1}+\alpha_{2}-1}\int_{0}^{1-\frac{q_{1}}{A}}\frac{dq_{3}}{B(\underline{\alpha})}q_{3}^{\alpha_{3}-1}(1-\frac{q_{1}}{A}-q_{3})^{\alpha_{4}-1}$

if you're still with me, we're now basically there. We do the $q_{3}$ integral with the substitution $u=\frac{q_{3}}{1=\frac{q_{1}}{A}}$, to get $(1-\frac{q_{1}}{A})^{\alpha_{3}+\alpha_{4}-1}B(\alpha_{3}, \alpha_{4})$

That leaves us needing to compute the integral $\int_{0}^{A}dq_{1}q_{1}^{\alpha_{1}+\alpha_{2}-1}(1-\frac{q_{1}}{A})^{\alpha_{3}+\alpha_{4}-1}$, which is done using the substitution $u=\frac{q_{1}}{A}$ to obtain $A^{\alpha_{1}+\alpha_{2}}B(\alpha_{1}+\alpha_{2}, \alpha_{3}+\alpha_{4})$

Putting this together, we have $(1-A)^{\alpha_{2}-1}A^{\alpha_{1}-1}\frac{B(\alpha_{1}+\alpha_{2}, \alpha_{3}+\alpha_{4})B(\alpha_{3}, \alpha_{4})}{B(\underline{\alpha})}$

and I'll leave to you to verify that the beta functions do indeed cancel to give $\frac{1}{B(\alpha_{1},\alpha_{2})}$, and so the distribution we're left with, we can recognise as normalised.

And thus, somewhat remarkably, the relative prevalence the first event being 0 and 1 does not have any bearing on how likely the second event is, given the first event. In Dirichlet inference, your $\alpha$s will usually be counts of events (combined with prior constants). If you want to know the relative probabilities of all four outcomes, you can use their respective counts and do Dirichlet inference and if you want to know the conditional probabilities, you can throw away the counts relating to first events that didn't go the right way, and then do binomial inference. The two are entirely compatible with one and other.

Now that I've done the analytics, I'd be interested/relieved to see your numerics give the same.