AJO_Ajustenolineal_ejemplo2

77 days ago by sevillad1819

Ejemplos de ajuste no lineal

En esta ejemplo se quiere calcular las coordenadas de un punto en distintas situaciones. En cada caso haremos una estimación previa de la precisión, y luego el ajuste en sí; solo en el primer caso está hecho el ajuste, ilustrando el código a utilizar. Queda como ejercicio hacer los ajustes de los otros casos.

 

E1: a partir de la distancia y dos ángulos coordenados.

Conocemos la distancia al origen y los ángulos que se forman con los ejes.

ld = 543 sd = 1 la = 71 sa = 0.5 lb = 30 sb = 0.3 gradrad = pi.n()/200 
       

Con tres datos tenemos redundancia: el punto habría quedado determinado con la distancia y uno de los dos ángulos. De ahí nos salen dos estimaciones para el punto.

p_est = [(ld*sin(la*gradrad),ld*cos(la*gradrad)), (ld*cos(lb*gradrad),ld*sin(lb*gradrad))] p_est; point(p_est,aspect_ratio=1,figsize=4)+point((0,0),size=0) 
       

                                
                            

                                

Planteamiento: los valores observados (si fueran exactos) en función de los parámetros, para los que elegimos obviamente las coordenadas buscadas. Para comprobar que está bien, tomamos un valor de X aproximado y vemos si lo que sale al aplicar F se parece a las observaciones.

var('x,y') F(x,y) = (sqrt(x^2+y^2), atan2(x,y)/gradrad, atan2(y,x)/gradrad) J = jacobian(F,(x,y)) F(488,239).n() 
       

                                
                            

                                

Calculamos la precisión del ajuste. Nótese la elección de un punto para calcular la matriz $A$. Por ejemplo, tomemos uno de los dos puntos estimados al principio.

  • Necesitamos información sobre el valor aproximado del punto para saber la precisión del ajuste; en el caso lineal no hacía falta saber nada de los valores, solo las precisiones de las observaciones.
  • Si cambiamos el punto elegido cambian los valores de la precisión; lo que hacemos a continuación solo es una estimación (una aproximación lineal de algo que no lo es).

Experimento: prueba a cambiar arriba la precisión de la distancia y observa el efecto que tienen los cambios en la precisión del punto ajustado.

Aprevia = J(ld*sin(la*gradrad), ld*cos(la*gradrad)) SL = diagonal_matrix([sd^2,sa^2,sb^2]) vdr = 1 QL = 1/vdr*SL P = QL.inverse() Nprevia = transpose(Aprevia)*P*Aprevia QXprevia = Nprevia.inverse().n() SXprevia = vdr*QXprevia # vemos las matrices cofactor inicial y final, y las desviaciones y correlación finales SXprevia 
       

                                
                            

                                

Ahora ajustamos iterativamente. Haremos unas cuantas iteraciones y veremos cómo va mejorando el resultado. Primero creamos unos conjuntos vacíos a los que iremos añadiendo los datos de cada iteración:

  • X son las coordenadas del punto (van mejorando en cada iteración)..
  • DeltaL y A son los "ingredientes" que se usan en cada iteración (recuerda que en cada paso se hace un ajuste lineal).
  • DeltaX es el resultado de una iteración (hay que sumarlo a X para ir mejorando el punto).
  • V y sumaV son respectivamente el vector de residuos y la suma ponderada de sus cuadrados.

Cuando vayamos iterando veremos que DeltaX va disminuyendo (nos vamos acercando al valor óptimo) y sumaV también (vamos acercándonos al mínimo).

X = {} DeltaL = {} DeltaX = {} A = {} V = {} sumaV = {} 
       

La variable n contará las iteraciones. Al principio necesitamos una estimación del punto, que puede ser la que queramos, pero cuanto más cerca del valor real empecemos menos tardaremos en "llegar". Por ejemplo tomamos uno de los dos puntos estimados arriba.

Experimento: prueba a empezar con un valor de X muy alejado y observa cómo progresan las iteraciones. Para ello tendrás que ejecutar primero la celda anterior, que borra lo que haya en los conjuntos de datos.

L = vector([ld,la,lb]).column() n = 0 X[0] = vector([ld*sin(la*gradrad),ld*cos(la*gradrad)]).column() #X[0] = vector([400,300]).column() V[n] = F(*X[n].column(0)).column() - L V[n] = V[n].n() sumaV[n] = transpose(V[n])*P*V[n] sumaV[n] = sumaV[n][0,0].n() n; X[n]; V[n]; sumaV[n] 
       

                                
                            

                                

