AJO_1415_Ejemplo_ajusteparametrico_nolineal1

959 days ago by sevillad1415

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 
       
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.

aprox1 = (l_d*sin(l_alfa*pi/200).n(),l_d*cos(l_alfa*pi/200).n()) aprox2 = (l_a,sqrt(l_d^2-l_a^2).n()) aprox3 = (l_a,l_a/tan(l_alfa*pi/200).n()) print aprox1, aprox2, aprox3 point([aprox1,aprox2,aprox3],size=40,figsize=4) 
       
(196.107237801441, 120.174669884668) (196, 120.349491066643) (196,
120.108954475427)
(196.107237801441, 120.174669884668) (196, 120.349491066643) (196, 120.108954475427)

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 tomamos valores razonablemente buenos, similares a los ya estimados más arriba: elegimos un punto más o menos en el centro del triángulo que dibujamos antes.

X0 = vector([196.05,120.2]).column() X0 
       

Podemos hacernos una idea de la bondad de este punto inicial calculando los residuos, $V=F(X_0)-L$. También la suma de cuadrados $V^T \cdot P \cdot V$, 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 
       
X1 = X0 + DeltaX0 X0; X1 
       


Ahora comparamos los residuos del punto original y del punto corregido, y sus 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] V0; V1 
       


sumaV0; sumaV1 
       


Usamos el nuevo punto para crear otro sistema lineal aproximado, la segunda iteración.

J1 = jacobian(F,(x,y))(*X1.column(0)) DeltaL1 = L - F(*X1.column(0)).column() 
       

Resolvemos este ajuste lineal para calcular la nueva corrección del punto.

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 las primeras cifras decimales 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:

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.

 

Bonadad 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) a partir del punto que tenemos.

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


Para calcular las matrices de covarianzas, necesitamos la varianza de referencia (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% (luego $c=2.447$).

e = SigmaX.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaX[0,1]/(max(e)-SigmaX[0,0])) c = 2.447 p1 = ellipse((X2[0,0],X2[1,0]), c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1) + point((X2[0,0],X2[1,0])) p2 = ellipse((X2[0,0],X2[1,0]), c*smax, c*smin, pi/2-alfa, figsize=5, xmin=190,xmax=202,ymin=115,ymax=125, aspect_ratio=1) + point((X2[0,0],X2[1,0])) p1.show(); p2.show()