Analytical solution to the covariance between a continuous and a categorical variable Let $X$ be a continuous variable with mean $\mu$ and $Y$ be a categorical variable with event probability vector $\mathbf{p}$. I am trying to calculate $\operatorname{Cov}(X, Y)$.
I have the solution if $\mathbf{p} = p_1$, i.e., $Y$ is binary. Using the law of total expectation:
$$
\begin{align}
\operatorname{Cov}(X, Y)
&= E(XY) - E(X)E(Y) \\
&= E(XY | Y = 0) p_0 + E(XY | Y = 1) p_1 - \mu p_1 \\
&= E(X0 | Y = 0) p_0 + E(X1 | Y = 1) p_1 - \mu p_1 \\
&= E(0 | Y = 0) p_0 + E(X | Y = 1) p_1 - \mu p_1 \\
&= 0 + E(X | Y = 1) p_1 - \mu p_1 \\
&= E(X | Y = 1) p_1 - \mu p_1
\end{align}
$$
The operations above assume $E(Y) = p_1$, $p_0 = 1 - p_1$ and $E(0 | Y = 0) p_0 = E(0) p_0 = 0 p_0 = 0$.
However, if $Y$ is polytomous, I think I am making some mistake in my calculations, because when I simulate data the observed covariance is very different from the analytical solution. For $\mathbf{p} = \{p_1, p_2\}$, I have
$$
\operatorname{Cov}(X, Y) = E(XY) - E(X)E(Y),
$$
where
$$
E(X)E(Y) = 
\begin{bmatrix}
\mu  p_1 \\
\mu  p_2
\end{bmatrix}
$$
and
$$
\begin{align}
E(XY) 
&= E(XY | Y = 0) p_0 + E(XY | Y = 1) p_1 + E(XY | Y = 2) p_2\\
&= E(X0 | Y = 0) p_0 + E(X1 | Y = 1) p_1 + E(X2 | Y = 2) p_2\\
&= E(0 | Y = 0) p_0 + E(X | Y = 1) p_1 + 2E(X | Y = 2) p_2\\
&= 0 + E(X | Y = 1) p_1 + 2E(X | Y = 2) p_2\\
&= E(X | Y = 1) p_1 + 2E(X | Y = 2) p_2
.\end{align}
$$
For a general case of $i$ categories, this generalizes to
$$
\operatorname{Cov}(X, Y) = \sum_{i} i E(X | Y = i) p_i.
$$
I can already see a problem here because the solution above is a scalar, whereas the $E(X)E(Y)$ above is a vector, so I can't just subtract one from the other. In any case, if I set $\mu = 0$, I get rid of $E(X)E(Y)$ but the results still don't match, which makes me think $E(XY)$ is problematic by itself.
I guess I am also getting a bit lost in the fact that for a binary $Y$ the categories are 0 and 1, whereas a polytomous $Y$ has categories 1, 2, 3... (no zero).
 A: This may not be the most straightforward solution, but intuitively it worked for me, and my simulations also show this works.
So if $Y$ is polynomial, it can be broken down into several binary variables. For example, if $\mathbf{p} = \{p1, p2\}$ instead of working out $Cov(X, Y)$, we can work out $\operatorname{Cov}(X, Y_1)$ and $\operatorname{Cov}(X, Y_2)$ separately, where $Y_1$ and $Y_2$ are Bernoulli-distributed with parameters $p_1$ and $p_2$, respectively. So, for $i = 1, 2$:
$$
\begin{align}
\operatorname{Cov}(X, Y_i)
&= E(XY_i) - E(X)E(Y_i) \\
&= E(XY_i | Y_i = 0) (1 - p_i) + E(XY | Y_i = 1) p_i - \mu p_i \\
&= E(X0 | Y_i = 0) (1 - p_i) + E(X1 | Y_i = 1) p_i- \mu p_i \\
&= E(0 | Y_i = 0) (1 - p_i) + E(X | Y_i = 1) p_i - \mu p_i \\
&= 0 + E(X | Y_i = 1) p_i - \mu p_i \\
&= E(X | Y_i = 1) p_i - \mu p_i
.\end{align}
$$
Note that $Y_i = \{0, 1\}$, so it doesn't matter how the categories are labeled (starting from zero, starting from one, centered in zero, etc.).
This means that $\operatorname{Cov}(X, Y)$ is actually a vector, composed as
$$
\operatorname{Cov}(X, Y) =
\begin{bmatrix}
\operatorname{Cov}(X, Y_1) \\
\operatorname{Cov}(X, Y_2)
\end{bmatrix},
$$
It should be noted that $Y_1$ and $Y_2$ are correlated with covariance equal to $-p_1 p_2$. Moreover, the extension of the above to more categories is straightforward.
