Sobre los ajustes por mínimos cuadrados en metrología

Ángel Mª Sánchez Pérez
Laboratorio de Metrología y Metrotecnia
ETSII
Universidad Politécnica de Madrid

Jesús de Vicente y Oliva
Laboratorio de Metrología y Metrotecnia
ETSII
Universidad Politécnica de Madrid

Resumen: Este artículo presenta algunas de las peculiaridades del método de mínimos cuadrados cuando se emplea en aplicaciones metrológicas. Un gran campo de utilización de los mínimos cuadrados se refiere al análisis de series de datos con objeto de poner de manifiesto tendencias o leyes de comportamiento que permitan modelar los fenómenos de forma más compacta para su utilización presente o futura.

Con frecuencia, los datos de partida que se utilizan para el ajuste en la literatura suelen considerarse como valores exactos, sin incertidumbre, de forma que la mayor o menor exactitud del resultado se debe, exclusivamente, al desajuste entre dichos datos y el modelo empleado para el ajuste.

Sin embargo, cuando se consideran los datos con verdadero sentido metrológico deben contemplarse como valores con incertidumbres, no siempre despreciables y, en muchos casos, con contribuciones comunes que determinan correlaciones entre aquellos datos.

Se incluyen tres ejemplos para ilustrar las aplicaciones del método en casos reales.

Palabras clave: mínimos cuadrados, metrología, ajuste, incertidumbre.

Abstract: This paper presents some characteristics of least squares methods (LSMs) when used in metrological applications. A large field where LSMs are used is related with the analysis of data sets in order to discover trends or behavioural laws. That allows the phenomena to be modelled in a more compact mathematical form for present or future use.

Often, in the literature, input data used for least squares fitting are regarded as exact values without uncertainty. Therefore, the uncertainty of output parameters coming from the least squares fit is due exclusively to the lackoffit between input data and values predicted by the model (once the least squares fit is performed).

However, input data in metrological application always have uncertainties, sometimes negligible, sometimes not. Even in many cases we should consider correlation between input data.

In this paper three examples are included to illustrate the applications of LSM in actual metrological applications.

Keywords: least squares, metrology, adjust, uncertainty.

1.Introducción

En muchas ocasiones, bien no es posible o bien no es interesante, medir directamente los parámetros característicos de un mensurando. En vez de ello, se miden otros parámetros de manera directa y utilizando un cierto modelo físico-matemático que relaciona los parámetros medidos con aquellos otros que se desean conocer se genera un sistema de ecuaciones a partir del cual es posible determinar estos últimos.

Sean x₁, x₂… x_m los resultados proporcionados directamente por el sistema de medida (variables de entrada) e y₁, y₂… y_n los parámetros que se desean conocer (variables de salida). La relación entre ambos conjuntos de variables se supone conocida y expresada de la siguiente forma:

El signo “≅” en vez del signo igual (=) pone de manifiesto la variabilidad de los resultados de medida.

Se trabaja con información superabundante (m > n ) por lo que el sistema (1) está sobredeterminado y no posee, en general, solución. Sin embargo, es posible encontrar una solución y₁, y₂… y_n tal que las diferencias entre los términos primero y segundo de las ecuaciones (1) sean mínimas. Para ello, el sistema de ecuaciones “aproximadas” (1) puede transformarse en un sistema “exacto” introduciendo los denominados residuos e₁, e₂… e_m:

El vector y₁, y₂… y_n se elige de tal manera que el vector de residuos resultante e₁, e₂… e_m sea “el que más se aproxime al vector nulo”. Cuando esta condición se impone minimizando la suma cuadrática de los residuos, el método se conoce como método de ajuste por mínimos cuadrados (MM. CC.). Véase, por ejemplo, [1] (apartado 1.1.4), [2] (apartado 5.3), [3] (apartado 4.3).

Además, y dado que las desviaciones que introduce el sistema de medida suelen ser muy pequeñas, el sistema de ecuaciones no lineal antes descrito es, en general, linealizable entorno al valor nominal del punto de medida. Adoptando el nominal como nuevo origen para las variables, el sistema (1) resulta:

evaluándose las derivadas parciales para el valor nominal del punto de medida.

