AJO_1617_Resumen_elipses_error

hace 191 días por sevillad1617

Elipses de error

 

Cuando trabajamos con puntos en el plano (pares de variables aleatorias), es muy habitual utilizar elipses de error para representar el error cometido.

Ejemplos: Figura 1 en https://www.researchgate.net/figure/238527902_fig1_Fig-1-Structural-map-of-the-North-Andes-Numbers-are-the-locations-of-geological; figuras 13 y 14 en http://www.mdpi.com/1424-8220/13/8/9747/htm

Más concretamente, supongamos que tenemos un punto de coordenadas exactas $(\mu_x,\mu_y)$ y conocemos la imprecisión con las que se miden las coordenadas, concretamente sus desviaciones típicas, $\sigma_x$, $\sigma_y$. Para una sola variable (por ejemplo la $x$), habitualmente tomamos un intervalo \[ [\mu_x-c\cdot\sigma_x , \mu_x+c\cdot\sigma_x] \] que contendrá nuestra medición con una cierta probabilidad; en la fórmula, $c$ es una constante que depende del valor de certidumbre elegido (por ejemplo $c=1$ corresponde a un 68%, $c=2$ a un 95%, etc.)

Cuando tenemos dos variables, si simplemente tomamos un intervalo para cada variable (digamos $c=1$), tenemos un rectángulo (el punto negro es el punto real, ignora los colores por el momento):

 

Pero esta representación tiene dos problemas:

  • Da una falsa impresión de uniformidad, porque las zonas de distintos colores tienen probabilidades muy distintas; por ejemplo la zona roja es más probable que las otras dos, porque solo hace falta un "gran" error en la $y$ (el valor de la $x$ no es difícil que esté cerca de $\mu_x$), y para la zonas azul o verde hacen falta "grandes" errores en ambas variables, algo más improbable.
  • No refleja la correlación entre variables: si la correlación entre $x,y$ es por ejemplo cercana a 1, los valores de la zona azul serían mucho menos probables que los de la zona verde. Situaciones así ocurren: imagina que las coordenadas polares de un punto son $(20,40^g)$, medidas independientes, donde la distancia es muy precisa pero el ángulo no; es fácil ver que la correlación entre $x,y$ será cercana a $-1$ (puedes calcularlo como una propagación no lineal).

Como puedes imaginar a partir del título de esta hoja, representaremos el error en las coordenadas de un punto con una elipse, la llamada elipse de error. Pero, al igual que en el caso de una sola variable, no hay una sola elipse de error, como no hay un solo intervalo de confianza: dependerá del valor de confianza que elijamos. Luego veremos los casos más habituales en la práctica.

 

Nuestro punto de partida es la matriz de varianzas y covarianzas del par de variables:

\[ \begin{pmatrix} \sigma_x^2 & \sigma_{xy} \\ \sigma_{xy} & \sigma_y^2 \end{pmatrix} \]

De manera equivalente, podemos partir de las desviaciones típicas $\sigma_x$, $\sigma_y$ y el coeficiente de correlación lineal entre $x$ e $y$.

Ejemplo: coordenadas polares a rectangulares

Dadas la distancia al origen y el ángulo con el semieje derecho, calculamos las coordenadas del punto buscado. Propagando obtenemos la matrix 2x2 del punto.

Empecemos por una simulación a modo de ilustración: una nube de puntos, dibujada a partir de errores en ángulo y distancia.

vd = 250; sd = 0.010 va = 40; sa = 0.0050 va = va*pi.n()/200; sa = sa*pi.n()/200 # conversión a radianes F(d,a) = (d*cos(a), d*sin(a)) F(vd,va) 
       
(202.254248593737, 146.946313073118)
(202.254248593737, 146.946313073118)
nube = sum([point(F(gauss(vd,sd),gauss(va,sa)),size=1,color="black") for i in range(5000)]) nube.show(aspect_ratio=1,figsize=4) 
       

