AJO_1617_Ajuste_Ejemplonolineal1

hace 158 días por etlopez

Ejemplo de ajuste no lineal con el método de relaciones de observación

Calculamos el ajuste de un problema no lineal sencillo con el método visto en clase: iteración de ajustes lineales.

 

Queremos conocer las coordenadas de un punto en el plano. Hemos conseguido medir su distancia al origen, acimut (también con el origen) y coordenada horizontal.

 

  • $d=230\ m$ con peso 10
  • $\alpha=65^g$ con peso 10
  • $a=196\ m$ con peso 1

El vector columna $L$ contiene las tres observaciones. Elegimos trabajar con ángulos en grados en lugar de radianes; tendremos que prestar atención a la hora de usar funciones trigonométricas, que siempre son en radianes.

l_d = 230 l_alfa = 65 l_a = 196 L = vector([l_d,l_alfa,l_a]).column() L 
       

                                
                            

                                

Usamos los pesos dados. En la práctica no se tiene directamente los pesos; ¿qué sentido tendría decir que una medida de un ángulo es diez veces más precisa que otra de una distancia? Pero para ilustrar después la bondad del ajuste utilizaremos estos datos.

P = diagonal_matrix([10,10,1]) P 
       

                                
                            

                                

A partir de dos de los tres datos podemos calcular las coordenadas del punto, de tres maneras distintas (ver a continuación). Por tanto el número de parámetros es dos, y la redundancia es uno.

Por cierto, existen métodos geométricos para resolver este tipo de problemas: teniendo tres puntos aproximados, tomamos uno "en medio". Pero los métodos de ajuste estadísticos son mejores porque podemos estimar la bondad del ajuste, las precisiones, etc.

Elegimos como parámetros las coordenadas del punto $(x,y)$, pues son lo que nos interesa al fin y al cabo. Podríamos haber elegido dos magnitudes cualesquiera, sólo con la siguiente condición: que las tres observaciones se puedan poner como funciones de los parámetros.

A continuación escribimos las tres observaciones en función de las coordenadas del punto (los dos parámetros elegidos). Usaremos esta función para varias cosas durante el cálculo. Tomamos las observaciones en el orden dado. Como queremos mantener las unidades del ángulo, y la función atan2 nos devuelve radianes, mutiplicamos por la constante de conversión de radianes a grados centesimales.

var('x,y') F(x,y) = (sqrt(x^2+y^2), atan2(x,y)*200/pi, x) 
       

Recordemos que necesitamos linealizar el problema. Esto significa que:

  • Empezaremos eligiendo una aproximación para los parámetros $X_0$.
  • No calcularemos los ajustes de las $X$, sino desviaciones $\Delta X_0$.
  • Para ello no usaremos las observaciones $L$ sino desviaciones $\Delta L_0$.
  • Linealizaremos la función $F$ tomando su jacobiano.

Como punto inicial podemos tomar por ejemplo las estimaciones de  x mediante d*sin(alpha), e y mediante d*cos(alpha)

X0 = vector([l_d*sin(l_alfa*pi/200).n(),l_d*cos(l_alfa*pi/200).n()]).column() X0 
       

                                
                            

                                

Podemos hacernos una idea de la bondad de este punto inicial calculando los residuos, $V_0=F(X_0)-L$. También la suma de cuadrados $V_0^t \cdot P \cdot V_0$, aunque no tenemos manera de decir si el valor estará cerca del mínimo.

V0 = L - F(*X0.column(0)).column().n() sumaV0 = transpose(V0)*P*V0 sumaV0 = sumaV0[0,0] V0; sumaV0 
       

                                
                            

                                

El punto elegido, $X_0$, se usa para escribir las matrices que componen el ajuste lineal.

  • $A$ se sustituye por el jacobiano de $F$ en $X_0$.
  • $\Delta L_0$ es la diferencia entre los $L$ observados y los valores de $F$ en $X_0$.
J0 = jacobian(F,(x,y))(*X0.column(0)) DeltaL0 = L - F(*X0.column(0)).column() 
       