Análogamente, el sistema (2) se linealiza mediante

Un planteamiento general para muchas aplicaciones metrológicas puede establecerse mediante las siguientes ecuaciones y variables

que cabe escribir todavía de forma más compacta mediante x ≅ f (y,q)

Figura 1. Modelo general de ajuste por MM.CC.

En la figura 1 se presenta un despliegue del sistema de ecuaciones (5) donde se trata de determinar n parámetros, y₁, y₂… y_n, (magnitudes de salida) a partir de m datos de medida, x₁, x₂… x_m, (magnitudes de entrada), conociendo, además, un conjunto de parámetros complementarios, q₁, q₂… q_p. Al minimizar la suma cuadrática de los residuos e_i, se obtienen las ecuaciones necesarias para determinar los parámetrosy₁, y₂… y_n. Veáse, por ejemplo, [5] (apartado 3) y [6].

De forma similar a la que permitió pasar de (1) y (2) a (3) y (4), la linealización del sistema (5) determina

donde

son matrices con valores conocidos que se han evaluado en el entorno del nominal en el punto de trabajo.

En el caso más general, se incluyen valores e incertidumbres en las tres clases de variables consideradas, aunque algunas veces las incertidumbres pueden ser muy pequeñas y tratarse como despreciables.

El modelo (5) es flexible y permite recoger diferentes posibilidades. Por ejemplo, las funciones f_i
pueden ser varias o reducirse sólo a una; los datos de entrada pueden corresponder a valores de diferentes magnitudes o a reiteraciones de la misma magnitud; etc.

Dependiendo de las hipótesis de trabajo adoptadas, se analizan en los siguientes apartados tres modelos de MM. CC. que se denominan:

• Mínimos cuadrados ordinarios (Ordinary Least Squares, OLS)
• Mínimos cuadrados generalizados (Generalized Least Squares, GLS)
• Mínimos cuadrados totales (Total Least Squares, TLS)

2 Mínimos cuadrados ordinarios (OLS)

Un problema de MM.CC. ordinarios verifica las siguientes hipótesis:

• Las funciones f_i ( y₁, y₂… y_n ) no contienen ningún parámetro que proceda de un resultado de medida con incertidumbre. Dicho de otro modo, las funciones f_i  (y₁, y₂… y_n ) sólo contienen las incógnitas y₁, y₂… y_n a determinar y parámetros o coeficientes conocidos sin incertidumbre.En la práctica, las incertidumbres de los parámetros presentes en f_i ( y₁, y₂… y_n ) no tienen por qué ser estrictamente nulas. Basta con que su contribución final a las incertidumbres de las variables de salida y₁, y₂… y_n sea despreciable.
• Las variables de entrada (medidas conocidas), es decir, las x₁, x₂… x_m cumplen con lo siguiente:
  • Todas ellas poseen la misma incertidumbre típica, estimada por
    u( x₁ ) = u( x₂ ) = … = u( x_m ) = u_x
  • Todas ellas son estadísticamente incorreladas. Es decir, todos los coeficientes de correlación r ( x_i,x_j ) son nulos (para i ≠ j ). O dicho de otro modo, todas las covarianzas son nulas,u( x_i,x_j ) = 0 .Es decir, la matriz de covarianzas C_x = cov(x) del vector de variables aleatorias x = [x₁, x₂… x_m]^T responde a
    
  • La estimación inicial de u_x será confirmada o modificada por el proceso de ajuste según se explica más adelante.

Por consiguiente, los sistemas de ecuaciones (6), (7) se reducen en este caso a

donde se ha simplificado la notación mediante A = J_y.

La suma cuadrática de los residuos es

La condición de mínimo sobre Q exige , ∇Q= 0 es decir,

o sea

sistema cuya solución es

A partir de la expresión anterior se obtiene la matriz de covarianza de las variables de salida:

Por otra parte, una medida de la calidad del ajuste es:

Puede demostrarse que s² es un estimador de la varianza del ajuste. Si la varianza del ajuste es mucho mayor que la estimada inicialmente para las variables de entrada, es que estas últimas han sido infravaloradas y, en consecuencia, se adopta el valor s² para u_x².

