Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

5.- Ajuste de curvas El ajuste de curvas es un proceso mediante el cual, dado un conjunto de N pares de puntos {xi, yi} (siendo x la variable independente e y la dependiente), se determina una función matemática f(x) de tal manera que la suma de los cuadrados de la diferencia entre la imagen real y la correspondiente obtenida mediante la función ajustada en cada punto sea mínima: ⎛ N ⎞ ε = min⎜ ∑ ( y i − f ( xi )) 2 ⎟ ⎝ i ⎠ Generalmente, se escoge una función genérica f(x) en función de uno o más parámetros y se ajusta el valor de estos parametros de la manera que se minimice el error cuadrático, ε. La forma más típica de esta función ajustada es la de un polinomio de grado M; obteniendose para M = 1 un ajuste lineal (o regresión lineal), f ( x) = a 0 + a1 x para M = 2 un ajuste parabólico,

f ( x) = a 0 + a1 x + a 2 x 2 etc..

Por otro lado podemos tener un conjunto de datos multidimensionales; es decir, un conjunto de N puntos en un espacio k+1-dimensional del tipo { xi(1), xi(2), ..., xi(k),... yi,}. La función que ajustaremos a estos puntos será una función de k variables y = f(x(1), x(2),..., x(k)) El ajuste multidimensional más sencillo es considerar una dependencia lineal de la función respecto a cada una de las variables de que depende; es decir, ajustando una funcion del tipo

f ( x (1) , x ( 2 ) ,..., x ( k ) ) = a0 + a1 x (1) + a2 x ( 2) + ... + ak x ( k ) de tal manera que se minimice el error cuadrático respecto al conjunto de parámetros {a0, a1,..,ak}. Es lo que se conoce como ajuste o regresión multilineal. En esta sección veremos que el ajuste lineal, el de un polinomio de grado M y el ajuste multilineal se pueden expresar dentro de un mismo formalismo de manera que las respectivas soluciones al problema se pueden determinar mediante algoritmos prácticamente análogos. 20

Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

5.1. Regresión lineal Supongamos que tenemos un conjunto de N puntos en el plano {xi, yi}. El objetivo es determinar la ecuación de la recta tal que minimiza el error cuadrático ⎛ N ⎞ ⎛ N ⎞ ε = min⎜ ∑ ( yi − yicalc ) 2 ⎟ = min⎜ ∑ ( yi − a0 − a1 xi ) 2 ⎟ ⎝ i ⎠ ⎝ i ⎠ respecto a los parámetros a0 (ordenada al origen) y a1 (pendiente). Matemáticamente: N N N ∂ε = ∑ 2( yi − a0 − a1 xi ) = 2∑ yi − 2 Na0 − 2a1 ∑ xi = 0 ∂a0 i i i N N N N ∂ε = ∑ 2( yi − a0 − a1 xi ) xi = 2∑ xi yi − 2a0 ∑ xi − 2a1 ∑ xi2 = 0 ∂a1 i i i i

Simplificando las ecuaciones anteriores vemos que se debe cumplir que N

N

i

i

Na0 + a1 ∑ xi = ∑ yi N

N

N

i

i

i

a0 ∑ xi + a1 ∑ xi2 = ∑ xi yi

o bien, dividiendo ambas ecuaciones por el numero total de puntos e introduciendo valores medios a0 + a1 x = y a0 x + a1 x 2 = xy

En forma matricial, podemos escribir

⎛ 1 x ⎞⎛ a0 ⎞ ⎛ y ⎞ ⎟⎜ ⎟⎟ = ⎜⎜ ⎟⎟ ⎜⎜ 2 ⎟⎜ ⎝ x x ⎠⎝ a1 ⎠ ⎝ xy ⎠ por lo que determinar los parámetros de la recta se resume a resolver el sistema de ecuaciones lineales de dos ecuaciones y dos incógnitas anterior.

21

Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

Algoritmo general matricial