Resolvemos la primera iteración, un ajuste lineal normal y corriente. Lo que obtendremos será $\Delta X_0$, es decir, cuánto tenemos que modificar $X_0$ para acercarnos a la solución. Tras calcular $\Delta X_0$ ya tenemos un punto mejor, $X_1$.

DeltaX0 = (transpose(J0)*P*J0)^(-1) * (transpose(J0)*P*DeltaL0) DeltaX0 = DeltaX0.n() DeltaX0 
       

                                
                            

                                

Observamos que hay una corrección a partir de las centésimas. Veamos el punto inicial y el nuevo.

X1 = X0 + DeltaX0 X0; X1 
       

                                
                            

                                

Ahora comparamos las sumas ponderadas de cuadrados; ha habido una clara disminución de esa cantidad al corregir el punto.

V1 = L - F(*X1.column(0)).column().n() sumaV1 = transpose(V1)*P*V1 sumaV1 = sumaV1[0,0] sumaV0; sumaV1 
       

                                
                            

                                

Usamos el nuevo punto para crear otro sistema lineal aproximado, la segunda iteración. Resolvemos este ajuste lineal para calcular la nueva corrección del punto.

J1 = jacobian(F,(x,y))(*X1.column(0)) DeltaL1 = L - F(*X1.column(0)).column() DeltaX1 = (transpose(J1)*P*J1)^(-1) * (transpose(J1)*P*DeltaL1) DeltaX1 = DeltaX1.n() DeltaX1 
       

                                
                            

                                

La nueva corrección es mucho más pequeña que la del paso anterior, comparémoslas.

DeltaX0; DeltaX1 
       

                                
                            

                                

Al añadir la corrección calculada en la segunda iteración, vemos que varias cifras decimales ya no cambian; seguramente las cuatro o cinco primeras cifras decimales ya sean correctas.

X2 = X1 + DeltaX1 X1; X2 
       

                                
                            

                                

A continuación comparamos las tres sumas de residuos obtenidas hasta ahora.

V2 = L - F(*X2.column(0)).column().n() sumaV2 = transpose(V2)*P*V2 sumaV2 = sumaV2[0,0] sumaV0; sumaV1; sumaV2 
       

                                
                            

                                

Vemos que la suma de residuos se está estabilizando; el último punto es una aproximación buena y podríamos parar aquí (ahora vemos también que el punto anterior era tan bueno como el nuevo, pero eso no podíamos saberlo). Calculemos una iteración más por ver qué pasa.

J2 = jacobian(F,(x,y))(*X2.column(0)) DeltaL2 = L - F(*X2.column(0)).column() DeltaX2 = (transpose(J2)*P*J2)^(-1) * (transpose(J2)*P*DeltaL2) DeltaX2 = DeltaX2.n() DeltaX2 
       

                                
                            

                                
X3 = X2 + DeltaX2 V3 = L - F(*X3.column(0)).column().n() sumaV3 = transpose(V3)*P*V3 sumaV3 = sumaV3[0,0] sumaV2; sumaV3 
       

                                
                            

                                

El punto ya no cambia apenas. De hecho, nos ha pasado algo curioso: estamos tan cerca del mínimo buscado, que "nos hemos pasado un poco de largo" y ha aumentado ligeramente la suma de cuadrados. Nos quedaremos con el punto anterior, que es un pelín mejor, aunque la diferencia es tan pequeña que nos daría igual.

En resumen, las coordenadas del punto y las observaciones ajustadas son las siguientes.

X2; F(*X2.column(0)).column().n() 
       

                                
                            

                                

Para terminar con el ajuste en sí, supongamos que no hubiéramos elegido un punto inicial relativamente bueno, como al principio, sino que lo elegimos realmente mal.

Xmalo = vector([400,30]).column() Xmalo; F(*Xmalo.column(0)).column().n() 
       

                                
                            

                                

Residuos y suma de cuadrados (es enorme comparada con el mínimo que ya hemos calculado antes):

Vmalo = L - F(*Xmalo.column(0)).column().n() sumaVmalo = transpose(Vmalo)*P*Vmalo sumaVmalo = sumaVmalo[0,0] Vmalo; sumaVmalo 
       

                                
                            

                                