Para ilustrar este razonamiento pueden analizarse los resultados de una comparación entre laboratorios que miden un mismo patrón. Una gran variabilidad de las medidas de desviación al patrón (s² elevada) no es consistente con unos valores mucho menores de la incertidumbre estimada por los laboratorios participantes y puede adoptarse el criterio

Una justificación estadística más detallada, adoptando normalidad para las variables implicadas, se basa en que una buena estimación de u_x determina que Q / u_x² se distribuya según una χ² con m-n grados de libertad. Es posible obtener un valor crítico χ_c² por debajo del cual la probabilidad de obtener valores de χ² sea muy alta (por ejemplo del 95 %). A este valor le corresponde otro, s_c², según (15). Cuando s ≤ s_c, se acepta como un suceso dentro del 95 % de sucesos esperables en la hipótesis de que u_x no ha sido infravalorada. En la mayor parte de los casos prácticos, los resultados obtenidos son muy similares a los del criterio (16), por lo que adoptaremos el criterio (16) cuya aplicación es inmediata.

3. Mínimos cuadrados generalizados (GLS)

Un problema de MM.CC. generalizados (ver [1], apartado 4.3, [13]) verifica las siguientes hipótesis:

• Las funciones f_i ( y₁, y₂… y_n ) no contienen ningún parámetro que proceda de un resultado de medida con incertidumbre. Dicho de otro modo, las funciones f_i ( y₁, y₂… y_n ) sólo contienen las incógnitas y₁, y₂… y_n a determinar y parámetros o coeficientes conocidos sin incertidumbre.En la práctica las incertidumbres de los parámetros presentes en f_i ( y₁, y₂… y_n ) no tienen por qué ser estrictamente nulas. Basta con que su contribución final a las incertidumbres de las variables de salida y₁, y₂… y_n sea despreciable.
• Las variables del término independiente, es decir, las x₁, x₂… x_m no suelen poseer la misma incertidumbre y es posible que, también, existan correlaciones no nulas entre ellas: para alguna pareja x_i, x_j  r ( x_i, x_j ) ≠ 0 lo cual es equivalente a u ( x_i, x_j ) ≠ 0 .
• Se supone conocida la matriz de covarianzas  del vector de variables aleatorias x = [x₁, x₂… x_m]^T que se presenta multiplicada por un factor positivo, s², desconocido, que permitirá corregir en el ajuste una posible infravaloración de  .

La idea básica para resolver un problema generalizado (GLS) es reducirlo a otro del tipo ordinario (OLS) mediante un cambio de variables a partir de una matriz de ponderación W. La nueva variable, z, viene dada por

Debiendo ser

La matriz W puede determinarse por varios procedimientos. Uno de ellos es mediante la descomposición de Cholesky (W = L^-1 siendo C_x = LL^T y L triangular inferior) sobre la matriz de covarianzas de los valores medidos, C_x = cov(x) , lo que puede implementarse directamente en MATLAB u Octave mediante la instrucción

W = inv(chol(Cx,’lower’))

Multiplicando (9) por la matriz de ponderación se obtiene

El nuevo sistema de ecuaciones, con Ã = W( A · y ), se puede resolver por MM.CC. Ordinarios (OLS) pues ya se cumplen las condiciones exigidas a las varianzas y covarianzas de las variables de entrada en dicho modelo (8).

Los residuos resultan:

Aplicando el mismo criterio que en el caso anterior (16)

La covarianza de las variables de salida, C_y = cov( y ), responde a la misma expresión general (14) por lo que

resultando afectada la estimación de la matriz inicial, C_x = cov( x ), por el mismo factor, es decir, s²C_x .

4. Mínimos cuadrados totales (TLS)

Un problema de MM. CC. totales [12] verifica las siguientes hipótesis:

• En este caso, las funciones f_i ( y₁, y₂… y_n ) contienen una serie de parámetros q₁, q₂… q_p, que se conocen con incertidumbres no despreciables por lo que las incertidumbres están presentes tanto en el término x₁, x₂… x_m como en las funciones f_i ( y₁, y₂… y_n ). Se trata de analizar el modelo general de la figura 1.
• Es conocida la matriz de covarianzas C_x = cov(x), quizá en el formato  comentado en el caso GLS.
• Asimismo, otro dato es la matriz de covarianzas  del vector de parámetros q = [q₁, q₂… q_p]^T conocidos que se presenta multiplicada por el mismo factor positivo anterior, s², desconocido, que permitirá también corregir en el ajuste una posible infravaloración de .

