Estimating values of a sequence from observed differences

Question

I have a sequence of random variables $S_1, S_2 \dots S_N$ that is guaranteed to satisfy

$$S_1 + S_2 + \cdots + S_N = 0$$

I can't observe any of these random variables directly, however I can observe their differences:

$$S_{ij} = S_i - S_j$$

I typically have many observations of each difference $S_{ij}$ - call the observations $s_{ij}^{(1)}\dots s_{ij}^{(n_{ij})}$. Each observation has some error - I'm happy to assume that the error is the same for each observation, and is independent of $i$ and $j$.

I would like to estimate the means of the $S_i$ (and possibly get a confidence interval for the means).

I can form a matrix $\hat{S}_{ij}$ of estimated values for the mean of the differences in the following way:

$$\hat{S}_{ij} = \frac{1}{n_{ij}} \sum_{t=1}^{n_{ij}} s_{ij}^{(t)}$$

The catch is that I don't observe every possible difference, so the matrix has some NA values (put another way, some of the $n_{ij}$ are zero). I can assume that the adjacency matrix of non-NA values has exactly one connected component.

How can I come up with an estimates $\hat{S}_j$ for the mean of the original, non-differenced values?

Answer 1

I agree with @whuber that this is a neat question.

I'm going to assume that each observation $s_{ij}^{(t)} = S_{ij} + \varepsilon_{ij}^{(t)}$, where $\varepsilon_{ij}^{(t)} \sim N(0, \sigma^2)$ (the noise is iid normal).

You can get an estimate for $\sigma^2$ based on the $s_{ij}^{(t)}$ values; I think the weighted average of the sample variance for each $ij$ should probably work okay. But we're going to work out the maximum likelihood estimate for $S$, for which we'll see in a moment that the value of $\sigma$ doesn't actually matter.

Now, of course $S_{ij} = -S_{ji}$. I'm going to assume that your observations already reflect that, so that $n_{ij} = n_{ji}$, and that $s_{ij}^{(t)} = -s_{ji}^{(t)}$. I'll also assume $n_{ii} = 0$, since we know $S_{ii} = 0$.

Then the likelihood function should look like

$$f(S) = I(\sum_i S_i = 0) \cdot \prod_{i=1}^n \prod_{j=i+1}^n \prod_{t=1}^{n_{ij}} N(s_{ij}^{(t)}; S_i - S_j, \sigma^2)$$

Thus, for an arbitrary vector $S$, the likelihood is 0 if the vector doesn't sum to 0, and otherwise is the product of a normal likelihood for each observation with mean at $S_i - S_j$. Taking the product over $j>i$ avoids double-counting observations, but for convenience, I'm going to instead multiply over all $j$ and take the square root below.

The problem of maximizing the likelihood function is then:

$$\max_{S \in \mathbb{R}^n} \sqrt{\prod_{i,j,t} \frac{1}{\sqrt{2 \pi} \sigma} \exp\left( - \tfrac{1}{2 \sigma^2} (s_{ij}^{(t)} - S_i + S_j)^2 \right)} \; \text{subject to } \sum_i S_i = 0,$$

or equivalently:

$$\min_{S \in \mathbb{R}^n} \sum_{i,j,t} (s_{ij}^{(t)} - S_i + S_j)^2 \; \text{subject to } \sum_i S_i = 0.$$

Expanding the square, dropping constant terms, and splitting up the sums, this is

$$\min_{S \in \mathbb{R}^n} \sum_{i,j,t} S_i^2 + \sum_{i,j,t} S_j^2 + 2 \sum_{i,j,t} s_{ij}^{(t)} S_j - 2 \sum_{i,j,t} s_{ij}^{(t)} S_i - 2 \sum_{i,j,t} S_i S_j \; \text{subject to } \sum_i S_i = 0.$$