Una iteración:

Jmalo = jacobian(F,(x,y))(*Xmalo.column(0)) DeltaLmalo = L - F(*Xmalo.column(0)).column() DeltaXmalo = (transpose(Jmalo)*P*Jmalo)^(-1) * (transpose(Jmalo)*P*DeltaLmalo) DeltaXmalo = DeltaXmalo.n() DeltaXmalo 
       

                                
                            

                                

Punto corregido y su suma de cuadrados:

Xmenosmalo = Xmalo + DeltaXmalo Vmenosmalo = L - F(*Xmenosmalo.column(0)).column().n() sumaVmenosmalo = transpose(Vmenosmalo)*P*Vmenosmalo sumaVmenosmalo = sumaVmenosmalo[0,0] Xmenosmalo; sumaVmenosmalo 
       

                                
                            

                                

Vemos que el primer parámetro (la coordenada horizontal) se ha acercado bastante al valor óptimo, pero el segundo parámetro (la coordenada vertical) se ha pasado mucho de largo. Eso se debe a la naturaleza del problema que tenemos: recordemos que una de las observaciones es la coordenada horizontal precisamente, por eso tendemos a su valor óptimo más rápido que la otra coordenada. Una iteración más a partir del punto corregido nos acercaría bastante al valor buscado, y una o dos más nos dejaría prácticamente en el mismo punto que obtuvimos cuando lo hicimos bien.

 

Bondad y precisión del ajuste

 

Como no conocemos las precisiones de las observaciones, no podemos estimar la bondad del ajuste (recuerda que se comparaba la varianza de referencia, que la elegiríamos, con la varianza de referencia estimada, también llamada a posteriori).

 

Veamos el cálculo de la precisión. Es tan sencillo como usar las fórmulas que ya conocemos, con las matrices aproximadas que ya tenemos de antes, y con la varianza de referencia a posteriori correspondiente. Usamos la mejor aproximación que hemos encontrado.

X2 
       

                                
                            

                                

Las matrices cofactor de las X y de las observaciones se calculan como siempre:

  • $Q_X = N^{-1} = (A^tPA)^{-1}$
  • $Q_{L+V} = AQ_XA^t = A N^{-1} A^t$

pero teniendo en cuenta que $A$ es ahora la aproximación (jacobiano) usando el punto tomado.

Recuerda que las matrices $Q$ no contienen las varianzas y covarianzas. Con lo que sabemos nos valdrían, como mucho, para comparar con las observaciones sin ajustar. Por ejemplo comparemos $Q_{L+V}$ con $Q_L$ (se confirma que mejoran las precisiones después de ajustar, y que las medidas más imprecisas por regla general son las que más mejoran):

QX = (transpose(J2)*P*J2)^(-1) QX = QX.n() QLV = J2*QX*transpose(J2) QLV = QLV.n() P^(-1); QLV 
       

                                
                            

                                

Para calcular las matrices de covarianzas ($\Sigma$) necesitamos la varianza de referencia, y como no la tenemos calculamos su estimación. Como dijimos antes, la redundancia es $3-2=1$.

redund = 3-2 vdr = sumaV2/redund vdr 
       

                                
                            

                                
SigmaX = vdr*QX SigmaLV = vdr*QLV SigmaX; SigmaLV 
       

                                
                            

                                

Ahora las desviaciones típicas salen de aplicar la raíz cuadrada a los elementos de la diagonal. Hagámoslo primero para las tres observaciones ajustadas:

map(sqrt,SigmaLV.diagonal()) 
       

                                
                            

                                

Y para las dos coordenadas:

map(sqrt,SigmaX.diagonal()) 
       

                                
                            

                                

Como además hemos ajustado las coordenadas de un punto, a partir de la matriz de covarianzas calculamos algo más interesante: la elipse de error. Tomamos probabilidad 95%.

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) punto2 = (X2[0,0],X2[1,0]) point(punto2,figsize=6)+elipse(punto2,SigmaX,0.95)