Suponiendo que trabajamos sobre un sistema linealizado, la expresión (6) puede escribirse

Ahora, el término de la izquierda (término independiente) es el que recoge toda la incertidumbre del problema. En el término de la derecha (matriz A) ya no hay incertidumbre. Por tanto, además de haber linealizado el problema, éste se ha reconvertido en un problema de MM.CC. Generalizados (GLS). Por tanto, es de aplicación la transformación (17) a partir de la matriz

Obteniéndose formalmente expresiones similares a (18) y (19) y, finalmente (21)

Aplicando el mismo criterio que en los casos anteriores (16) y (20), resulta

5. Ejemplo de ajuste por mínimos cuadrados ordinarios (OLS)

Los ejemplos que se presentan se refieren a un patrón de ranura utilizado en metrología dimensional como patrón de calibración de la escala vertical de rugosímetros y microscopios confocales y de efecto túnel. La altura o profundidad de la ranura es del orden de unos pocos micrómetros pero su valor suele expresarse en nanómetros con incertidumbres típicas del orden de la decena de nanómetros.

Enunciado 1

Se desea analizar la variación en función del tiempo (deriva temporal) de un patrón de ranura.

En la tabla siguiente se muestra el histórico de calibraciones externas (en Institutos Nacionales de Metrología, INM) de dicho patrón:

Los tiempos en los cuales se realizaron las respectivas calibraciones se miden en años respecto de la fecha actual.

Las calibraciones se han realizado en cuatro INMs diferentes, habiéndose repetido tres calibraciones en un mismo INM.

Se pretende realizar un ajuste del tipo: 

Donde a representaría la estimación de la altura de la ranura en la fecha actual (expresada en nm) y b representaría la deriva anual (expresada en nm/año)

Despreciando la incertidumbre de calibración, se desea realizar un ajuste por MM.CC. Ordinarios (OLS) proporcionando los resultados siguientes:

• La altura a de la ranura en la fecha actual junto con su incertidumbre típica
• La deriva anual b junto con su incertidumbre típica
• El coeficiente de correlación r;(a,b)

El enunciado anterior satisface las condiciones del análisis OLS del apartado 2, si bien se han modificado algunos símbolos ( y₁ por a; y₂ por b ) y existe un parámetro, t, cuyos valores se conocen sin incertidumbre. De acuerdo con todo ello, el sistema (1) se escribe:

La primera de las expresiones (9) , x ≅ A·y, se concreta en:

La solución (13) determina los resultados: a = 2396,0 nm y b = -2,9 nm/año .

Aunque el enunciado recoge las incertidumbres típicas de los valores certificados por los INM, el enunciado indica que, en este primer enunciado, se prescinda de aquellas, por lo que la incertidumbre final procede exclusivamente del desajuste que se deduce de (13) con valor: s = 8,15 nm. Nótese, sin embargo, que el valor obtenido para s es del mismo orden que las incertidumbres u( x ) contenidas en el enunciado nº1.

La matriz de covarianzas de las variables de salida responde a (14) con u_x² = s² . El resultado es:

y con dos cifras significativas u( a ) = 6,4 nm ; u( b ) = 1,2 nm/año; r ( a,b ) = +0,85.

La figura 2 reproduce el programa Ej1 que realiza estos cálculos en MATLAB/Octave.

6. Ejemplo de ajuste por mínimos cuadrados generalizados (GLS)

Este segundo ejemplo amplía los datos del enunciado 1 con más información.

Enunciado 2

Se trata de volver a realizar un ajuste del tipo x
(t) = a + b ·t esta vez por MM.CC. Generalizados (GLS) teniendo en cuente la siguiente información adicional:

