AJO_1415_Ejemplo_ajusteparamétrico_lineal

959 days ago by sevillad1415

Un ejemplo sencillo de ajuste lineal para ver todos los conceptos.

Queremos saber los tres ángulos de un triángulo plano. Mediremos los tres ángulos y ajustaremos. Antes de eso haremos un estudio de la precisión del ajuste usando datos sobre la precisión de las mediciones.

Las observaviones serán $l_1$, $l_2$ y $l_3$. El número de parámetros es dos, llamados $x,y$, correspondiendo a los valores reales de los dos primeros ángulos. Las ecuaciones serán:

\[ l_1+v_1 = x; \quad l_2+v_2=y; \quad l_3+v_3=200-(x+y)\]

Empezamos por escribir el planteamiento, es decir, las matrices $A$ y $L$. Todavía no tenemos valores para $L$ así que la dejamos indicada.

var('l1,l2,l3') L = vector([l1,l2,l3-200]).column() A = matrix([[1,0],[0,1],[-1,-1]]) A; L 
       


Precisión del ajuste

Para la precisión necesitamos $A$ (todavía no tenemos observaciones), para el ajuste necesitaremos $L$ también, pero eso viene luego.

La precisión puede venir dada explícitamente (desviaciones típicas, varianzas, correlaciones...) o implícitamente (matriz de pesos). En el primer caso conoceremos la varianza de referencia (de hecho la elegiremos nosotros), en el segundo caso no la conoceremos.

Resolveremos un ejemplo de cada caso.

Desde la matriz de pesos

Supongamos que solo tenemos los pesos de las observaciones. Supongamos que la primera es más precisa que las otras dos.

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

La matriz cofactor de las observaciones, llamada $Q_L$, es la inversa de $P$. Las varianzas de las observaciones serán las entradas de esa matriz multiplicadas por la varianza de referencia $\sigma^2$ (desconocida). Aunque no sepamos cuanto valen las varianzas exactamente, sí sabemos lo siguiente:

  • Si el peso es 1, la varianza es igual a la varianza de referencia.
  • Podemos saber las relaciones entre las varianzas. Por ejemplo la primera varianza es igual a la mitad de la segunda.
QL = P^(-1); QL 
       

Propagamos a los parámetros para saber sus precisiones (salvo el factor de escala desconocido, como siempre), y después propagamos a las observaciones corregidas. En algunos problemas nos interesará lo primero, en otros lo segundo. En este caso nos interesan los tres ángulos ($L+V$), no solo dos ($X$); en otros problemas nos interesarán otras cosas.

N = transpose(A)*P*A QX = N^(-1) QLV = A*QX*transpose(A) QLV 
       

Comparando las matrices $Q_L$ y $Q_{L+V}$ observamos:

  • Los tres ángulos tienen mejor precisión después de ajuste (sus varianzas han disminuido). En general, siempre mejora la precisión (una medida aislada tiene menos información que al juntarla con otras).
  • El primer ángulo pasó de tener varianza $\frac12\sigma^2$ a tener $\frac25\sigma^2$: una pequeña mejora. Los otros dos ángulos bajaron más en comparación.
  • Las observaciones ajustadas tienen correlaciones (lógico si pensamos que todas se usan juntas para compensar, es decir, se afectan unas a otras).

No podemos decir mucho más.

Desde la matriz de varianzas y covarianzas

Probemos ahora de otra manera. Supongamos que conocemos las desviaciones típicas de las observaciones: $0.01$, $0.05$ y $0.02$. Fíjate en que al definir la matriz ya la aproximamos a pocos decimales de precisión.

var('l1,l2,l3') L = vector([l1,l2,l3-200]).column() A = matrix([[1,0],[0,1],[-1,-1]]) SigmaL = diagonal_matrix([.01^2, .05^2, .02^2]).n(digits=4) SigmaL 
       

Ahora elegimos la varianza de referencia. Elijamos la de la segunda observación.

vdr = SigmaL[1,1] QL = 1/vdr * SigmaL vdr; QL 
       


La propagación se hace como antes.

P = QL^(-1) N = transpose(A)*P*A QX = N^(-1) QLV = A*QX*transpose(A) QLV 
       

Y como antes podemos comparar cuánto mejora la precisión de cada observación al ajustar.

Pero ahora podemos llegar más lejos. Como conocemos la varianza de referencia, podemos usarla para pasar de la matriz cofactor ($Q$) a la matriz de varianzas y covarianzas ($\Sigma$). Con esto tendremos los valores de las desviaciones típicas, por ejemplo. Las calculamos a continuación.