Hagamos una propagación para obtener la matriz del punto.

Sda = matrix([[sd^2,0],[0,sa^2]]) J = jacobian(F,(d,a))(vd,va) Sxy = J*Sda*transpose(J) Sxy 
       
[ 0.000198648680058493 -0.000135778259710855]
[-0.000135778259710855  0.000286882741859060]
[ 0.000198648680058493 -0.000135778259710855]
[-0.000135778259710855  0.000286882741859060]

Estas son las d.t. y la correlación.

sx = Sxy[0,0]^0.5 sy = Sxy[1,1]^0.5 rxy = Sxy[0,1]/sx/sy sx; sy; rxy 
       
0.0140942782737710
0.0169376132279333
-0.568768039915406
0.0140942782737710
0.0169376132279333
-0.568768039915406
 
       

Cálculo de la elipse

Ahora aprenderemos a calcular la elipse de error a partir de los datos anteriores. Supondremos que está centrada en el origen para desarrollar las fórmulas. Más abajo tienes el código que calcula las fórmulas que explicamos a continuación.

Una elipse queda determinada por tres valores: las longitudes de los semiejes mayor $\sigma_{max}$ y menor $\sigma_{min}$, y el ángulo $\alpha$ de rotación del semieje mayor. Este ángulo lo tomaremos desde el norte (eje positivo de las $y$) en sentido horario, como el acimut.

Dados $\sigma_x$, $\sigma_y$, $\rho$, calcularemos esos tres valores. Para empezar, la ecuación de una elipse de error se deduce de la distribución normal bidimensional, y resulta

\[ \left(\frac{x}{\sigma_x}\right)^2 - 2 \rho \left(\frac{x}{\sigma_x}\right) \left(\frac{y}{\sigma_y}\right) + \left(\frac{y}{\sigma_y}\right)^2 = c^2 \cdot (1-\rho^2) \]

donde $c$ es un parametro que decide la probabilidad asociada al interior de la elipse (análogamente al coeficiente por el que se multiplica $\sigma$ en una distribución unidimensional una vez elegido el valor de probabilidad). Cuando $c=1$ la elipse se llama estándar; la probabilidad de que un punto de la distribución esté en esa elipse es 39.4% aproximadamente. A modo de comparación, el rectángulo $\pm\sigma_x$, $\pm\sigma_y$ en el caso incorrelado corresponde a una probabilidad del 68.3% al cuadrado, es decir, 46.6%. Claramente este valor de probabilidad es muy bajo (sólo dos de cada cinco mediciones caen dentro de la elipse estándar) así que, de la misma manera que en el caso unidimensional, se usa un coeficiente para cambiar la escala. Eso es precisamente lo que hace $c$ en la ecuación. Al final diremos algo más sobre $c$.

 La idea del cálculo que haremos es bastante simple: si giramos las coordenadas hasta que la elipse tenga sus semiejes en posición horizontal y vertical, eso se corresponde con que la covarianza tras el giro sea cero (en otro caso, geométricamente se ve que hay correlación). De ahí sacamos ecuaciones que nos permiten calcular todos los elementos. Más concretamente, si giramos la elipse un ángulo $\alpha$ indeterminado, la fórmula del giro es

\[ A = \begin{pmatrix} \sin \alpha & \cos \alpha \\ -\cos \alpha & \sin \alpha \end{pmatrix}. \]

Al aplicarlo a la matriz $\Sigma$ original, mediante la fórmula $A\cdot\Sigma\cdot A^T$ que ya conocemos, obtenemos ecuaciones para $\sigma_{max}$, $\sigma_{min}$ y $\alpha$.

Hay varias fórmulas para estos cálculos, todas derivadas de la técnica descrita. Lo más sencillo es:

  • Al diagonalizar la matriz, se obtienen los cuadrados de los semiejes. Por tanto, los semiejes $\sigma_{max}$ y $\sigma_{min}$ son las raíces cuadradas de los autovalores de la matriz.
  • El ángulo $\alpha$ queda determinado por la fórmula $\tan\alpha = \displaystyle\frac{\sigma_{xy}}{\sigma_{max}^2-\sigma_x^2}$.

