# Problem Deriving Expected Sufficient Statistic

Any exponential family distribution can be expressed as $$p(x|\theta) = g(\theta) f(x) e^{\phi(\theta)^T T(x)} = f(x) e^{\phi(\theta)^T T(x) - A(\theta)}$$ where $$A(\theta) = -\log{g(\theta)}$$.

We know that this integrates to 1 since it's a valid pdf:

$$\int f(x) e^{\phi(\theta)^T T(x) - A(\theta)} dx = 1$$

The expected value of the sufficient statistic is obtained by differentiating both sides with respect to $$\theta$$. Assuming this is a parameter vector, we will need a vector derivative:

\begin{align} \frac{\partial}{\partial \theta_i} \int f(x) e^{\phi(\theta)^T T(x) - A(\theta)} dx &= \frac{\partial}{\partial \theta_i} (1) \\ \int \frac{\partial}{\partial \theta_i} [ f(x) e^{\phi(\theta)^T T(x) - A(\theta)}] dx &= 0 \\ \int p(x|\theta) \left[ \frac{\partial \phi(\theta)^T}{\partial \theta_i} T(x) - \frac{\partial A(\theta)}{\partial \theta_i } \right] dx &= 0 \\ \nabla_i \phi(\theta)^T \int p(x|\theta) T(x) dx - \nabla_i A(\theta) \int p(x|\theta) dx &= 0 \\ \nabla_i \phi(\theta)^T E \left[ T(x) \right] - \nabla_i A(\theta) &= 0\end{align}

But this is not quite correct since the first term is a scalar (dot product between two vectors?) and the second term is a vector.

The expectation should be the vector i.e. I should have something like $$E[T(x)]_i$$ with the vector index $$i$$ on the expected value. What has gone wrong?

