Encontrar la permutación más probable

Question

[Esperando que este sea el sitio correcto de Stackexchange; inspirado en una historia real vista en el trabajo].

Joe tiene un instrumento de medición y $n$ objetos a medir (por ejemplo, una balanza y $n$ pesos). Mide cada uno, obteniendo una lista de medidas $X=\left[x_1 \dots,x_n\right] \in\mathbb R^n$ .

Más tarde, me envía los objetos. Quiero encontrar la correspondencia entre cada objeto y su respectivo valor medido $x_i$ pero Joe se olvidó de numerar los objetos o de ordenarlos de forma que me permitieran encontrar cuál es el $i$ -ésimo objeto. Por lo tanto, vuelvo a medirlos con un instrumento similar, obteniendo una lista de valores $Y=\left[y_1 \dots,y_n\right] \in\mathbb R^n$ .

Si nuestros instrumentos fueran perfectamente precisos, entonces $Y$ sería una permutación de $X$ . Sin embargo, nuestros instrumentos no son perfectos. veracidad tienen imperfecciones precisión . En otras palabras, si medimos el mismo objeto muchas veces, el valor medio de las mediciones repetidas tiende al valor verdadero, pero los resultados tienen desviaciones típicas (conocidas) $\sigma_J$ y $\sigma_I$ (para el instrumento de Joe y para el mío, respectivamente). Por lo tanto, los valores de $X$ serán en general diferentes de los valores de $Y$ .

En el caso límite de que todos los valores sean distintos entre sí (es decir, $\displaystyle\min_{x_i,x_j\in X}\{|x_i-x_j|\}\gg\sigma_J$ y de forma similar $\displaystyle\min_{y_i,y_j\in Y}\{|y_i-y_j|\}\gg\sigma_I$ ), encontrar la permutación correcta (es decir, la correspondencia entre un valor en $X$ y el valor correspondiente en $Y$ ) es trivial. Sin embargo, cuando éste no es el caso, ¿cómo se puede encontrar la permutación más probable a partir de $X$ a $Y$ a partir de los datos disponibles?

Preguntas adicionales: ¿cambia la respuesta si dejo de suponer una veracidad perfecta? ¿Es el caso $\sigma_J=\sigma_I$ ¿Más fácil?

EDITAR Olvidé preguntar: ¿cómo se calcula el probabilidad de una permutación dada, es decir, la probabilidad de que sea la correcta entre el espacio de $n!$ posibles permutaciones? ¿Existe una expresión sencilla (preferiblemente de forma cerrada) para la probabilidad de la permutación óptima (que parece ser la correspondiente a la ordenación de ambos vectores, véase la solución de whuber más abajo - al menos si los errores se distribuyen normalmente)?

EDITAR 2 Según la observación de Aksakal (véanse los comentarios a la pregunta): supongamos que todas las verdadero son estrictamente distintos (el medidas tanto para mí como para Joe, pueden ser valores no distintos debido a un error de medición).

Answer 1

Siempre que los errores de medición sean independientes e idénticamente distribuidos normalmente para cada instrumento, la solución consiste en hacer coincidir los dos conjuntos de mediciones en orden ordenado. Aunque esto es intuitivamente obvio (los comentarios publicados poco después de la publicación de la pregunta exponen esta solución), queda por pruebe eso.

Para ello, dejemos que el primer conjunto de mediciones en orden de clasificación sea $x_1\le x_2\le \cdots \le x_n$ y que el segundo conjunto de medidas en orden de clasificación sea $y_1\le y_2\le \cdots \le y_n.$ Que las distribuciones de error tengan medias y varianzas nulas $\sigma^2$ para el instrumento X y $\tau^2$ para el instrumento Y. (Esta notación me parece un poco más congenial que la subinscripción de la pregunta).

Para encontrar la permutación más probable, resolvemos el problema de máxima verosimilitud. Sus parámetros son (a) el $n$ pesos verdaderos $\theta_i$ correspondientes a los objetos medidos por cada $x_i$ y (b) la permutación $s$ que hace $y_{s(i)}$ la segunda medición del objeto $i.$ En la medida en que la probabilidad depende de $(\theta)$ y $s,$ la probabilidad de estas observaciones es proporcional a la exponencial de

$$\mathcal{L}(\theta,s) = -\frac{1}{2}\sum_{i=1}^n \left(\frac{x_i-\theta_i}{\sigma}\right)^2 + \left(\frac{y_{s(i)}-\theta_i}{\tau}\right)^2.$$

Para cualquier $s,$ esta expresión (y por tanto su exponencial) se maximiza término a término tomando

$$\hat\theta_i = \frac{\tau^2 x_i + \sigma^2 y_{s(i)}}{\sigma^2 + \tau^2}.$$

Para estos valores óptimos de $\theta,$ el valor de $-2\mathcal{L}$ (que deseamos minimizar ) es

$$-2\mathcal{L}(\hat\theta,s) = \frac{1}{\sigma^2+\tau^2}\sum_{i=1}^n \left(x_i - y_{s(i)}\right)^2.$$

Al expandir cada expresión al cuadrado obtenemos (a) una suma de las $x_i^2,$ (b) una suma de $y_{s(i)}^2$ (que es igual a la suma de los $y_i^2$ porque $s$ es una permutación), y (c) los términos cruzados,

$$-2\sum_{i=1}^n x_i y_{s(i)}.$$

