Metrología y enseñanza

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

Angel Sanchez
Ángel Mª Sánchez Pérez
Laboratorio de Metrología y Metrotecnia, ETSII
Universidad Politécnica de Madrid
Jesus de Vicente y Oliva
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 x1, x2... xm los resultados proporcionados directamente por el sistema de medida (variables de entrada) e y1, y2... yn 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:

(1)

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 y1, y2... yn 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 e1, e2... em:

(2)

El vector y1, y2... yn se elige de tal manera que el vector de residuos resultante e1, e2... em 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:

(3)

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

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

(4)

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

(5)

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

Figura 1 Modelo general de ajuste por MM.CC
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, y1, y2... yn, (magnitudes de salida) a partir de m datos de medida, x1, x2... xm, (magnitudes de entrada), conociendo, además, un conjunto de parámetros complementarios, q1, q2... qp. Al minimizar la suma cuadrática de los residuos ei, se obtienen las ecuaciones necesarias para determinar los parámetrosy1, y2... yn. 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

(6)
(7)

donde

(7a)

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 fi 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 fi ( y1, y2... yn ) no contienen ningún parámetro que proceda de un resultado de medida con incertidumbre. Dicho de otro modo, las funciones f(y1, y2... yn ) sólo contienen las incógnitas y1, y2... yn a determinar y parámetros o coeficientes conocidos sin incertidumbre.

    En la práctica, las incertidumbres de los parámetros presentes en fi ( y1, y2... yn ) no tienen por qué ser estrictamente nulas. Basta con que su contribución final a las incertidumbres de las variables de salida y1, y2... yn sea despreciable.

  • Las variables de entrada (medidas conocidas), es decir, las x1, x2... xm cumplen con lo siguiente:

    • Todas ellas poseen la misma incertidumbre típica, estimada por
      u( x1 ) = u( x2 ) = ... = u( xm ) = ux

    • Todas ellas son estadísticamente incorreladas. Es decir, todos los coeficientes de correlación r ( xi,xj ) son nulos (para ij ). O dicho de otro modo, todas las covarianzas son nulas,u( xi,xj ) = 0 .

      Es decir, la matriz de covarianzas Cx = cov(x) del vector de variables aleatorias x = [x1, x2... xm]T responde a

      (8)
    • La estimación inicial de ux 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

(9)

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

La suma cuadrática de los residuos es

(10)

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

(11)
o sea
(12)

sistema cuya solución es

(13)

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

(14)

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

(15)

Puede demostrarse que s2 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 s2 para ux2.

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 (s2 elevada) no es consistente con unos valores mucho menores de la incertidumbre estimada por los laboratorios participantes y puede adoptarse el criterio

(16)

Una justificación estadística más detallada, adoptando normalidad para las variables implicadas, se basa en que una buena estimación de ux determina que Q / ux2 se distribuya según una χ2 con m-n grados de libertad. Es posible obtener un valor crítico χc2 por debajo del cual la probabilidad de obtener valores de χ2 sea muy alta (por ejemplo del 95 %). A este valor le corresponde otro, sc2, según (15). Cuando s sc, se acepta como un suceso dentro del 95 % de sucesos esperables en la hipótesis de que ux 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 fi ( y1, y2... yn ) no contienen ningún parámetro que proceda de un resultado de medida con incertidumbre. Dicho de otro modo, las funciones fi ( y1, y2... yn ) sólo contienen las incógnitas y1, y2... yn a determinar y parámetros o coeficientes conocidos sin incertidumbre.

    En la práctica las incertidumbres de los parámetros presentes en fi ( y1, y2... yn ) no tienen por qué ser estrictamente nulas. Basta con que su contribución final a las incertidumbres de las variables de salida y1, y2... yn sea despreciable.

  • Las variables del término independiente, es decir, las x1, x2... xm no suelen poseer la misma incertidumbre y es posible que, también, existan correlaciones no nulas entre ellas: para alguna pareja xi, xj  r ( xi, xj ) ≠ 0 lo cual es equivalente a u ( xi, xj ) ≠ 0 .

  • Se supone conocida la matriz de covarianzas (1) del vector de variables aleatorias x = [x1, x2... xm]T que se presenta multiplicada por un factor positivo, s2, desconocido, que permitirá corregir en el ajuste una posible infravaloración de (1) .

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

(17)

Debiendo ser

(17)

La matriz W puede determinarse por varios procedimientos. Uno de ellos es mediante la descomposición de Cholesky (W = L-1 siendo Cx = LLT y L triangular inferior) sobre la matriz de covarianzas de los valores medidos, Cx = 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

(18)

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:

(19)

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

(20)

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

(21)

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