• Las incertidumbres de calibración U(x_i)
• La presencia de correlación entre correcciones de calibración:
  • Si las correcciones de calibración han sido proporcionadas por el mismo INM se supondrá que el coeficiente de correlación es +0,5
  • Si las correcciones de calibración han sido proporcionadas por distintos INMs se supondrá que el coeficiente de correlación es nulo

Con estos datos, se cumplen las condiciones del análisis GLS del apartado 3 por lo que el proceso se inicia determinando la matriz de covarianzas, C_x , del término independiente (vector de medicionesx). En el programa de la figura 3, se trabaja con la matriz  del apartado 3, que se designa Cx, para mayorarla al final si  (20). Para ello, se obtiene primero la matriz R de coeficientes de correlación, que satisface:

• Valor unidad en cada elemento de su diagonal principal
• Valores iguales a 0,5 en las posiciones (1,3) (1,4) (3,4) (3,1) (4,1) y (4,3). Nótese que el patrón se ha recalibrado de una vez en el INM#3 (las calibraciones 1, 3 y 4). En el resto de INMs no se ha repetido la calibración.

Introduciendo las incertidumbres de medición de x se obtiene, en nm²,

La matriz de ponderación se determina mediante la transformación de Choleski, y resulta

Con la matriz Ã = W·A se llega al sistema (18) que determina la solución

En el programa Ej2 (donde A2 es la matriz Ã) los valores concretos son:

es decir: a = 2395,0 nm y b = -3,5 nm/año .

Obsérvese que ha habido una pequeña variación en a (1 nm) pero una significativa variación en b (0,6 nm/año) respecto de la solución obtenida en el ejemplo anterior (ajuste OLS).

El residuo medio cuadrático (19) es s = 1,1833 lo que implica que la matriz  está ligeramente infravalorada, y de acuerdo con el criterio (20),  , y

y con dos cifras significativas .

La figura 3 reproduce el programa Ej2 que realiza estos cálculos en MATLAB/Octave.

Figura 3. Código Ej2 en MATLAB/Octave para ajuste GLS.

7. Ejemplo de ajuste por mínimos cuadrados totales (TLS)

Este último ejemplo amplía los datos de los enunciados 1 y 2 con más información.

Enunciado 3

Se analiza el mismo tipo de ajuste que en los dos ejemplos anteriores, x
(t) = a + b·t, ahora por MM.CC. Totales (TLS) al considerar la siguiente información adicional sobre las medidas de tiempo:

• tienen una incertidumbre de u( t_i ) = 1/ √12 mes
• no existe correlación entre ellas: u( t_i,t_j ) = 0 ∀ i ≠ j

La razón por la cual se toma u( t_i ) = 1/ √12 mes es porque en los certificados de calibración la duración de la calibración del patrón es de un mes. Asumiendo una distribución uniforme con una amplitud de un mes, la incertidumbre típica resulta con el valor indicado. Expresada en años, es
u ( t_i ) = 1 / (12√12) año = 0,024 año.

En estas condiciones hay que tener en cuenta parámetros con incertidumbres en las funciones (5) por lo que debe aplicarse un tratamiento con MM.CC. totales (TLS) expresando  de acuerdo con (23).

El programa Ej3 emplea las matrices  donde  es la misma del ejemplo anterior. Las otras dos matrices son:

Correspondiendo la matriz  a la matriz de derivadas parciales respecto de los parámetros q_i(en este caso t_i) de las funciones f_i = a + b·t_i (apartado 1).

El valor de b (que es una incógnita) no se conocería hasta después del ajuste. No obstante, como en el ejercicio anterior, se ha obtenido una estimación para b de -3,5 nm/año. Se utiliza este valor en el cálculo de . A partir de esta matriz, se trabaja ahora como en el ejercicio anterior, obteniéndose

Los resultados finales, con dos cifras significativas, son:

que, en el nivel de redondeo adoptado, coinciden con los obtenidos en el ejemplo anterior dado que el coeficiente es muy pequeño y la contribución  es despreciable frente a C_x en (23). En consecuencia, C_b ≅ C_x .

La figura 4 reproduce el programa Ej3 que realiza los cálculos de este último ejemplo en MATLAB/Octave.

Figura 4. Código Ej3 en MATLAB/Octave para ajuste TLS.