# How do I characterise a distribution, when each measurement comes paired with an erroneous measurement whose identity is ambiguous?

I have a set of measurements $$x_1 ... x_n$$, which I believe are normally distributed, $$x\sim N(\mu, \sigma^2)$$. I would like to estimate their mean and variance.

Here's the catch: Due to the nature of how the measurements were made, each measurement $$x_i$$ actually comes as a pair $$(x_{i_{lower}}, x_{i_{higher}})$$. One member of each of these pairs is the true measurement, the other is the true measurement +/- some error $$E_i = x_{i_{higher}} - x_{i_{lower}}$$. The identity of the members of the pairs is ambiguous (I don't know a priori which measurement is real and which contains the error).

The probability distribution of the size of the errors can obviously be empirically measured. Probably best not to assume the errors have any particular distribution if possible.

Is there a way I can get meaningful summary statistics from these data? Bonus points for a method of estimating the respective probabilities of $$x_{i_{lower}}$$ or $$x_{i_{higher}}$$ being the real measurement, for each $$i$$.

If the errors are independent of the $$x_i$$ and symmetric, so that errors of $$y$$ and $$-y$$ are equally likely, then:

Let $$A$$ be the average and let $$D$$ be the difference of the higher and lower observations. If the error is $$Y$$, then $$D=|Y|$$ and $$Var[Y]=E[Y^2]-E[Y]^2=E[D^2]-0.$$

Similarly, since $$A=X+\frac12Y$$, \begin{align} \mu &= E[A],\\ \sigma &= Var[A]-\frac14E[D^2]. \end{align}

Probably best not to assume the errors have any particular distribution if possible.

With no prior information about the distribution of $$x$$ or $$E$$ other than a belief that the $$x$$ are independent draws from a normal distribution with unknown parameters, the only criterion we can use to assess whether a candidate set of points, $$\hat x$$, is the "real measurement" is how well it fits a normal distribution.

A reasonable approach may be to find the $$\hat x$$ that minimizes the statistic of a goodness of fit test. This can be done by nested optimization. The inner optimization uses the Hungarian algorithm to find the $$\hat x$$ that minimizes the test statistic for a given $$\mu$$ and $$\sigma$$. The outer optimization uses a general-purpose optimizer to find the $$\mu$$ and $$\sigma$$ that minimize the test statistic returned by the inner optimizer.

## Simulation

The following R code demonstrates this approach with simulated data and the Cramér–von Mises test statistic. The process is repeated 15 times.

A function to perform the simulations:

library(RcppHungarian)

f0 <- function(...) {
(seed <- sample(.Machine$integer.max, 1)) set.seed(seed) # simulate data n <- 500 # Sample E from a distribution with an erratic PDF between -1 and 1 f <- splinefun(seq(-1, 1, length.out = 20), runif(20) + (1:20)/40) E <- runif(n*10, -1, 1) E <- E[f(E) > runif(length(E), 0, 2*max(E))][1:n] # the "real" data is in the first column, and the error-shifted data is in the # second x <- cbind(x <- rnorm(n), x + E) # get rough upper and lower bounds on the parameters of x i <- seq(1/n/2, 1 - 1/n/2, 1/n) xx <- sort(x) rsd <- range(TTR::runSD(c(x), n, 0)[n:length(xx)]) lower <- c(mean(xx[1:n]), log(rsd[1])) upper <- c(mean(xx[(n + 1):(2*n)]), log(rsd[2])) # function to be optimized--uses the Hungarian algorithm to find the optimum # set of values that minimizes the Cramer-von Mises test statistic for a given # mean and standard deviation f <- function(params) { y <- pnorm(x, params[1], exp(params[2])) HungarianSolver( pmin( outer(i, y[,1], "-")^2, outer(i, y[,2], "-")^2 ) )$cost
}

# function to perform the Gaussian parameter optimization and return the
# candidate set that globally minimizes the Cramer-von Mises test statistic
normOptim <- function(x) {
# method "L-BFGS-B" with bounds does well at finding a global minimum in the
# presence of multiple local minima
params <- optim(c(mean(x), log(sd(x))), f, lower = lower, upper = upper,
method = "L-BFGS-B")$$par # retrieve the candidate set y <- pnorm(x, params[1], exp(params[2])) A <- A1 <- outer(i, y[,1], "-")^2 A2 <- outer(i, y[,2], "-")^2 bln <- A1 < A2 A[bln] <- A1[bln] A[!bln] <- A2[!bln] pairs <- HungarianSolver(A)$$pairs
out <- integer(n)
out[pairs[,2]] <- 2L - bln[pairs]
out
}

# for each row, identify the column of x that contains x-hat
idx <- normOptim(x)
xhat <- x[cbind(1:n, idx)]
y <- pnorm(sort(x[,1]), mean(x[,1]), sd(x[,1]))
yhat <- pnorm(sort(xhat), mean(xhat), sd(xhat))

list(
seed = seed,
p = mean(idx == 1L),
# Cramer-von Mises test statistic for x and x-hat
T1 = sum((i - y)^2) + 1/12/n,
T2 = sum((i - yhat)^2) + 1/12/n
)
}


Perform the simulation 15 times in parallel and summarize the results.

library(parallel)

cl <- makeCluster(15, type = "PSOCK")
clusterEvalQ(cl, library(RcppHungarian))
data.table::rbindlist(parLapply(cl, 1:15, f0))

#>          seed     p         T1           T2
#> 1   382014738 0.502 0.05428307 0.0006133394
#> 2   908878520 0.506 0.05151929 0.0008763339
#> 3  1583518653 0.482 0.11607714 0.0007895081
#> 4  2066117082 0.472 0.07424786 0.0007267635
#> 5  1626352239 0.500 0.03307783 0.0005275291
#> 6  2119421263 0.486 0.02932221 0.0005932570
#> 7   626401422 0.496 0.04912518 0.0010036447
#> 8    66718651 0.492 0.10045322 0.0006823435
#> 9  1107913526 0.502 0.02884935 0.0007366666
#> 10 1990492963 0.514 0.02569923 0.0010301989
#> 11   50074952 0.534 0.02930599 0.0006010961
#> 12  165422105 0.550 0.02351056 0.0008490050
#> 13  621963233 0.478 0.08712079 0.0011667846
#> 14  129867844 0.484 0.07650852 0.0008400564
#> 15  843946982 0.520 0.03574757 0.0012544344


## Conclusion

In each of the 15 repetitions, $$\hat x$$ results in a much smaller test statistic (column T2) than the test statistic computed on $$x$$ (column T1). However, notice from column p that, in aggregate, $$\hat x_i$$ is no more likely to be $$x_i$$ than $$x_i + E_i$$. With $$2^n$$ possible sets to form $$\hat x$$, it is unsurprising that a large portion the values in $$\hat x$$ are from $$x_i + E_i$$.

There is insufficient information about the distribution of $$x$$ or $$E$$ in order to reliably distinguish $$x_i$$ from $$x_i + E_i$$.