I think (but not sure) that the answer should be $$E[T(x)]_i = \frac{1}{\phi'(\theta)} \nabla_iA(\theta)$$

Let $$\theta \in \Theta \subset \mathbb R^m$$ and $$\phi : \Theta \to \mathbb R^d$$ so $$T(x) \in \mathbb R^d$$ too. I'll use without proof the fact that differentiation and integration can be exchanged for exponential families.

First I'm going to treat this as a natural exponential family and see what happens (this is like if $$\phi$$ is just the identity function). I'll use $$\theta$$ as the parameter so right now I have $$p(x|\theta) = \exp(\theta^TT(x) - A(\theta))f(x)$$.

In this case the differentiation is easy and we have $$0 = \int p(x|\theta) \left[T_j(x) - \nabla _jA(\theta)\right]\,\text dx$$ so $$\nabla A(\theta) = \text E(T(X)).$$

So if you're content to work with $$\phi := \phi(\theta)$$ as the parameter, then this is a tidy result.

Now I'll consider what happens if $$m = d = 1$$ to build our intuition for the general case. We'll have $$0 = \int p(x|\theta)\left[T(x) \cdot \phi'(\theta) - A'(\theta) \right]\,\text dx \\ \implies \text E(T(X)) \cdot \phi'(\theta) = A'(\theta)$$ so $$\text E(T(X)) = \frac{A'(\theta)}{\phi'(\theta)}$$ as you thought.

But the situation is more complicated if $$m,d > 1$$. The derivates in particular are more complicated because now $$\phi$$ is vector-valued so we'll get a matrix of first derivatives.

We need $$\frac{\partial}{\partial \theta_j} \phi(\theta)^TT(x).$$ Writing this as a sum, we have $$\frac{\partial}{\partial \theta_j} \phi(\theta)^TT(x) = \sum_{i=1}^d T_i(x) \frac{\partial}{\partial \theta_j}\phi_i(\theta).$$ $$\phi : \Theta\to\mathbb R^d$$ with $$\Theta\subset\mathbb R^m$$ so for each component we have $$\phi_i : \Theta\to\mathbb R$$ with $$i=1,\dots,d$$. Each $$\phi_i$$ is just a typical scalar-valued function so $$\frac{\partial}{\partial \theta_j} \phi_i = (\nabla \phi_i)_j$$.

I'll collect these gradients into a matrix $$\Phi(\theta) \in \mathbb R^{m\times d}$$ where $$\Phi(\theta) = \left[\begin{array}{c|c|c|c}\nabla \phi_1(\theta) & \nabla \phi_2(\theta) & \cdots & \nabla \phi_d(\theta)\end{array} \right].$$

This means $$\frac{\partial}{\partial \theta_j}\phi_i(\theta) = \Phi(\theta)_{ji}$$ so $$\frac{\partial}{\partial \theta_j} \phi(\theta)^TT(x) = \left(\Phi(\theta)T(x)\right)_j.$$

All together this means $$\mathbf 0 = \int p(x|\theta) \left[ \Phi(\theta)T(x) - \nabla A(\theta) \right] \,\text dx \\ \implies \Phi(\theta)\text E [T(X)] = \nabla A(\theta)$$ which is a linear system that we're trying to solve for $$\text E [T(X)]$$. The dimensions in question and the specific properties of $$\phi$$ now will determine how many solutions this system has.

Here's one example. Suppose $$X \sim \mathcal N(\mu, \sigma^2)$$ so $$p(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).$$ Rearranging and expanding the quadratic we have $$p(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi}}\exp\left(\frac{\mu}{\sigma^2}x - \frac{1}{2\sigma^2} x^2 - \left[\frac{\mu^2}{2\sigma^2} + \frac 12 \log \sigma^2\right]\right)$$ which means we'll have $$\theta = (\mu,\sigma^2)$$, $$\Theta = \mathbb R \times (0,\infty)$$, $$T(x) = {x \choose x^2}$$, $$\phi(\theta) = \left(\begin{array}{c}\frac{\theta_1}{\theta_2} \\ \frac{-1}{2\theta_2}\end{array}\right),$$ and $$A(\theta) = \frac 12 \theta_1^2\theta_2^{-1} + \frac 12 \log\theta_2.$$

This means $$\Phi(\theta)= \left(\begin{array}{cc}\theta_2^{-1} & 0 \\ -\theta_1\theta_2^{-2} & \frac 12 \theta_2^{-2}\end{array}\right).$$ Noting that $$\det \Phi(\theta) = \frac 12 \theta_2^{-3} > 0$$ we know this matrix is invertible, and it's just $$2\times 2$$ so the inverse is easy enough to compute: $$\Phi^{-1}(\theta) = 2\theta_2^3 \left(\begin{array}{cc}\frac 12 \theta_2^{-2} & 0 \\ \theta_1\theta_2^{-2} & \theta_2^{-1}\end{array}\right) = \left(\begin{array}{cc}\theta_2 & 0 \\ 2\theta_1\theta_2 & 2\theta_2^2\end{array}\right).$$ We also have $$\nabla A(\theta) = \left(\begin{array}{cc}\theta_1\theta_2^{-1} \\ -\frac 12\theta_1^2\theta_2^{-2} + \frac 12 \theta_2^{-1}\end{array}\right)$$ so $$\text E(T(X)) = \Phi^{-1}(\theta)\nabla A(\theta) \\ = \left(\begin{array}{c}\theta_1 \\ \theta_1^2 + \theta_2\end{array}\right) = {\mu \choose \mu^2 + \sigma^2}$$ which we know to be the first two raw moments of $$X$$. This was a lot of work to get there but it shows the unavoidable role of the matrix $$\Phi$$.

Update: here's the multivariate normal case. Let $$X \sim \mathcal N_n(\mu, \Omega^{-1})$$ where $$\Omega$$ is the precision matrix. I'll use $$\Omega_{(i)}$$ for the $$i$$th row of $$\Omega$$, and for a vector $$v$$ and $$1\leq i\leq j$$ I'll use $$v_{i:j}$$ to be $$(v_i,v_{i+1},\dots,v_j)$$.

The density is $$p(x \mid \mu, \Omega) = \exp\left[-\frac 12 x^T\Omega x + x^T \Omega \mu - \frac 12 \left(\mu^T\Omega\mu - \log\det \Omega\right)\right] \cdot(2\pi)^{-n/2}.$$

The parameters are $$\mu$$ and $$\Omega$$ but I want to deal with these as a flat vector so I'll combine $$\mu$$ with $$\Omega$$ "unrolled" as $$\theta = (\mu, \Omega_{11}, \Omega_{12}, \dots, \Omega_{1n}, \Omega_{22}, \dots, \dots, \Omega_{n-1,n}, \Omega_{nn}).$$ I'm only taking the diagonal and one triangle of $$\Omega$$ since it's symmetric. This means $$\Theta\subset \mathbb R^m$$ where $$m = n + {n+1\choose 2} = 2n + {n\choose 2}$$.

I now need to get this as $$\phi(\theta)^TT(x) - A(\theta)$$. For $$\phi$$ and $$T$$ I'll use the fact that $$x^T\Omega x = 2\sum_{i < j} x_ix_j \Omega_{ij} + \sum_{i=1}^n x_i^2\Omega_{ii}$$ so I'll take $$\phi(\theta) = \left(\Omega\mu, -\frac 12 \Omega_{11}, -\Omega_{12}, \dots, -\Omega_{1n}, -\frac 12 \Omega_{22}, -\Omega_{23}, \dots, \dots, -\Omega_{n-1,n}, -\frac 12 \Omega_{nn}\right)$$ (I'm writing this in terms of $$\mu$$ and $$\Omega$$ instead of $$\theta$$ just because it's way clearer, but these can all be written in terms of $$\theta$$ if desired) and $$T(x) = (x, x_1^2, x_1x_2, \dots, x_1x_n, x_2^2, x_2x_3, \dots,\dots, x_n^2).$$

We can now construct $$\Phi$$. I'll build it as a $$2\times 2$$ block matrix where $$\Phi(\theta) = \left(\begin{array}{c|c}A&B\\\hline C&D\end{array}\right)$$ with $$A$$ being $$n\times n$$ and $$D$$ being $${n+1c\choose 2}\times{n+1\choose 2}$$. $$A$$, $$B$$, and $$D$$ are easy to work out.

For $$A$$, I need $$\frac{\partial \phi_i}{\partial \theta_j}$$ where $$i,j=1\dots,n$$. Since $$\phi_i(\theta) = \Omega_{(i)}^T\mu$$ for this range of $$i$$ values, and $$\theta_j = \mu_j$$, we'll get $$(\nabla \phi_i)_{1:n} = \Omega_{(i)}$$ (using the symmetry of $$\Omega$$ for the fact that $$\Omega_{k\ell} = \Omega_{\ell k}$$). This means that $$A = \Omega$$.

Now for $$i > n$$, $$\phi_i(\theta)$$ is proportional to some $$\Omega_{j\ell}$$ so derivatives w.r.t. $$\mu$$ are zero. This means $$B = \mathbf O$$, the appropriately sized zero matrix.

Next, for $$i,j>n$$ if $$i\neq j$$ then $$\frac{\partial \phi_i}{\partial \theta_j} = 0$$ so $$D$$ is diagonal. If $$i$$ corresponds to a diagonal entry of $$\Omega$$ then we have a derivative of $$-\frac 12$$, otherwise it'll be a derivative of $$-1$$. The first element of the diagonal of $$D$$ corresponds to $$\Omega_{11}$$ so it begins with $$-\frac 12$$. The next value of $$-\frac 12$$ is for $$\Omega_{22}$$ which happens at element $$1+n$$ since all of the first row of $$\Omega$$ has been counted and $$\Omega_{21}$$ is skipped (since it's accounted for by $$\Omega_{12}$$). There are $$n-2$$ remaining elements of $$\Omega_{(2)}$$ to account for, so the 3rd value of $$-\frac 12$$ in $$D$$ is at position $$1+n+(n-1)$$. The next $$-\frac 12$$ is at $$1+n+(n-1)+(n-2)$$ and so on, so the locations of $$-\frac12$$s are given by the sequence $$1+\sum_{j=0}^{i-2}(n-j)$$ with $$i=1,\dots,n$$. As an example, if $$n=5$$ then the $$-\frac 12$$s will be at locations $$1,6,10,13,15$$. As a sanity check, $$D$$ is supposed to have a diagonal of length $${n+1\choose 2}$$, and the final $$-\frac 12$$ is for $$\Omega_{nn}$$. We have $$1+\sum_{j=0}^{n-2}(n-j) = \sum_{j=0}^{n-1}(n-j) \\ = \sum_{j=1}^n j = {n+1\choose 2}$$ so this agrees with that.

