AJO_1617_Propagacion_ejemplo1

hace 194 días por sevillad1617

Ejemplo de propagación no lineal

Problema: dados un acimut y una distancia con sus precisiones (d.t.), calcular las coordenadas del punto resultante y sus precisiones (d.t. de cada coordenada, correlación, y dibujo de elipses de error).

 

Primero escribimos la función que convierte los datos a coordenadas. Como queremos trabajar en grados, pero las funciones trigonométricas están en radianes, dentro metemos el factor de conversión.

f(alfa,d) = ( d*sin(pi*alfa/200) , d*cos(pi*alfa/200) ) f 
       
(alfa, d) |--> (d*sin(1/200*pi*alfa), d*cos(1/200*pi*alfa))
(alfa, d) |--> (d*sin(1/200*pi*alfa), d*cos(1/200*pi*alfa))

A continuación tenemos los datos. Evaluamos la función en esos datos para asegurarnos de que está bien.

ld, sd = 50, 0.01 lalfa, salfa = 70, 0.0020 f(lalfa, ld).n() point(f(lalfa, ld),size=30,aspect_ratio=1,figsize=4) + point((0,0)) 
       

Con las precisiones construimos la matriz de covarianzas inicial. El orden de las variables debe ser el que usamos con la función.

S1 = diagonal_matrix([salfa^2, sd^2]) S1 
       
[ 4.00000000000000e-6    0.000000000000000]
[   0.000000000000000 0.000100000000000000]
[ 4.00000000000000e-6    0.000000000000000]
[   0.000000000000000 0.000100000000000000]

Calculamos el jacobiano y lo guardamos evaluado en los datos.

J = jacobian(f,(alfa,d)) J = J(lalfa,ld) J 
       
[ 1/4*pi*cos(7/20*pi)         sin(7/20*pi)]
[-1/4*pi*sin(7/20*pi)         cos(7/20*pi)]
[ 1/4*pi*cos(7/20*pi)         sin(7/20*pi)]
[-1/4*pi*sin(7/20*pi)         cos(7/20*pi)]

Ahora simplemente hacemos la propagación. La matriz de covarianzas resultante contiene la información estadística de las coordenadas, en el orden en que escribimos la función.

S2 = (J*S1*transpose(J)).n() S2 
       
[0.0000798978121756447 0.0000394527650077175]
[0.0000394527650077175 0.0000225695889246277]
[0.0000798978121756447 0.0000394527650077175]
[0.0000394527650077175 0.0000225695889246277]

A continuación calculamos las desviaciones típicas y la correlación entre las coordenadas.

sx, sy = sqrt(S2[0,0]), sqrt(S2[1,1]) rxy = S2[0,1]/sx/sy sx; sy; rxy 
       
0.00893855761158615
0.00475074614398914
0.929069445421225
0.00893855761158615
0.00475074614398914
0.929069445421225

Pasamos a la elipse de error. En otra hoja tienes el resumen teórico con algún ejemplo de cálculo.

Los semiejes mayor y menor son las raíces cuadradas de los autovalores de la matriz (aunque no sepas nada de autovalores, basta con usar el código siguiente para calcularlos). El semieje mayor es un indicador importante de la precisión.

e = map(sqrt,S2.eigenvalues()) smax = max(e).n() smin = min(e).n() smax; smin 
       
0.0100000000000000
0.00157079632679490
0.0100000000000000
0.00157079632679490

A continuación se 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.

def elipse(centro,matriz,prob,color="blue"): c = sqrt(r.qchisq(prob,2))._sage_() e = map(sqrt,matriz.eigenvalues()) 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, c*smax, c*smin, pi/2-alfa, color=color) elipse 
       
<function elipse at 0x7fd9189a7050>
<function elipse at 0x7fd9189a7050>

La usamos con nuestro ejemplo para dibujar al 95% y al 99.9%. El centro es la estimación puntual, que se obtiene directamente evaluando la función como hicimos arriba. Es mejor guardar el gráfico en una variable y dibujarlo después con opciones adicionales (tamaño, relación entre los ejes, etc.)

punto = f(lalfa, ld) gelipse1 = elipse(punto, S2, 0.95) gelipse2 = elipse(punto, S2, 0.999, color="red") (gelipse1+gelipse2).show(aspect_ratio=1, figsize=4) 
       

Dibujamos la elipse en referencia al origen de coordenadas.

gorigen = point([0,0]) (gelipse1+gorigen).show(aspect_ratio=1, figsize=7) 
       

Ahora cambia los datos iniciales y ejecuta el código de nuevo para ver cómo cambia la forma de la elipse. Puedes ejecutarlo todo de una vez en el menú "Acción" -> "Evaluar todo".