N.B.: Sorry for cross-posting from https://stats.stackexchange.com/posts/347804/edit (I realized it was the wrong venue for the question, but couldn't find an easy way to transfer the question here).
So, consider a set of infintiely-differentiable convex functions $f_i: \mathcal X \rightarrow \mathbb R$, where $i$ varies from $1$ to $m$, and suppose we know all the moments of $f_i(x)$ (and of all the derivatives of $f_i$) for all $i$, as $x$ is sampled from some distribution $P$ on a Hilbert space $\mathcal X$.
Question:
What is a low-variance estimate for the quantities
$$\mathbb E_{x \sim P}\left[\frac{\exp(f_i(x))}{\sum_{j=1}^m\exp(f_j(x))}\right] ?$$
I mean to replace the integrands $\exp(...)/\sum_j \exp(...)$ with other random quantities of same expectation (or approx the same), but with smaller variance (the smaller the better).
Particular case: Affine functons. For simplicity, take $f_i(x) \equiv \langle a_i,x\rangle + b_i$, for some vectors $a_1,\ldots,a_m \in \mathcal X$ and scalars $b_1,\ldots,b_m \in \mathbb R$. Note that in this case, the above sought-for quantity can be rewritten in the form
$$\mathbb \nabla_{b_i} R(b),$$ where $R(b) := E_{x \sim P}\left[\log\left(\sum_{j=1}^m\exp(f_j(x))\right)\right]$. People in finance refer to $R$ as "logarithmic log-returns".
Important note: I should precise that i don't want Monte Carlo (or other black-box simulation technique). I need something more principled which exploits the structure of the problem..