Finally, we need $$C$$. I'm going to skip $$C$$ for now until the final step.

For $$\Phi^{-1}$$, $$\Omega$$ and $$D$$ are always invertible so $$\Phi$$ is too, and we can use the formula for the inverse of a $$2\times 2$$ block matrix to find $$\Phi^{-1}(\theta) = \left(\begin{array}{c|c}\Omega^{-1} & \mathbf O \\ \hline -D^{-1}C\Omega^{-1} & D^{-1}\end{array}\right).$$

Now we can get $$\nabla A$$ to finish this off. I'll also partition it into the first $$n$$ elements, which lines up the upper two blocks of $$\Phi^{-1}$$, and the remaining $${n+1\choose 2}$$ which I'll call $$V$$. We have $$A(\theta) = \frac 12 \mu^T\Omega\mu - \frac 12 \log\det \Omega$$ so $$\frac{\partial A}{\partial \theta_{1:n}} = \frac{\partial A}{\partial \mu} = \Omega\mu.$$ This means $$\nabla A(\theta) = {\Omega\mu \choose V}$$ so $$\Phi^{-1}(\theta)\nabla A(\theta) = {\mu \choose -D^{-1}(C\mu - V)}.$$ This shows that $$\text E(T_{1:n}(X)) = \text E(X) = \mu$$ which is nice to see.