El siguiente código calcula iteraciones, basta ejecutarlo tantas veces como se quiera. Hace lo siguiente:

  • Calcula A y DeltaL, necesarios para un paso de ajuste lineal.
  • Calcula DeltaX, el resultado de ese ajuste lineal.
  • Suma ese resultado a X para mejorar el punto.
  • Calcula los nuevos residuos y su suma de cuadrados ponderada.

Tras ejecutar una iteración se muestran n, la corrección calculada, el nuevo punto y la nueva suma de cuadrados.

Podemos parar después de un número predefinido de iteraciones, o cuando DeltaX sea tan pequeño que conozcamos el punto con la precisión que buscamos, o cuando sumaV apenas deje de decrecer...

A[n] = J(*X[n].column(0)) DeltaL[n] = L - F(*X[n].column(0)).column() DeltaX[n] = (transpose(A[n])*P*A[n])^(-1) * (transpose(A[n])*P*DeltaL[n]) DeltaX[n] = DeltaX[n].n() X[n+1] = X[n] + DeltaX[n] n = n+1 V[n] = F(*X[n].column(0)).column() - L sumaV[n] = transpose(V[n])*P*V[n] sumaV[n] = sumaV[n][0,0].n() n; DeltaX[n-1]; X[n]; sumaV[n] 
       

                                
                            

                                

Para comparar, vemos la corrección y suma de cuadrados del paso anterior.

n-1; X[n-1]; DeltaX[n-2]; sumaV[n-1] 
       

                                
                            

                                

O podemos ver todos los resultados anteriores, porque con el código usado se guardan.

X; sumaV 
       

                                
                            

                                

Teniendo una buena aproximación del punto podemos repetir ahora el cálculo de la precisión y obtener la matriz de covarianzas de las dos coordenadas.

A[n] = J(*X[n].column(0)) N = transpose(A[n])*P*A[n] QX = N.inverse().n() SX = vdr*QX SX 
       

                                
                            

                                
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 
       

                                
                            

                                
elipse(X[n].column(0),SX,0.95)+point((0,0)) 
       
 
       
 
       

E2: a partir de las coordenadas mismas y la distancia

lx = 234 sx = 1 ly = 123 sy = 2 ld = 264 sd = 1 
       
# estimaciones con pares de datos l3 = [ (lx,ly), (lx,(ld^2-lx^2)^.5), ((ld^2-ly^2)^.5,ly) ] l3; point(l3,aspect_ratio=1,figsize=4); point(l3,size=1,figsize=5)+point((0,0),size=0) 
       

                                
                            

                                

# planteamiento var('x,y') F(x,y) = (x, y, sqrt(x^2+y^2)) J = jacobian(F,(x,y)) # notese la eleccion del punto para calcular A, esto afectara a todos los calculos A = J(lx,ly) 
       
# precisión del ajuste # probar a mejorar o empeorar mucho la precisión de la distancia para ver el efecto en las precisiones de las coordenadas SL = diagonal_matrix([sx^2,sy^2,sd^2]) vdr = 1 QL = 1/vdr*SL P = QL.inverse() N = transpose(A)*P*A QX = N.inverse().n() QL; QX; sx,sy; sqrt(QX[0,0]), sqrt(QX[1,1]), QX[1,0]/sqrt(QX[0,0])/sqrt(QX[1,1]) 
       






# ahora ajustamos, iterativamente # ... # ... # ... 
       

E3: a partir de los cinco datos

# simulamos observaciones a partir de un punto exacto arbitrario gradrad = pi.n()/200 sd = 15 sx = 2 sy = 4 sa = 10 sb = 5 xx = 121. yy = 635. lx = gauss(xx,sx) ly = gauss(yy,sy) ld = gauss(sqrt(xx^2+yy^2),sd) la = gauss(atan2(yy,xx)/gradrad,sa) lb = gauss(atan2(xx,yy)/gradrad,sb) 
       
# estimaciones con pares de datos: ahora tendríamos 10 pares, no lo hacemos 
       
# planteamiento var('x,y') F(x,y) = (x, y, sqrt(x^2+y^2), atan2(y,x)/gradrad, atan2(x,y)/gradrad) J = jacobian(F,(x,y)) # notese la eleccion del punto para calcular A, esto afectara a todos los calculos A = J(lx,ly) 
       
# precisión del ajuste # probar a mejorar o empeorar mucho la precisión de la distancia para ver el efecto en las precisiones de las coordenadas SL = diagonal_matrix([sx^2,sy^2,sd^2,sa^2,sb^2]) vdr = 1 QL = 1/vdr*SL P = QL.inverse() N = transpose(A)*P*A QX = N.inverse().n() QL; QX; sx,sy; sqrt(QX[0,0]), sqrt(QX[1,1]), QX[1,0]/sqrt(QX[0,0])/sqrt(QX[1,1]) 
       






# ahora ajustamos, iterativamente # ... # ... # ...