Consideremos ahora el mismo problema desde otra perspectiva. Vamos a suponer que los N puntos pueden pasar exactamente por la recta que buscamos. En este caso, plantearíamos el siguiente sistema N de ecuaciones con 2 incógnitas ⎧a0 + a1 x1 = y1 ⎪a + a x = y ⎪ 0 1 2 2 ⎨ ⎪... ⎪⎩a0 + a1 x N = y N

o bien, en forma matricial ⎛ y1 ⎞ ⎛ 1 x1 ⎞ ⎜ ⎟ ⎜ ⎟ a ⎛ ⎞ x 1 ⎜y ⎟ 0 ⎜ 2 ⎟ ⎜⎜ ⎟⎟ = ⎜ 2 ⎟ ⎜ ... ... ⎟⎝ a1 ⎠ ... ⎜ ⎟ ⎜ ⎟ ⎜y ⎟ ⎝ 1 xN ⎠ ⎝ N⎠

⎛ 1 x1 ⎞ ⎟ ⎜ 1 x2 ⎟ ⎜ A= ⎜ ... ... ⎟ ⎟ ⎜ ⎝ 1 xN ⎠

r r A⋅ x = y

r ⎛a ⎞ x = ⎜⎜ 0 ⎟⎟ ⎝ a1 ⎠

⎛ y1 ⎞ ⎜ ⎟ r ⎜ y2 ⎟ y =⎜ ⎟ ... ⎜ ⎟ ⎜y ⎟ ⎝ N⎠

Como ya hemos visto, una manera directa de resolver los sistemas de ecuaciones expresados en forma matricial es la de multiplicar por la izquierda a ambos lados de la igualdad por la inversa de la matriz de coeficientes. sin embargo, en este caso, al tener mas ecuaciones que incógnitas, la matriz de coeficientes, A, no es una matriz cuadrada (tendrá dimensión N × 2) por lo que no esta definida su inversa. Una posible estrategia a seguir es multiplicar la ecuación matricial anterior por la transpuesta de la matriz de coeficientes (el producto de una matriz por su transpuesta es siempre una matriz cuadrada y simétrica) de manera que tendremos r r AT ⋅ A ⋅ x = AT ⋅ y

S = AT ⋅ A

r zr = AT ⋅ y

El sistema de N ecuaciones y dos incógnitas inicial lo hemos condensado en otro sistema de dos ecuaciones y dos incógnitas dado por

r r S⋅x = z

Ahora bien, la matriz S es una matriz cuadrada de dimensión 2×2 y simétrica por lo que es invertible así que podemos escribir r r S −1 ⋅ S ⋅ x = S −1 ⋅ AT ⋅ y

y por tanto, el vector que buscamos nos quedaría de la forma

r r x = S −1 ⋅ AT ⋅ y

22

Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

Ahora bien, representa esta estrategia una solución diferente al problema de la que hemos deducido anteriormente mediante minimización del error cuadrático? Veremos a continuación que no es el caso.

A partir de la estructura de la matriz A, el producto de su transpuesta por ella misma resulta en ⎛1 AT ⋅ A = ⎜⎜ ⎝ x1

1 x2

⎛ 1 x1 ⎞ ⎛ ⎟ ⎜ N ⎜ ... 1 ⎞⎜ 1 x2 ⎟ ⎜ ⎟ = ... x N ⎟⎠⎜ ... ... ⎟ ⎜ N ⎟ ⎜ ∑ xi ⎜ ⎝ 1 xN ⎠ ⎝ i =1

⎞ ⎟ i =1 ⎟ N ⎟ xi2 ⎟ ∑ i =1 ⎠ N

∑x

i

Del mismo modo, el producto de AT por el vector de términos independientes queda r ⎛1 r z = AT ⋅ y = ⎜⎜ ⎝ x1