En Desigualdad de reordenación establece que tales sumas de productos son maximizado (maximizando así $\mathcal{L}(\hat\theta, s)$ ) cuando el $y_{s(i)}$ están en orden creciente, QED.

---

Este análisis se basa en la hipótesis de normalidad. Aunque puede relajarse, se necesita algún supuesto de distribución, como señala perspicazmente @fblundun en un comentario a la pregunta.

Answer 2

@whuber (+1) ha respondido a la pregunta en su título sobre la búsqueda de la permutación más probable. Mi propósito aquí es explorar brevemente mediante simulación si se puede esperar que esa permutación más probable sea exactamente correcta. A grandes rasgos, el orden para el segundo pesaje es más probable que sea correcto si la diferencia mínima en los pesos de los objetos es mayor que $3\sigma,$ donde $\sigma$ es la desviación típica de la balanza utilizada para realizar el pesaje.

Digamos que hay $n = 6$ objetos con pesos verdaderos $\mu = 10, 20, 30, 40, 50, 60$ medido inicialmente en una balanza con $\sigma=1.$ Entonces los pesos podrían ser como se muestra en el vector x abajo. [Muestreo y cálculos en R.]

set.seed(112)
n = 6; mu = c(10,20,30,40,50,60)
x = numeric(n)
for(i in 1:n) { x[i] = rnorm(1, mu[i], 1) }
x
[1]  9.685768 22.403375 29.281736 38.239389 48.874719 59.280459

Dado que las diferencias más pequeñas (10) en pesos verdaderos son considerablemente mayores que $\sigma=1$ podemos esperar que una segunda pesada, los ponga en el mismo orden que en y .

set.seed(113)
n = 6; mu = c(10,20,30,40,50,60)
y = numeric(n)
for(i in 1:n) { y[i] = rnorm(1, mu[i], 1) }
y
[1] 10.13335 21.37522 30.74872 38.70615 49.44123 58.26755

Así pues, las clasificaciones de los seis objetos en las dos pesadas coinciden:

sum(rank(x)==rank(y))
[1] 6

Sin embargo, si las diferencias entre los pesos verdaderos no son varias $\sigma$ s separados, entonces las filas pueden estar revueltas. (Abajo, de nuevo con espacios de 10, pero $\sigma=4,$ muestras 2 y 3 no son coherentes).

set.seed(112)
n = 6; mu = c(10,20,30,40,50,60)
x = numeric(n)
for(i in 1:n) { x[i] = rnorm(1, mu[i], 4) }
x
[1]  8.743073 29.613500 27.126944 32.957556 45.498875 57.121838

set.seed(113)
n = 6; mu = c(10,20,30,40,50,60)
y = numeric(n)
for(i in 1:n) { y[i] = rnorm(1, mu[i], 4) }
y
[1] 10.53342 25.50089 32.99486 34.82460 47.76492 53.07020

sum(rank(x)==rank(y))
[1] 4

Así pues, para saber cuál puede ser la mejor permutación, habría que tener en cuenta la diversidad de los pesos de los objetos y la precisión de la balanza.

Answer 3

Inspirada en la respuesta de BruceET, he aquí una simulación que calcula la distribución de la diferencia de rango. Mientras que la respuesta de Whuber muestra que la identidad es la permutación más probable, no es la diferencia más probable en el rango (sólo hay una permutación, la identidad, sin cambio en el rango, pero hay muchas permutaciones con la misma diferencia en el rango y estos se suman).

La simulación utiliza alguna distribución aleatoria de los pesos verdaderos (con la línea mu <- c(1:n)*d se distribuyen uniformemente)

El histograma puede aproximarse con una distribución binomial de Poisson, que a su vez puede aproximarse con una distribución gaussiana. La imagen siguiente muestra una comparación entre la simulación (el histograma) y la aproximación (la línea con puntos).

### settings
set.seed(1)
n <- 200
d <- 2
sig <- 1
mu <- runif(n,0,d*n) 
mu <- mu[order(mu)]
#mu <- c(1:n)*d

### simulation of drawing X and Y
### the output is the average difference in rank
simulate <- function(n,d,sig,mu) {
  x <- rnorm(n,mu,sig)
  y <- rnorm(n,mu,sig)
  return(mean(abs(rank(x)-rank(y))))
}

### perform the simulation 10 000 times and plot historgram
s <- replicate(10^4, simulate(n,d,sig,mu))
hist(s, breaks = seq(-1/n,max(s)+1/n,2/n), freq = 0,
     main = "average difference in rank",
     xlim = c(0.5,1))

### compute probabilities that ranks get swapped 1 or 2 places
###    the 2*p*(1-p) relates to the probability that 
###    the swap occurs in X but not in Y or vice versa
p <- 1-pnorm(mu[-1]-mu[-n],0,sig*sqrt(2))
p <- 2*p*(1-p)
p2 <- 1-pnorm(mu[-c(1:2)]-mu[-c(n-1,n)],0,sig*sqrt(2))
p2 <- 2*p2*(1-p2)

### normal approximation of the Poisson-binomial
dvar <- sum(p*(1-p))+sum(p2*(1-p2))
dmu <- sum(p)+sum(p2)

### plotting
ds <- 0:100
lines(ds*(2/n), dnorm(ds,dmu,sqrt(dvar))*(n/2))
points(ds*(2/n), dnorm(ds,dmu,sqrt(dvar))*(n/2), pch = 21 , col = 1, bg = 0, cex = 0.7)