Estimating the actual distance by using a series of measurements in between An object is moving along a single axis, back and forth. We take $N$ snapshots of the object. For all pairs of snapshot positions $i$ and $j$, we can make signed measurements of distance $x_{i,j} = -x_{j,i}$, but only once for each pair. Each of these measurements has some normally distributed error, and we can estimate its variance $\sigma^2_{i,j}$ for each measurement.
What would be the best estimate of the actual distance between points $1$ and $N$?
It's simple to just take $x_{1,N}$, but we obviously have more information than that. I was thinking of calculating all possible paths from $1$ to $N$, summing distances and variances along those paths, and then taking the best one. I believe I could write a program to do this in $O(n^2)$, but it doesn't feel like it combines all the information -- it just takes the 'shortest' path and ignores the rest.
As always, I have no idea which tags to use on this Stack Exchange so I'm open to suggestions.
 A: Let $z_i$ be the true positions for $i=1,2,\ldots,N.$ For each $i=1,2,\ldots, N-1$ define
$$\beta_i = z_{i+1}-z_i$$
to be the true signed distance.  In these terms, $x_{i,j}$ measures $\beta_i + \beta_{i+1} + \cdots + \beta_{j-1}$ with a precision of $1/\sigma^2_{i,j}.$ This is a linear model and Weighted Least Squares is a good choice for estimating its parameters $\beta_i.$

Specifically, the model matrix $X$ has one row for each distance measurement and $n-1$ columns for the $n-1$ betas, indexed by $1$ through $N-1.$. The row for the measurement of $x_{i,j}$ has ones in positions between $i$ and $j$ inclusive and otherwise has zeros. Putting the measurements $x_{i,j}$ into a vector $\mathbf{x}$ in parallel with the rows of $X,$ and putting the variances $\sigma_{i,j}^2$ into another such vector $\sigma,$ the model is
$$\mathbf{x} = X\beta + \epsilon \circ \sigma$$
where $\epsilon$ is a vector of independent random errors of mean zero and unit variance.  "$\circ$" denotes the component-by-component scaling of those vectors by the components of $\sigma.$
Software typically fits such a model using a function that accepts $X,$ $\mathbf{x},$ and $\sigma$ as its arguments.  The function outputs (at a minimum) the estimated values of $\beta,$ written $\hat\beta,$ and a covariance matrix $\hat V$ of those estimates (which is essential for obtaining standard errors).  In particular, the estimated distance between $z_1$ and $z_N$ is $\hat\beta_1+\cdots+\hat\beta_{N-1} = \mathbf{1}_{N-1}^\prime \hat\beta$ where $\mathbf{1}$ is the $N-1$-vector of ones, whence its estimation variance is
$$\operatorname{Var}(\mathbf{1}_{N-1}^\prime \hat\beta) = \mathbf{1}_{N-1}^\prime \hat V \mathbf{1}_{N-1}$$
(which simply is the sum of all the entries of $\hat V$).  The square root of this number is the standard error of the distance estimate.

For example, I generated a random sequence of $N=200$ numbers on the interval $[0,1]$ and added independent Gaussian noise to each of their differences, using variances ranging from $0.0005$ to $0.055$ and averaging $0.01.$  (The associated standard deviations $\sigma_{i,j}$ were therefore, on average, one-third the distance moved between each snapshot: that's a pretty big measurement error.)
I used weighted least squares to fit the $N-1=199$ parameters $\beta_1, \ldots, \beta_{199},$ with the following results for the signed distance from $z_1$ to $z_{200},$ which is estimated as $\hat\beta_1 + \cdots + \hat\beta_{N-1}$ (the hats represent the individual parameter estimates):
      Estimate         Actual Standard error              Z 
   0.491944302    0.502549159    0.008167935   -1.298352188 

Z is the number of standard errors between the estimate and the true distance.  Notice that the standard error of the estimate is only about $1/120$ times the standard deviation of each individual measurement: this reflects the use of all $\binom{200}{2} = 19900$ measurements in the estimate.
For more details, check out this R code to simulate data, fit this model, run some diagnostics, and display the results.
#
# Generate data in a data frame `df`.
#
n <- 200
sigma <- 0.1

set.seed(17)
z <- runif(n)
x <- outer(z, z, `-`)
j <- rep(seq_len(n), n)[lower.tri(x)]
i <- rep(seq_len(n), each=n)[lower.tri(x)]
x <- x[lower.tri(x)]
sigma2 <- rgamma(length(c(x)), 3, 3) * sigma^2
df <- data.frame(dz = x,
                 x = x + rnorm(length(sigma2), 0, sqrt(sigma2)),
                 sigma2 = sigma2,
                 i=i,
                 j=j)
#------------------------------------------------------------------------------#
#
# Create a model matrix.
#
n <- with(df, max(c(i,j)))
zero <- rep(0, n-1)
B <- t(apply(as.matrix(df[c("i", "j")]), 1, function(ij) {
  b <- zero
  b[min(ij):(max(ij)-1)] <- 1
  b
}))
colnames(B) <- paste0("beta.", seq_len(n-1))
X <- cbind(df["x"], B)
#
# Conducted weighted least squares.
#
w <- 1 / df$sigma2
fit <- lm(x ~ . - 1, X, weights=w)
summary(fit)
#
# Compare the estimates to the actual values.
#
beta.hat <- coefficients(fit)
dz <- diff(z)
se <- sqrt(diag(vcov(fit)))
Z <- (beta.hat - dz) / se
Result <- data.frame(Fit = beta.hat,
                     Actual = dz,
                     `Standard Error` = se,
                     `Z score` = Z)
if (n <= 10) print(Result, digits=3)

par(mfrow=c(1,2))
with(Result, {
  qqnorm(Z, main="Z scores")
  abline(0:1, col="Red", lwd=2)
  
  H <- hist(Z, plot=FALSE)
  plot(H, ylim=c(0, max(c(H$density, 1/sqrt(2*pi)))),
   freq=FALSE, main="Histogram of Z scores", col="#f0f0f0")
  curve(dnorm(z), xname="z", add=TRUE, lwd=2, col="Red")
})
par(mfrow=c(1,1))
#
# Follow-on tests.
#
ks.test(Result$Z, pnorm) # Are the Z-scores Normal?
SE <- sqrt(rep(1,n-1) %*% vcov(fit) %*% rep(1,n-1))
c(Estimate = sum(beta.hat), Actual = z[n] - z[1], `Standard error` = SE, 
  Z = (sum(beta.hat) - (z[n] - z[1])) / SE)