⎛ y1 ⎞ ⎛ N ⎜ ⎟ ⎜ ∑ yi ⎞⎟ ... 1 ⎞⎜ y2 ⎟ ⎜ i=1 ⎟ ⎟ = ⎟ ... x N ⎟⎠⎜ ... ⎟ ⎜ N ⎜ ⎟ ⎜ ∑ xi yi ⎟ ⎜ y ⎟ ⎝ i =1 ⎠ ⎝ N⎠

1 x2

por lo que la ecuación matricial 2×2 anterior la podemos escribir mas explicitamente como ⎛ ⎜ N ⎜ ⎜ N ⎜ ∑ xi ⎝ i =1

⎛ N ⎞ ⎜ ∑ yi ⎟ ⎞ i =1 ⎟. ⎜⎜ 0 ⎟⎟ = ⎜ Ni =1 N ⎜ ⎟ 2 ⎟⎝ a1 ⎠ xi ⎟ ∑ ⎜ ∑ xi yi ⎟ i =1 ⎝ i =1 ⎠ ⎠ ⎞

N

∑ x ⎟⎟⎛ a i

Dividiendo cada ecuación por N y utilizando la notación típica para el valor medio llegamos a

⎛1 ⎜⎜ ⎝x

x ⎞⎛ a0 ⎞ ⎛ y ⎞ ⎟⎜ ⎟ = ⎜ ⎟ x 2 ⎟⎠⎜⎝ a1 ⎟⎠ ⎜⎝ xy ⎟⎠

que es exactamente el mismo sistema de ecuaciones al que habíamos llegado anteriormente imponiendo la condición de mínimo error cuadrático.

Así pues, podemos plantear el siguiente algoritmo para el ajuste lineal a) Lectura de los N pares de valores {xi, yi} y construcción de las matriz A y vector y. b) Construcción de la transpuesta de la matriz A: AT c) Construcción de la matriz S mediante el producto matricial ATA d) Construcción del vector z mediante el producto matriz por vector ATy e) Inversión de matriz S: S-1 f) Producto matriz por vector S-1 z para obtener el vector de soluciones final.

23

Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

A priori puede parecer un algoritmo demasiado complicado para un problema tan simple como la regresión lineal, para el que existen fórmulas directas de implementación sencilla. Sin embargo, vamos a ver que podemos extender este algoritmo de manera trivial a otro tipo de ajustes.

24

Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

5.2. Ajuste polinómico por mínimos cuadrados. De manera análoga al caso lineal, el objetivo es determinar la ecuación del polinomio de grado M que minimiza el error cuadrático ⎛

N

⎞

⎛

N

⎞

⎝

i

⎠

⎝

i

⎠

ε = min⎜ ∑ ( yi − yicalc ) 2 ⎟ = min⎜ ∑ ( yi − a0 − a1 xi − a2 xi2 ... − aM xiM ) 2 ⎟ respecto a los parámetros M +1 parámetros a0 , a1 ,... aM . Por ejemplo, para un ajuste parabólico (M = 2), la condición de mínimo del error cuadrático lleva a las ecuaciones siguientes: N N N N ∂ε = ∑ 2( yi − a0 − a1 xi − a2 xi2 ) = 2∑ yi − 2 Na0 − 2a1 ∑ xi − 2a2 ∑ xi2 = 0 ∂a0 i i i i N N N N N ∂ε = ∑ 2( yi − a0 − a1 xi − a2 xi2 ) xi = 2∑ xi yi − 2a0 ∑ xi − 2a1 ∑ xi2 − 2a2 ∑ xi3 = 0 ∂a1 i i i i i N N N N N ∂ε = ∑ 2( yi − a0 − a1 xi − a2 xi2 ) xi2 = 2∑ xi2 yi − 2a0 ∑ xi2 − 2a1 ∑ xi3 − 2a2 ∑ xi4 = 0 ∂a2 i i i i i

Procediendo de manera análoga al caso lineal llegamos a que la determinación de los parámetros del polinomio pasa por la resolución de un sistema de ecuaciones de la forma:

⎛1 ⎜ ⎜x ⎜⎜ 2 ⎝x