(La manera más habitual en la literatura parece ser calcular $\alpha$ mediante la fórmula $\tan 2\alpha =\displaystyle\frac{2\sigma_{xy}}{\sigma_x^2 - \sigma_y^2}$, lo que involucra estudiar en qué cuadrante está el ángulo, ver por ejemplo esta página; y usar fórmulas explícitas para los semiejes al cuadrado, que en realidad no son más que las soluciones de la ecuación de segundo grado que define los autovalores).

http://books.google.es/books?id=2gB7w9XlNJAC&pg=PA412&lpg=PA412&dq=error+ellipse+standard+probability&source=bl&ots=7NvzL2levX&sig=em2yul1XnPw2dSyzAJRIjcbMXgM&hl=en&sa=X&ei=lJhmUeWgM8GrhQeRh4DQBA&redir_esc=y#v=onepage&q=error%20ellipse%20standard%20probability&f=false

El último ingrediente es la probabilidad con la que queremos trabajar. Así como para intervalos de confianza usamos como coeficiente el un cuantil de la $N(0,1)$ o bien de la $t_{n-1}$, para la elipse de error el coeficiente es la raíz cuadrada de un cuantil de $\chi_2^2$. Los valores más habituales son:

 

c P
1 39.4 %
1.177 50 %
2.447 95 %
3.035 99 %

En resumen, para calcular la elipse de error:

  • Calculamos los semiejes mayor y menor y el ángulo de giro
  • Calculamos el coeficiente a partir del cuantil deseado
  • Dibujamos la elipse, centrada en el punto de interés y con los semiejes multiplicados por el coeficiente

El código siguiente que define una función con la que podrás dibujar las elipses que quieras. Toma tres argumentos: el centro, la matriz 2x2 y la probabilidad. Opcionalmente se puede cambiar el color y añadir un cambio de escala (por ejemplo una elipse de pocos centímetros en un mapa de kilómetros no se vería así que la querríamos hacer 10000 veces más grande).

def elipse(centro,matriz,prob,escala=1,color="blue"): c = sqrt(r.qchisq(prob,2))._sage_() e = map(sqrt,matriz.eigenvalues()) e = map(abs,e) smax = max(e).n() smin = min(e).n() try: alfa = arctan(matriz[0,1]/(smax^2-matriz[0,0])) except ZeroDivisionError: alfa = 0 return ellipse(centro, escala*c*smax, escala*c*smin, pi/2-alfa, color=color) elipse 
       
<function elipse at 0x7f671564b668>
<function elipse at 0x7f671564b668>
 
       

Ejemplo (continuación)

El centro y la matriz.

F(vd,va); Sxy 
       
(202.254248593737, 146.946313073118)
[ 0.000198648680058493 -0.000135778259710855]
[-0.000135778259710855  0.000286882741859060]
(202.254248593737, 146.946313073118)
[ 0.000198648680058493 -0.000135778259710855]
[-0.000135778259710855  0.000286882741859060]

Dibujamos dos elipses con esos datos, junto con la nube de puntos de antes.

punto = F(vd,va) gelipse1 = elipse(punto, Sxy, 0.95) gelipse2 = elipse(punto, Sxy, 0.999, color="red") (nube+gelipse1+gelipse2).show(aspect_ratio=1, figsize=4) 
       

Si dibujamos la elipse en relación al origen no vemos nada porque es muy pequeña.

(nube+gelipse1+point((0,0))).show(aspect_ratio=1, figsize=4) 
       

Así que vamos a dibujarla con un factor de escala.

gelipse1grande = elipse(punto, Sxy, 0.95,escala=100) (nube+gelipse1grande+point((0,0))).show(aspect_ratio=1, figsize=4)