4. Mínimos cuadrados totales (TLS)

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

  • En este caso, las funciones fi ( y1, y2... yn ) contienen una serie de parámetros q1, q2... qp, que se conocen con incertidumbres no despreciables por lo que las incertidumbres están presentes tanto en el término x1, x2... xm como en las funciones fi ( y1, y2... yn ). Se trata de analizar el modelo general de la figura 1.

  • Es conocida la matriz de covarianzas Cx = cov(x), quizá en el formato (21a) comentado en el caso GLS.

  • Asimismo, otro dato es la matriz de covarianzas (21b) del vector de parámetros q = [q1, q2... qp]T conocidos que se presenta multiplicada por el mismo factor positivo anterior, s2, desconocido, que permitirá también corregir en el ajuste una posible infravaloración de (21c).

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

(22)

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

(23)

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

(24)

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

(25)

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:

(Tabla_1)

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: (25a)

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 ( y1 por a; y2 por b ) y existe un parámetro, t, cuyos valores se conocen sin incertidumbre. De acuerdo con todo ello, el sistema (1) se escribe:

(26)

La primera de las expresiones (9) , xA·y, se concreta en:

(27)

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 ux2 = s2 . El resultado es:

(28)

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.

Figura 2
Figura 2. Código Ej1 en MATLAB/Octave para ajuste OLS.

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(xi)

  • 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, Cx , del término independiente (vector de medicionesx). En el programa de la figura 3, se trabaja con la matriz (28) del apartado 3, que se designa Cx, para mayorarla al final si (28) (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 nm2,

(28)

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

(28)

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

(28)

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

(28)

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 (28) está ligeramente infravalorada, y de acuerdo con el criterio (20), (28) , y

(28)

y con dos cifras significativas .

(28)

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

Figura 2
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( ti ) = 1/ √12 mes

  • no existe correlación entre ellas: u( ti,tj ) = 0 ∀ i j

La razón por la cual se toma u( ti ) = 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 ( ti ) = 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 (28) de acuerdo con (23).

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

(28)

Correspondiendo la matriz (28) a la matriz de derivadas parciales respecto de los parámetros qi (en este caso ti) de las funciones fi = a + b·ti (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 (28). A partir de esta matriz, se trabaja ahora como en el ejercicio anterior, obteniéndose

(28)

Los resultados finales, con dos cifras significativas, son:

(28)

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 (28) es despreciable frente a Cx en (23). En consecuencia, CbCx .

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

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

Referencias

[1]Björck, A.: Numerical Methods for Least Squares Problems. SIAM, 1996. ISBN 08987713609

[2]Golub, G.H. ; Van Loan, C.F. : Matrix Computations. John Hopkins University Press. 1996. 3ª Edición. ISBN 080185414-8

[3]Strang, G.; Borre, K. : Linear Algebra, Geodesy and GPS. Ed. WellesleyCambridge Press. ISBN 0-9614088-6-3 (1997)

[4]Forbes, A.B. : Parameter Estimation Based on Least Squares Problems. Sección 5 del libro Data Modelling for Metrology and Testing in Metrology Science. Ed. Birkhaüser. ISBN 978-0-8176-4592-2 (2009)

[5]Cox, M.G. ; Forbes, A.B. ; Harris, P.M. : The classification and solution of regression problems for calibration. NPL Report CMSC 24/03.

http://www.npl.co.uk/publications/the-classification-and-solution-of-regression-problems-for-calibration

[6]Nielsen, L. : Evaluation of measurements by the method of least squares. Actas del congreso Algorithms for Approximation IV. Ed J Leversley et al (Huddersfield: University of Huddersfield) pp 170–86.

Disponible en la web del BIPM:

http://www1.bipm.org/utils/common/pdf/JCGM/nielsen_final.pdf

[7]Van Huffel, S.; Lemmerling, P. : Total Least Squares an Error-in-Variables Modelling.. Ed. Springer-Science + Business Media, B.V.

ISBN 978-90-481-5957-4 (2002)

[8]Morkovsky, I.; Van Huffel, S. : Overview of TLS Methods. Signal Processing, vol. 87, pp. 2283–2302, 2007 .

 Disponible en internet: http://eprints.ecs.soton.ac.uk/13855/1/tls_overview.pdf

[9]de Vicente, J. : Ajustes por MM.CC.: aproximación al problema general. Capítulo 8 de los apuntes de la asignatura Modelos para Calibraciones y Medidas del Master en Metrología del CEM-UPM.

[10]de Vicente, J. : Guía para la realización de ajustes por MM.CC. en Metrología. Capítulo 14 de los apuntes de la asignatura Modelos para Calibraciones y Medidas del Master en Metrología del CEM-UPM.

[11]JCGM 102:2011. Evaluation of measurement data – Supplement 2 to the. “Guide to the expression of uncertainty in measurement” –. Extension to any number of output quantities.

Disponible en la web del BIPM (en ingles): http://www.bipm.org/utils/common/documents/jcgm/JCGM_102_2011_E.pdf

[12] Wikipedia. Total Least Squares.

https://en.wikipedia.org/wiki/Total_least_squares

[13] Wikipedia. Generalized Least Squares.

https://en.wikipedia.org/wiki/Generalized_least_squares

Descargar PDF