x 2 ⎞⎛ a0 ⎞ ⎛ y ⎞ ⎟⎜ ⎟ ⎜ ⎟ x3 ⎟⎜ a1 ⎟ = ⎜ xy ⎟ ⎟⎜ ⎟ ⎜ 2 ⎟ x 4 ⎟⎠⎝ a2 ⎠ ⎝ x y ⎠

x x2 x3

Para el caso general de un polinomio de grado M ya podemos intuir que la solución vendrá dada por un sistema de ecuaciones lineales de dimensión (M+1) ×(M+1) de la forma

⎛ 1 ⎜ ⎜ x ⎜ ⎜ ... ⎜ xM ⎝

x x2 ... x M +1

x M ⎞⎛ a0 ⎞ ⎛ y ⎞ ⎟⎜ ⎟ ⎜ ⎟ M +1 a xy ... x ⎟⎜ 1 ⎟ ⎜ ⎟ ⎟⎜ ... ⎟ = ⎜ ... ⎟ ... ... ⎟⎜ ⎟ ⎜ ⎟ M ... x 2 M ⎟⎠⎜⎝ aM ⎟⎠ ⎜⎝ x y ⎟⎠ ...

Una manera de implementar el ajuste polinomial general es la de escribir un

25

Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

subprograma (función) general de la forma ⎞ 1⎛ N r r f ( N , p, q, x , y ) → x p y q ≡ ⎜ ∑ xip yiq ⎟ N ⎝ i=1 ⎠ para determinar todos los promedios de potencias de x y productos de éstas por y que aparecen en la matriz de coeficientes y de términos independientes. r r f ( N ,1,0, x , y ) → x r r f ( N ,1,1, x , y ) → xy r r f ( N ,3,0, x , y ) → x 3

Sin embargo, es mucho mas sencillo aplicar el algoritmo general matricial descrito anteriormente para el caso lineal. Ahora, para el ajuste de un polinomio de grado M plantearíamos un sistema N de ecuaciones con M+1 incógnitas (incluyendo el término independiente) ⎧a0 + a1 x1 + a2 x12 + ... + aM x1M = y1 ⎪ 2 M ⎪a0 + a1 x2 + a2 x2 + ... + aM x2 = y2 ⎨ ⎪... ⎪a + a x + a x 2 + ... + a x M = y 2 N M N N ⎩ 0 1 N o bien, en forma matricial ⎛1 x1 ⎜ ⎜ 1 x2 ⎜1 ... ⎜ ⎜1 x N ⎝

x12 x22 ... x N2

... x1M ⎞⎛ a0 ⎞ ⎛ y1 ⎞ ⎟⎜ ⎟ ⎜ ⎟ ... xNM ⎟⎜ a1 ⎟ ⎜ y2 ⎟ = ... ... ⎟⎟⎜ ... ⎟ ⎜ ... ⎟ ⎜ ⎟ ⎜ ⎟ ... xNM ⎟⎠⎜⎝ aM ⎟⎠ ⎜⎝ y N ⎟⎠

r r A⋅ x = y

Así pues, una vez construida la matriz A de dimensión N × M + 1 y el vector de coeficientes y a partir de los N pares de puntos {xi, yi}, simplemente aplicamos los pasos b) a f) descritos anteriormente.

Un algoritmo general de ajuste de una función polinómica debería solicitar únicamente: a) el numero de total de pares de puntos de que disponemos, N, y b) el grado del polinomio que se pretende ajustar, M

A partir de aquí, deberemos definir en nuestro programa la matriz de coeficientes A de dimensión N × M + 1 y el vector de coeficientes y de dimensión N. Los diferentes pasos

26

Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

del algoritmo requieren el uso de la matriz transpuesta (de dimensión M + 1 × N), la matriz resultante del producto de la transpuesta de A por ella misma y su inversa (de dimensión M+1 × M+1 ambas) y el vector transformado por la transpuesta de A (de dimensión M + 1), así como el vector de soluciones, de dimensión M + 1.

