3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ It depends on the structure of the matrix $S_{ij}$. Considered as the adjacency matrix of an (edge-weighted undirected) graph, does the graph have a single component or not? Also, do you observe the $S_{ij}$ with any error? How, specifically, do you estimate $\widehat{S}$ from the observations? $\endgroup$ – whuber May 20 '13 at 17:55
  • $\begingroup$ @whuber I edited to answer some of the questions you asked - specifically, to indicate how I estimate $\hat{S}$ and some more info about the struct of the matrix $S_{ij}$. $\endgroup$ – Chris Taylor May 20 '13 at 22:17
2
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks, this was just what I needed. I'll have a crack at the confidence interval stuff myself. $\endgroup$ – Chris Taylor May 21 '13 at 11:24
  • 1
    $\begingroup$ (+1) This is clearly a good approach. But you ought to double-check this work. Starting at the beginning, your expression for the likelihood function causes each observation to appear twice: once for $(i,j)$ and again for $(j,i)$. Thus it seems you have written down the square of the likelihood. The solution will still be the same, but any confidence intervals (from the information matrix) will be wrong. The quadratic program is very rapidly solved. Incidentally, I believe the characterization of positive-definiteness may be equivalent to the graph having a unique connected component :-). $\endgroup$ – whuber May 21 '13 at 16:43
  • 1
    $\begingroup$ @whuber Good point about the double-counting; it should be square-rooted or changed so the products are over $j > i$ for anything that needs an actual likelihood. I just looked over everything else and didn't see anything else wrong, though. $\endgroup$ – djs May 21 '13 at 21:57
  • $\begingroup$ I think the solution is unique. Since the row sums of $A-N$ are zero, the vector ${\bf 1}$ is an eigenvector of $A-N$ with eigenvalue zero. By definition we have ${\bf 1}^T b = 0$ as well, so the zero eigenspace is perpendicular to $b$, and the constraint ${\bf 1}^T S = 0$ fixes a unique solution. $\endgroup$ – Chris Taylor May 23 '13 at 14:24
  • $\begingroup$ @ChrisTaylor I agree that $\mathbf 1$ is an eigenvector with value 0 (thus $A-N$ is never strictly positive definite; oh well...) and is perpendicular to $b$. But how do you know there aren't other eigenvectors with value 0? Consider the degenerate case where $A = N = 0$; the MLE is clearly nonunique but your argument still applies. $\endgroup$ – djs May 23 '13 at 15:49
0
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Did you actually read Dougal's reply? It covers everything you suggest that has any theoretical justification: he describes the "noise," writes the system of linear equations, and shows how to use least squares regression ("polynomial" is unnecessary and inappropriate). $\endgroup$ – whuber May 23 '13 at 22:00
  • $\begingroup$ Yes, Dougal's answer is way more complete and rigorous. Mine is a simple hack for the case where the noise is really small. I'll correct the regression part, this should be "least squares regression"; thanks! $\endgroup$ – eldodo42 May 24 '13 at 17:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.