But then $\sum_{i,j,t} S_i^2 = \sum_{i,j} n_{ij} S_i^2 = \sum_{j,i} n_{ji} S_j^2 = \sum_{i,j,t} S_j^2 = \sum_i \left( \sum_j n_{ij} \right) S_i^2$, where we just swapped the names of $i$ and $j$ for the second equality.

Similarly, $\sum_{i,j,t} s_{ij}^{(t)} S_j = - \sum_{i,j,t} s_{ji}^{(t)} S_j = - \sum_{i,j,t} s_{ij}^{(t)} S_i$.

Thus our problem has become (dividing the objective by 2):

$$\min_{S \in \mathbb{R}^n} \sum_i \left( \sum_j n_{ij} \right) S_i^2 - 2 \sum_i \left( \sum_{j,t} s_{ij}^{(t)} \right) S_i - \sum_{i,j} n_{ij} S_i S_j \; \text{subject to } \sum_i S_i = 0.$$

Now, define $A$ to be the matrix with $A_{ii} = \sum_j n_{ij}$ and off-diagonal elements 0, $N$ to be the matrix with $N_{ij} = n_{ij}$, and $b$ to be the vector with $b_i = \sum_j \sum_{t=1}^{n_{ij}} s_{ij}^{(t)}$. We can rewrite the problem in matrix form as

$$\min_{S \in \mathbb{R}^n} S^T A S - 2 b^T S - S^T N S \; \text{subject to } \sum_i S_i = 0.$$

Dividing the objective by 2 again and doing a trivial rearrangement, we have $$\min_{S \in \mathbb{R}^n} \tfrac{1}{2} S^T \left( A - N \right) S - b^T S \; \text{subject to } \mathbf{1}^T S = 0.$$

This is the standard form of a quadratic program with just one equality constraint.

Now, $A - N$ is the Laplacian of the undirected graph where the edge between $i$ and $j$ has weight $n_{ij}$. It's therefore positive semidefinite, and the multiplicity of 0 in its spectrum is equal to the number of connected components (one, by assumption).

As with all Laplacians, since the row sums are 0, the vector $\bf 1$ is an eigenvector with eigenvalue 0. We also have that $b^T 1 = 0$, since that just sums all the (paired) observations. Thus, ignoring the constraint, the maxima of the objective function form a line. Of course, our constraint is a hyperplane normal to that line, so the constrained maximum is a unique point.

The solution satisfies this linear system:

$$\begin{bmatrix} A-N & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \begin{bmatrix} \hat S \\ \lambda \end{bmatrix} = \begin{bmatrix} b \\ 0 \end{bmatrix}$$

If you run into numerical difficulties, there are many general QP solvers out there.

---

To get a confidence region on $S$, note that the feasible set is a hyperplane, and the likelihood function is (a monotone function of) a quadratic on that hyperplane. You should thus be able to find an ellipsoid about the MLE that contains any arbitrary amount of probability mass. I don't have time to work that out right now, but maybe I'll give it a shot later....

---

If you're not happy with the normal noise model, you can try the same thing for any other distribution. I don't know if it'll come out so (relatively) nicely, though.

Answer 2

A really easy way to figure out the means of the $S_i$'s, probably simpler than Dougal's method (but maybe not as accurate), is as follows:

Because of linearity of expectation, we know that $E[S_i-S_j]=E[S_i]-E[S_j]$. Since we get a lot of samples of $S_i-S_j$ we can get a good estimate for $E[S_i-S_j]$, so we have a good estimate for $E[S_i]-E[S_j]$. Since we get these estimates for a lot of pairs $(i,j)$, we can simply write a system of linear equations and solve them.

The problem with this is that there is some noise involved, so the system of linear equations won't really have a solution. There are various ways to solve this: One is to just keep exactly $N$ equations. This will make sure the system has only one solution. If the noise is small enough, the solution you get will be close to the correct solution. (Note that the error propagates so it can get multiplied by as much as $n$ in this process).

If the noise is a bit too big for this then you can run least squares regression: your variables are the errors, and you're trying to minimize the sum of squares of the errors.