5.3. Regresión multilineal El algoritmo matricial anterior se puede adaptar de manera sencilla para el caso de la regresión multilineal. En este caso disponemos de un conjunto de datos multidimensionales; es decir, un conjunto de N puntos en un espacio k+1-dimensional del tipo { xi(1), xi(2), ..., xi(k),... yi,}. La función que ajustaremos a estos puntos será del tipo f ( x (1) , x ( 2 ) ,..., x ( k ) ) = a0 + a1 x (1) + a2 x ( 2) + ... + ak x ( k )

Si planteamos el sistema de ecuaciones N ecuaciones con k incógnitas para la solución exacta tendremos ⎧a0 + a1 x1(1) + a2 x1( 2 ) + ... + ak x1( k ) = y1 ⎪ (1) ( 2) (k ) ⎪a0 + a1 x2 + a2 x2 + ... + ak x2 = y2 ⎨ ⎪... ⎪a + a x (1) + a x ( 2 ) + ... + a x ( k ) = y 2 N k N N ⎩ 0 1 N o bien, en forma matricial ⎛1 x1(1) ⎜ ⎜1 x2(1) ⎜1 ... ⎜ ⎜1 x (1) N ⎝

x1( 2) ... x1( k ) ⎞⎛ a0 ⎞ ⎛ y1 ⎞ ⎟⎜ ⎟ ⎜ ⎟ x2( 2) ... x N( k ) ⎟⎜ a1 ⎟ ⎜ y2 ⎟ = ... ... ... ⎟⎟⎜ ... ⎟ ⎜ ... ⎟ ⎜ ⎟ ⎜ ⎟ x N( 2) ... x N( k ) ⎟⎠⎜⎝ ak ⎟⎠ ⎜⎝ y N ⎟⎠

r r A⋅ x = y

donde las columnas de la matriz de coeficientes A corresponden simplemente a los N valores de entrada de cada dimensión o categoría (incluyendo una columna extra de “1” relativa al término independiente). En este punto podemos aplicar el algoritmo matricial general exactamente de la misma manera que lo haríamos para ajustar un polinomio de grado k a un conjunto de N pares

de puntos.

27

Técnicas Computacionales, Curso 2007-2008. Pedro Salvador

5.3. Coeficiente de determinación El coeficiente de determinación, R2, definido entre 0 y 1, nos da una idea de la bondad del ajuste, de manera que para valores cercanos a 1 el ajuste es perfecto mientras que para valores cercanos a cero indica inexistencia de relación entre x e y con el modelo de ajuste propuesto. El coeficiente R2 viene dado por la relación entre la varianza de los datos explicada con el modelo y la varianza de los datos experimentales. En concreto

∑ (y N

R2 =

calc i

−y

)

2

i

∑ (y N

− y)

2

i

i

donde y representa el valor medio de los valores de la variable independiente e yicalc los valores calculados para cada punto usando el modelo ajustado a los datos. La implementación computacional de este índice es muy sencilla una vez ajustado el modelo tras resolver el sistema de ecuaciones que plantea el algoritmo matricial general.

En el caso de la regresión lineal, el coeficiente de determinación tiene la misma expresión que el coeficiente de regresión r2, que indica también cómo de correlacionadas estadísticamente están las variables aleatorias x e y . Es importante ver que ambos coeficientes tiene significados e interpretaciones diferentes y que, salvo en el caso de la regresión lineal, no coinciden. Así, se puede comprobar que para el caso de la regresión lineal este índice coincide con el coeficiente de regresión, definido a partir de la relación entre la covarianza de las variables aleatorias x e y y el producto de la raíz cuadrada de las varianzas individuales (desviación típica) de ambas variables, con el fin de obtener un parámetro adimensional r=

(x

xy − x y 2

−x

2

)(y

2

−y

2

)

=

Cov( XY )

σ xσ y

Como se puede ver, el valor de este coeficiente es independiente del modelo ajustado, ya que únicamente indica la relación estadística entre el conjuntos de datos.

28