SigmaLV = QLV*vdr SigmaLV; [sqrt(SigmaLV[i][i]) for i in [0,1,2]] 
       


Para calcular intervalos de confianza necesitamos los valores ajustados y sus precisiones. Acabamos de calcular lo segundo. Aunque todavía no tengamos los valores ajustados, ya sabemos la anchura de los intervalos. Por ejemplo, calculemos las anchuras de los intervalos al 95% de los ángulos ajustados.

[4*sqrt(SigmaLV[i,i]) for i in [0,1,2]] 
       

Ajuste de observaciones

Hagamos dos ajustes para este problema, con distintas informaciones sobre la precisión como más arriba. Con esto veremos la utilidad de trabajar con la varianza de referencia.

Desde la matriz de pesos

A continuación daremos unos valores a las observaciones y ajustaremos. Supongamos que las observaciones tienen los pesos de arriba.

A = matrix([[1,0],[0,1],[-1,-1]]) P = diagonal_matrix([2,1,1]) L = vector([61,72,66-200]).column() L 
       

El ajuste calcula $X$, que en este problema no es lo que necesitamos: estamos buscando $L+V$. Los residuos son el vector columna $V$, que se calcula con la fórmula habitual. Si los sumamos a $L$ estamos ajustando las observaciones, pero como había una constante matemática metida en $L$ (200 grados) la compensamos sumándola.

N = transpose(A)*P*A T = transpose(A)*P*L X = N^(-1)*T V = A*X-L LmasV = L+V+vector([0,0,200]).column() V.n(digits=4); LmasV.n(digits=4) 
       


Finalmente, estimaremos la varianza de referencia (recuerda que en este ajuste no la conocemos, solo sabemos que los pesos son iguales).

redundancia = 3-2 vdr_estim = transpose(V)*P*V/redundancia vdr_estim = vdr_estim[0,0] vdr_estim 
       

Usando este valor y haciendo el cálculo de $Q_{L+V}$, una vez más podemos calcular las desviaciones típicas finales.

QX = N^(-1) QLV = A*QX*transpose(A) SigmaLV = QLV*vdr_estim dtipicas = [sqrt(SigmaLV[i][i]).n(digits=5) for i in [0,1,2]] SigmaLV; dtipicas 
       


De aquí podemos calcular intervalos de confianza para las tres observaciones ajustadas. Hay que tener en cuenta que la varianza no es conocida sino estimada. Por eso usamos no el factor 2 (que sale de la normal estándar) sino un cuantil de una $t$ de Student, que será mayor.

c = r.qt(0.975,redundancia).n() c 
       
[(LmasV[i,0]-c*dtipicas[i],LmasV[i,0]+c*dtipicas[i]) for i in [0,1,2]] 
       

Desde la matriz de varianzas y covarianzas

Usemos los datos anteriores.

A = matrix([[1,0],[0,1],[-1,-1]]) L = vector([61,72,66-200]).column() SigmaL = diagonal_matrix([.01^2, .05^2, .02^2]).n(digits=4) L; SigmaL 
       


Elijamos la misma varianza de referencia que antes.

vdr = SigmaL[1,1] QL = 1/vdr * SigmaL P = QL^(-1) vdr; QL; P 
       




Dejemos calculado el resultado del ajuste aunque eso no sea lo que nos importa ahora.

N = transpose(A)*P*A T = transpose(A)*P*L X = N^(-1)*T V = A*X-L LmasV = L+V+vector([0,0,200]).column() V.n(digits=4); LmasV.n(digits=4) 
       


Ahora viene lo nuevo: estimemos la varianza de referencia. ¡Pero si ya la sabemos! Sí, y eso nos va a permitir comparar el valor real con el valor estimado.

redundancia = 3-2 vdr_estim = transpose(V)*P*V/redundancia vdr_estim = vdr_estim[0,0] vdr_estim.n(); vdr 
       


Observamos que la varianza de referencia estimada es mucho mayor que la real (varios órdenes de magnitud). Esto significa que hay una seria disparidad entre nuestras observaciones y nuestras precisiones: la estimación nos dice que nuestras observaciones no son tan precisas como creíamos. En efecto, fijémonos en lo siguiente:

  • Las desviaciones típicas son $0.01$, $0.05$ y $0.02$. Por tanto los intervalos de confianza al 95% de esos tres ángulos son $\pm0.02$, $\pm0.1$ y $\pm0.04$.
  • Sin embargo, el error acumulado en los tres ángulos es $61+72+66-200=1$ grado.

Es extremadamente improbable que con esas precisiones se haya acumulado tanto error. La conclusión es que las observaciones no son coherentes con las precisiones.