AJO_1415_Resumen_elipses_error

932 days ago by sevillad1415

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.

Ejemplo: Figura 6 en http://www.fig.net/commission7/reports/events/penang97/penang974.htm

Más concretamente, supongamos que tenemos un punto de coordenadas reales $(\mu_x,\mu_y)$ y conocemos la imprecisión de las coordenadas (error máximo, etc.). En nuestro caso trabajaremos con 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, donde $c$ es una constante que depende del valor de certidumbre elegido (por ejemplo $c=1$ corresponde a un 68% de las medidas cayendo en el intervalo, $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 un momento):

 

Pero esta representación tiene dos problemas:

  • Da una falsa impresión de uniformidad: las zonas de distintos colores tienen probabilidades muy distintas, pues para la zona roja 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$.

# tomamos el ejemplo mencionado arriba: coordenadas polares a rectangulares vr = 20; sr = 1 va = 40; sa = 1 va = va*pi.n()/200; sa = sa*pi.n()/200 # conversión a radianes Sra = matrix([[sr^2,0],[0,sa^2]]) Sra 
       
[    1.00000000000000    0.000000000000000]
[   0.000000000000000 0.000246740110027234]
[    1.00000000000000    0.000000000000000]
[   0.000000000000000 0.000246740110027234]
# calculamos la propagación a las coordenadas x, y var('r,a') F(r,a) = (r*cos(a), r*sin(a)) J = jacobian(F,(r,a))(vr,va) Sxy = J*Sra*transpose(J) Sxy 
       
[0.688607141754449 0.428595500253020]
[0.428595500253020 0.410088902256445]
[0.688607141754449 0.428595500253020]
[0.428595500253020 0.410088902256445]
# finalmente, los valores y sus precisiones sx = Sxy[0,0]^0.5 sy = Sxy[1,1]^0.5 rxy = Sxy[0,1]/sx/sy F(vr,va); sx; sy; rxy 
       
(16.1803398874989, 11.7557050458495)
0.829823560616622
0.640381840979618
0.806534292365161
(16.1803398874989, 11.7557050458495)
0.829823560616622
0.640381840979618
0.806534292365161
 
       

Ahora aprenderemos a calcular la elipse de error a partir de los datos anteriores. Supondremos que $\mu_x=\mu_y=0$, en otro caso sólo hay que hacer una traslación. Por tanto nuestras elipses estarán centradas en el origen.

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
# ejemplo: coordenadas polares a rectangulares var('r,a') F(r,a) = (r*cos(a), r*sin(a)) vr = 20; sr = 0.01 va = 70; sa = 0.05 va = va*pi.n()/200; sa = sa*pi.n()/200 # conversión a radianes Sra = matrix([[sr^2,0],[0,sa^2]]) F(vr,va); Sra # estimación puntual y matriz de varianzas de las mediciones 
       
(9.07980999479094, 17.8201304837674)
[0.000100000000000000    0.000000000000000]
[   0.000000000000000  6.16850275068085e-7]
(9.07980999479094, 17.8201304837674)
[0.000100000000000000    0.000000000000000]
[   0.000000000000000  6.16850275068085e-7]
# calculamos la propagación a las coordenadas x, y J = jacobian(F,(r,a))(vr,va) Sxy = J*Sra*transpose(J) sx = Sxy[0,0]^0.5 sy = Sxy[1,1]^0.5 rxy = Sxy[0,1]/sx/sy Sxy; sx; sy; rxy 
       
[  0.000216495891310509 -0.0000593576213842410]
[-0.0000593576213842410   0.000130244218716725]
0.0147137993499473
0.0114124589250838
-0.353486174582213
[  0.000216495891310509 -0.0000593576213842410]
[-0.0000593576213842410   0.000130244218716725]
0.0147137993499473
0.0114124589250838
-0.353486174582213
e = [sqrt(w).n() for w in Sxy.eigenvalues()] smax = max(e) smin = min(e) alfa = arctan(Sxy[0,1]/(smax^2-Sxy[0,0])) e; smax; smin; alfa*200/pi.n() 
       
[0.00999999999999999, 0.0157079632679490]
0.0157079632679490
0.00999999999999999
-70.0000000000001
[0.00999999999999999, 0.0157079632679490]
0.0157079632679490
0.00999999999999999
-70.0000000000001

Finalmente, el coeficiente $c$ que multiplica a ambos semiejes para cambiar la probabilidad sigue básicamente una distribución chi cuadrado. Los valores más habituales son:

c P
1 39.4 %
1.177 50 %
2.447 95 %
3.035 99 %
# en el ejemplo anterior, dibujamos la elipse estándar y la elipse de confianza al 95% c = 2.447 plot(ellipse(F(vr,va), smax, smin, pi/2-alfa),figsize=6,aspect_ratio=1) + plot(ellipse(F(vr,va), c*smax, c*smin, pi/2-alfa),figsize=6,aspect_ratio=1) + point(F(vr,va)) 
       

El siguiente es el ejemplo lineal de los apuntes de propagación (media y diferencia de dos variables).

M = matrix([[0.00065,-0.0012],[-0.0012,0.0026]]) sqrt(M[0,0]); sqrt(M[1,1]); M[0,1]/sqrt(M[0,0])/sqrt(M[1,1]) 
       
0.0254950975679639
0.0509901951359278
-0.923076923076923
0.0254950975679639
0.0509901951359278
-0.923076923076923
e = [sqrt(w).n() for w in M.eigenvalues()] smax = max(e) smin = min(e) alfa = arctan(M[0,1]/(smax^2-M[0,0])) e; smax; smin; alfa*200/pi.n() 
       
[0.0563130944772761, 0.00887892957475039]
0.0563130944772761
0.00887892957475039
-28.2811895076503
[0.0563130944772761, 0.00887892957475039]
0.0563130944772761
0.00887892957475039
-28.2811895076503
plot(ellipse((0,0), smax, smin, pi/2-alfa),figsize=6,aspect_ratio=1) + point((0,0)) 
       

Finalmente, el último ejemplo de los apuntes de propagación.

# ejemplo 2 (coordenadas por acimut) M = matrix([[0.00343195,0.00345251],[0.00345251,0.00433372]]) sqrt(M[0,0]); sqrt(M[1,1]); M[0,1]/sqrt(M[0,0])/sqrt(M[1,1]) 
       
0.0585828473189892
0.0658309957390893
0.895228799552997
0.0585828473189892
0.0658309957390893
0.895228799552997
e = [sqrt(w).n() for w in M.eigenvalues()] smax = max(e) smin = min(e) alfa = arctan(M[0,1]/(smax^2-M[0,0])) e; smax; smin; alfa*200/pi.n() 
       
[0.0858176116979262, 0.0200251722205817]
0.0858176116979262
0.0200251722205817
45.8663851561901
[0.0858176116979262, 0.0200251722205817]
0.0858176116979262
0.0200251722205817
45.8663851561901
plot(ellipse((0,0), smax, smin, pi/2-alfa),figsize=6,aspect_ratio=1) + point((0,0))