Prelims
In the article Sequential updating of conditional probabilities on directed graphical structures by Spiegelhalter and Lauritzen they give an approximation to the expectation of a logistic transformed Gaussian random variable $\theta \sim N(\mu, \sigma^2)$. This uses the Gaussian cdf function $\Phi$ in the approximation
$$ \exp(\theta)/(1 + \exp(\theta)) \approx \Phi(\theta \epsilon) $$
for an appropriately chosen $\epsilon$ (in their case they chose $\epsilon = 0.607$). Hence
$$ \mathbb{E} \left [ \exp(\theta)/(1 + \exp(\theta))\right ] \approx \int_{- \infty}^{\infty} \Phi(\theta \epsilon) \phi(\theta | \mu, \sigma^2) d \theta$$
where $\phi$ is a Gaussian pdf function. The integral can be written as
$$ \int_{\infty}^{\infty} \Pr(U < 0 | \theta) \phi(\theta|\mu, \sigma^2) d\theta $$
where $U \sim N(-\theta, \epsilon^{-2})$ and the integral is then simply the marginal $\Pr(U < 0)$. Note that as $\theta \sim N(\mu, \sigma^2)$, we have $U \sim N(-\mu, \sigma^2 + \epsilon^{-2})$. Hence
$$ \mathbb{E} \left [ \exp(\theta)/(1 + \exp(\theta))\right ] \approx \Pr(U < 0) = \Phi(\frac{\mu}{\sqrt{\sigma^2 + \epsilon^{-2}}})$$
We can then use the initial approximation in the reverse direction to get
$$ \mathbb{E} \left [ \exp(\theta)/(1 + \exp(\theta))\right ] \approx \exp(c \mu)/(1 + \exp(c \mu)) $$
where $c = (1 + \epsilon^2 \sigma^2)^{-1/2}$.
Question
My question is, are there any approximations to the expectation of a softmax transformation of Gaussian multivariate variables. In particular, let
$$ \boldsymbol{Z} \sim MVN(\boldsymbol{\mu}, \Sigma) \in \mathbb{R}^{n} $$
Define the $k$ activations for each discrete outcome as
$$ f_i(\boldsymbol{Z}, \boldsymbol{w}_i) = \boldsymbol{w}_i^T \boldsymbol{Z} $$
Finally define our softmax transformed activations as $$ P_i(\boldsymbol{Z}) = \frac{\exp(f_i(\boldsymbol{Z}, \boldsymbol{w}_i))}{\sum_{j=1}^k \exp(f_j(\boldsymbol{Z}, \boldsymbol{w}_j))} $$
What I want is an estimate to the expectation $$ \mathbb{E}[P_i(\boldsymbol{Z})] $$
Note that in the case $k=2$, we have
$$ P_1(\boldsymbol{Z}) = \frac{\exp(f_1(\boldsymbol{Z}, \boldsymbol{w}_1))}{ \exp(f_1(\boldsymbol{Z}, \boldsymbol{w}_1)) + \exp(f_2(\boldsymbol{Z}, \boldsymbol{w}_2))} $$
Therefore
$$ P_1(\boldsymbol{Z}) = \frac{\exp(f_1(\boldsymbol{Z}, \boldsymbol{w}_1) - f_2(\boldsymbol{Z}, \boldsymbol{w}_2))}{ \exp(f_1(\boldsymbol{Z}, \boldsymbol{w}_1)- f_2(\boldsymbol{Z}, \boldsymbol{w}_2)) + 1} $$
and as $f_1(\boldsymbol{Z}, \boldsymbol{w}_1) - f_2(\boldsymbol{Z}, \boldsymbol{w}_2)$ is simply the sum of correlated Gaussian random variables, it is also Gaussian distributed. Hence we can use the initial approximation.
Can we generalise for $k > 2$?