For the rest, I won't work it out in general but I'll just show what we get for the elements corresponding to diagonal terms of $$\Omega$$ and off-diagonal terms.

For $$V$$, I'll need derivatives of $$\log\det \Omega$$ w.r.t. $$\Omega_{ij}$$ which I can get via this result: $$\frac{\partial}{\partial \Omega_{ij}} \log\det \Omega = (2 \Sigma - \text{diag}(\Sigma))$$ due to symmetry. This means that for the elements of $$V$$ corresponding to a diagonal element of $$\Omega$$ we'll have $$\frac 12 \mu_i^2 - \frac 12 \Sigma_{ii}$$ while the off-diagonal will be $$\mu_i\mu_j - \Sigma_{ij}$$.

For $$C\mu$$, the diagonal elements contribute $$\mu_i^2$$ because for that row the only non-zero element is $$\mu_i$$ in the $$i$$th column. $$-D^{-1}$$ will multiply the result by $$2$$, and then $$V_i$$ is $$2\left(\mu_i^2 - \frac 12 \mu_i^2 + \frac 12 \Sigma_{ii} \right) \\ = \mu_i^2 + \Sigma_{ii}.$$ For the off-diagonal, you can show that $$C\mu$$ will contribute $$2\mu_i\mu_j$$ so we get $$2\mu_i\mu_j - \mu_i\mu_j + \Sigma_{ij} = \mu_i\mu_j + \Sigma_{ij}.$$

Thus after all this work we have confirmed that $$\text E(X_i) = \mu_i$$ and $$\text E(X_iX_j) = \mu_i\mu_j + \Sigma_{ij}$$.

• Wow. Thanks for an awesome reply! But what about if I want to consider e.g. a multivariate normal distribution in D dimensions. In this case the matrix of derivatives $\Phi(\theta)$ will be DxD meaning I won't be able to invert it explicitly. And yet we should still be able to obtain algebraic solution for E[T(X)] in this same way? Sep 7, 2019 at 22:28
• @user11128 I think in general $\Phi(\theta)^{-1}\nabla A(\theta)$ might be as explicit as you can get if you take this approach unless $\Phi$ happens to have some special structure. I'll go through the multivariate normal example and update with that. I'm also looking more into general conditions for $\Phi$ to be invertible.
– jld
Sep 8, 2019 at 1:01
• @user11128 I just updated with a multivariate normal example, $\Phi$ had enough structure that the inverse was easy to do
– jld
Sep 12, 2019 at 16:48