Mínimas Desviaciones Absolutas

Least Absolute Deviations

28 de mayo de 2015

Resumen Se presentar´ a el m´etodo Least Absolute Deviations tambi´en conocido como L1 norm el cual es una alternativa al bien conocido m´etodo de M´ınimos Cuadrados (Least Squares) para estimar par´ ametros en un modelo de regresi´on en el que se presentan puntos at´ıpicos que no pueden ser eliminados dado que posiblemente estos formen parte del conjunto de datos. Se justificar´a el uso de este m´etodo, exponiendo las ventajas y desventajas que existen en el momento de la aplicaci´on en el modelado. Adem´as se revisar´ a la parte te´ orica, esto es, los pasos a seguir para estimar los par´ametros de un modelo bajo dicha metodolog´ıa. Dentro del art´ıculo, se plantear´a un ejemplo de pr´ actica para estimar una funci´on de regresi´on de un conjunto de datos, haciendo uso del m´etodo descrito, y se llevar´a a cabo una simulaci´on para comparar el m´etodo de m´ınimas desviaciones absolutas y el m´etodo de M´ınimos Cuadrados. Finalmente, se analizar´ an los resultados para extraer conclusiones y mostrar al lector la efectividad del m´etodo en cuesti´ on.

1.

Introducci´ on

La estad´ıstica inferencial es la rama de la estad´ıstica que estudia el conjunto de procedimientos y metodolog´ıas que por medio de inducci´on determina propiedades de una poblaci´ on a partir de una muestra de ´esta. En Regresi´on, se usa la estad´ıstica inferencial debido a que los modelos necesitan estimar los par´ametros del mismo a partir de un conjunto de datos u observaciones que ser´an tomados como la muestra de alguna poblaci´on. El m´etodo de Least Absolute Deviations o L1 norm es una alternativa al m´etodo cl´asico de M´ınimos Cuadrados (Least Squares) en el an´alisis estad´ıstico de los modelos de regresi´on. En vez de minimizar la suma de cuadrados del error (SSE), el m´etodo en cuesti´on se encarga de minimizar la suma de los valores absolutos del error (SAE). Para muchas aplicaciones en donde el conjunto de datos no contiene puntos at´ıpicos, el m´etodo de M´ınimos Cuadrados es una soluci´ on r´ apida y conveniente. La desventaja de M´ınimos Cuadrados es que se basa en la diferencia de cuadrados entre la observaci´on y su valor esperado y esto hace que haya problemas con los valores at´ıpicos. Pese a que existen diferentes definiciones estad´ısticas para un punto at´ıpico, la m´as pr´actica es que un punto at´ıpico es aquel valor que no deber´ıa de estar contenido en el conjunto de datos que ser´a modelado. Una raz´ on por la que se tienen estos valores en los datos es que pudo existir un error al momento de realizarse el experimento. En muchos casos, el conocimiento de los supuestos en los procesos da lugar a la eliminaci´ on de dichos valores. Sin embargo, esto no siempre es posible y lo que aparentemente ser´ıa un valor at´ıpico en realidad podr´ıa ser parte de la poblaci´on que ser´a analizada. En estos casos, es preferible usar el m´etodo de Least Absolute Deviations. Este m´etodo fue introducido 50 a˜ nos antes que el de M´ınimos Cuadrados, en 1757 por Roger Joseph Boscovich. A pesar de su larga historia y los variados trabajos innovadores, este m´etodo no ha sido objeto de estudio formal en contraste con el de m´ınimos cuadrados. 1

Esto debido a la complejidad para el c´alculo de las estimaciones. No fue hasta 1955 cuando Charles, Cooper y Ferguson redujeron la complejidad a un algoritmo computacional y fue entonces que empez´ o la popularidad de este m´etodo en la aplicaci´on. El prop´osito de este art´ıculo es que el lector comprenda que en la pr´actica no siempre se usar´an los m´etodos m´as conocidos para estimar par´ ametros de los modelos y que existen alternativas para ello. En el caso particular descrito anteriormente, se presentar´a el m´etodo de Least Absolute Deviations.

2.

Modelo

En esta secci´ on, se revisar´ a la parte te´orica del modelo. Comenzaremos con el estudio del modelo lineal simple. Consideremos el modelo de regresi´on simple: Yi = β0 + β1 Xi + ξi con i=1,2,...,n (1) donde: (a) Yi es el valor de respuesta en el i-´esimo intento. (b) β0 , β1 son par´ ametros. (3) Xi es una constante conocida. (4) ξi es el t´ermino del error tal que ξi ∼ Laplace(0,2λ 2 ). (5) cov(ξi ,ξj )=0 con i 6= j. Algunos resultados importantes son los siguientes. Notemos que β0 + β1 Xi es una constante y ξi es una variable aleatoria, por lo que Yi es una variable aleatoria. Puesto que E(ξi ) = 0 para toda i, se tiene que E(Yi ) = E(β0 + β1 Xi + ξi ) = β0 + β1 Xi + ξi . De la misma manera, haciendo uso de propiedades b´asicas de la varianza, llegamos a que var(Yi ) = 2λ 2 . Adem´ as, se define a los residuos del modelo como la diferencia entre el valor observado Yi y la estimaci´ on E(Yi ) = Yˆi , esto es, ei = Yi - Yˆi . Se ha considerado que los t´erminos del error siguen una distribuci´on Laplaciana debido a que se asume que el conjunto de datos contiene puntos at´ıpicos. En el m´etodo de M´ınimos Cuadrados, se estima β0 y β1 considerando la suma de los cuadrados de los residuos Σ e2i tan peque˜ na como sea posible. En el m´etodo de Least Absolute Deviations, los par´ ametros se estiman considerando la suma de los valores absolutos de los residuos Σ |ei |, tan peque˜ na como sea posible. Esto es, los par´ametros β0 y β1 se encuentran al minimizar la expresi´on

F (a, b) =

n X

|yi − (b + axi )|

(2)

i=1

A la diferencia entre yi y (b + axi ) se le conoce como desviaci´on del punto (xi ,yi ) a la recta Yˆi = b + aX. Desde este punto en adelante abreviaremos LAD para denotar el m´etodo Least Absolute Deviations y como LS al m´etodo de m´ınimos cuadrados. El concepto de estimaci´ on por LAD no es para nada complejo; en realidad, es m´as simple que el m´etodo LS pues |ei | es una medida m´ as directa del tama˜ no de un residuo que e2i . Sin embargo, cuando es necesario calcular los par´ametros, el m´etodo LAD es mucho m´as complicado. No existe manera anal´ıtica de hacer el c´alculo de estos; en su lugar, se presentar´a un algoritmo para calcularlos. 2

Antes de comenzar con este algoritmo, se analizar´an ciertos aspectos que ser´an considerados al momento de proceder con el m´etodo LAD. Notablemente, la minimizaci´on de la funci´ on (2) est´ a relacionada con el problema de estimaci´on definido como sigue. Sea la variable aleatoria V definida como V = ξ + U donde U es una constante desconocida de locaci´ on y ξ sigue la distribuci´on de Laplace. El estimador de m´axima verosimilitud (EMV) en el conjunto (V1 , V2 , ..., Vn ) es P (3) U ∗ = arg min ni=1 = |Vi − U | La soluci´ on del problema anterior es conocida como M ediana muestral: U ∗ = MED(Vi | i= 1,2,...,n)

(4)

La similitud entre (2) y (3) infiere que, para alg´ un valor a = a0 , el minimizador de (2), digamos b∗a0 que depende de a0 , es un EMV bajo la distribuci´on Laplaciana. Por razones que ser´ an explicadas en la siguiente secci´on, el minimizador de (2) a∗b0 que depende de b0 , dado b = b0 es un EMV para locaci´on bajo la suposici´on del modelo Laplaciano y algunas extensiones. Por lo tanto, una manera intuitiva de resolver el problema de regresi´on por el m´etodo LAD puede ser construido como sigue: primero, mantener uno de los par´ametros b o a constantes, optimizar el otro usando el concepto de EMV, luego alternar el rol de los par´ametros y repetir este proceso hasta que ambos par´ametros converjan.

3.

Algortimo

Consideremos el modelo lineal (1). Si el valor de a es fijo, digamos a0 , la ecuaci´on (2) consta ahora de un solo par´ ametro:

F (b) =

n X

= |yi − (b + a0 xi )|

(5)

i=1

Asumiendo que el t´ermino del error ξi sigue una distribuci´on de Laplace, la ecuaci´on (5) se reduce a un estimador de m´ axima verosimilitud para b. Esto es, se observa la secuencia aleatoria {Yi - a0 Xi } y la meta es estimar el valor dado pero desconocido del par´ametro b. Por ende, de acuerdo con (4), el par´ametro b∗ en este caso puede ser obtenido por la ecuaci´ on b∗ = MED(Yi - a0 Xi | i= 1,2,...,n)

(6)

Si por otro lado, se fija el valor b = b0 , la ecuaci´on 4 se reduce a:

F (a) =

n X

= |yi − (b0 + axi )|

(7)

i=1 0 =|xi || yix−b − a| i

Una vez m´ as, si el t´ermino del error ξi sigue una distribuci´on de Laplace, los ejemplos 0 observados { yix−b } seguir´ an una distribuci´on de Laplace, pero con la diferencia que en i cada observaci´ on en este conjunto tendr´a varianzas diferentes. La raz´on de esto es que para cada constante conocida Xi y que el t´ermino del error tiene media cero, la variable aleatoria ξi /Xi sigue una distribuci´on de Laplace con media cero y la varianza tiene una escala de X12 . Por ello, el par´ ametro a∗ que minimiza (7) puede ser tomada como un i estimador de m´ axima verosimilitud de a y puede ser calculada con la mediana ponderada: 3

a∗ = MED(|xi | ♦

yi −b0 xi

| i= 1,2,...,n)

(8)

donde ♦ es el operador de r´eplica. Este operador funciona de la siguiente manera: para alg´ un entero |Xi |, |Xi | ♦ Yi significa que Yi es replicado una cantidad de |Xi | veces. Si se cae en el caso en que la mediana ponderada no es un entero, sino un n´ umero real, ´esta se calcula de la siguiente manera. Consideremos la mediana ponderada de la forma: Y = MED(Wi ♦ Xi | i=1,2, ...,n) donde Xi son n´ umeros reales, ser´a calculada de la siguiente forma: (i) Calcular W0 como n X W0 = (1/2) Wi i=1

(ii) Ordenar las observaciones como X(1) , X(2) , ..., X(n) con sus pesos correspondientes W[1] , ..., W[n] . (iii) Sumar los pesos comenzando por W[1] y continuando en orden ascendente. (iv) La mediana ponderada es la primera observaci´on X(j) de la cual su peso hace que se satisfaga la siguiente desigualdad: j X

W[i] ≥ W0

i=1

Considerando la mediana ponderada, se construye el siguiente algoritmo para calcular los par´ametros por el m´etodo LAD: (1) Sea k=0. Encontrar un valor inicial a0 para a, tal y como el m´etodo LS lo hace, usando la f´ ormula: a0 = (2)

Pn (X −X)(Yi −Y ) i=1 Pn i 2 i=1 (Xi −X)

Sea k=k + 1. Obtener una nueva estimaci´on de b0 para un valor b1k−1 fijo usando bk = MED(Yi - ak−1 Xi | i= 1,2,...,n)

(3)

Obtener una nueva estimaci´on de a para el valor fijo bk usando ak = MED(|xi | ♦

yi −bk xi

| i= 1,2,...,n)

(4) Una vez que ak y bk no tengan una desviaci´on considerable con respecto a ak+1 y bk+1 , terminar con la iteraci´ on. De otro modo, regresar al paso (2). Dado que la mediana y la mediana ponderada son estimadores de m´axima verosimilitud bajo el criterio de m´ınimas desviaciones absolutas, las funciones obtenidas ser´an decrecientes conforme avance el algoritmo, esto es, F (ak−1 , bk−1 ) ≥ F (ak−1 , bk ) ≥ F (ak , bk ) Por lo tanto, el algoritmo converge.

4.

Ejemplo

Consideraremos los siguientes datos donde se busca estimar una funci´on de regresi´on entre la raz´on de cambio de los nacimientos y el porcentaje de urbanizaci´on. esto se muestra en la siguiente tabla: 4

Se inicia el algoritmo calculando el valor fijo a0 con la f´ormula anal´ıtica de M´ınimos Cuadrados: a0 =

Pn (X −X)(Yi −Y ) i=1 Pn i 2 i=1 (Xi −X)

= -0.4167

Luego, encontraremos b1 usando la mediana del arreglo de valores Yi - a0 Xi . Esto es, se encuentran dichos valores para toda i, se ordenan de menor a mayor y se usa el valor que se encuentre en la mitad de ese arreglo. Ese valor ser´a nuestro b1 , en este caso, b1 = 44.3349. Procederemos a encontrar un nuevo valor para a, llamado a1 en funci´on de b1 , usando la mediana una vez m´ as, pero ahora con el operador r´eplica como se muestra en k la f´ormula |xi | ♦ yix−b . Con esto, llegamos a que a1 = -0.3998. De aqu´ı en adelante lo que i haremos es repetir el proceso tantas veces sea necesario hasta que ak = ak+1 y bk = bk+1 . En la siguiente tabla se muestran los valores obtenidos.

Notemos que tuvieron que pasar 4 iteraciones para llegar al ´optimo.

5.

Simulaci´ on

En esta secci´ on se mostrar´ a la forma en la que se program´o el m´etodo LAD y LS para ver cu´ales son las ventajas y desventajas en la pr´actica. Se utiliz´o la programaci´on en C, donde se generaron 12 mil r´eplicas. Resumiremos los pasos que se siguieron para esta simulaci´on. 1) Se consider´ o el par´ ametro verdadero β0 = 10, mientras que β1 tom´o tres valores: 1, 5, 8.

5

2) Para cada uno de los valores de β1 , se consideraron 50 y 100 observaciones. Se consider´o que las constantes Xi se distribu´ıan de forma Uniforme con rango de 1 a 10. Los errores se distribuyen de forma Normal con media cero y varianza 1. 3) De esas observaciones, se tomaron aleatoriamente 5 y 10 errores y se le sum´o 6 para que estos fueran nuestros puntos at´ıpicos. 4) Despu´es se generaron las observaciones Yi usando las variables ξi y Xi . 5) Luego, se program´ o el m´etodo LS y el m´etodo LAD usando programaci´on lineal, debido a que la idea era minimizar (2). 6) Luego se encontraron los par´ ametros bajo los m´etodos antes mencionados en cada una de las 12 combinaciones (observaciones, puntos at´ıpicos y valores de β1 ) y se corri´o mil veces, es decir, se obtuvieron mil pares (a, b) para LAD y LS en cada combinaci´on. 7) Despu´es de esto se obtuvieron la media y la varianza de cada combinaci´on para β0 y β1 de LS y LAD para poder hacer comparaciones. Como se mencion´ o en lo anterior, se program´o el m´etodo LAD usando programaci´on lineal, usando el siguiente c´ odigo.

Luego, las mil corridas se hicieron con el siguiente c´odigo. Es un caso de los 12 que se consideraron, particularmente son 50 observaciones, 10 puntos at´ıpicos y b1 =1:

La siguiente imagen muestra un ejemplo de la salida para el caso con β1 =8, 50 observaciones y 5 puntos at´ıpicos.

6

Las primeras dos columnas corresponden a b y a obtenidos bajo el m´etodo LAD, mientras que las u ´ltimas dos corresponden a b y a respectivamente, obtenidas con el m´etodo LS. Despu´es de esto, se obtuvo la media y desviaci´on est´andar de los par´ametros bajo las mil repeticiones. Esto para obtener el sesgo y as´ı poder comparar. Esto se resume en la siguiente tabla:

En cuanto a comparaciones, para b1 o b como lo denotamos en este art´ıculo, el sesgo muestra que en 10 de 12 casos es mejor el par´ametro obtenido con LS que el de LAD. Sin embargo, el sesgo es muy cercano a cero, por lo que se podr´ıa concluir que el m´etodo LAD encuentra una buena estimaci´on de dicho par´ametro bajo las condiciones de puntos at´ıpicos en el conjunto de datos. Sin embargo, respecto a la desviaci´on est´andar, se tiene que en 12 de 12 casos la desviaci´ on del par´ametro b es menor usando LS que LAD. Para a, el sesgo fue menor usando el m´etodo LAD que LS en 12 de 12 combinaciones, por lo que se intuye que es mejor dicho m´etodo para su estimaci´on. En cuanto a la varianza, se tiene lo contrario, se obtuvo que la desviaci´on fue menor para LS que para LAD en 12 de 12 casos.

6.

Conclusiones

Para finalizar, en esta secci´ on daremos nuestra opini´on acerca de ambos m´etodos de estimaci´on. Si bien el m´etodo m´ as utilizado en regresi´on para la estimaci´on de par´ametros de una funci´ on de regresi´ on lineal es el de M´ınimos Cuadrados, en la pr´actica no siempre se podr´a implementar. Se analizaron casos cuando en un conjunto de datos existen puntos fuera de lo com´ un, poniendo en peligro la estimaci´on bajo LS. Para estos casos es mejor utilizar el m´etodo LAD. En el ejemplo que utilizamos en este art´ıculo, se vio que si bien es un poco tedioso y complicado realizar el algoritmo, las estimaciones son adecuadas, 7

mientras que en la parte de simulaci´on, se not´o que los errores disminu´ıan si se usaba este m´etodo. Aunque no se pudo completar la simulaci´on, lo realizado nos da pie a pensar que bajo estas circunstancias el m´etodo a elegir es LAD. Supondr´ıamos que durante un problema de implementaci´ on en la vida real este m´etodo tendr´a mejores resultados, pero como se mencion´ o, es una hip´ otesis a la que llegamos por medio de observaci´on. A manera m´as personal, este art´ıculo nos ayud´o a ver que no siempre se cumple lo que un libro dice, en la aplicaci´ on nos podremos encontrar con problemas que surgen por errores humanos o por factores que no est´ an en nuestras manos. Para eso tendremos diferentes opciones de c´omo manejar dichas situaciones y creemos que el haber analizado este m´etodo nos deja las puertas abiertas para analizar diferentes m´etodos para optimizar una operaci´ on definida. Es todo por nuestra parte.

Referencias [1] Dodge, Birkes, Alternative Methods of Regression, Second Edition, Wiley, United States of America, New York, 1993. [2] Arce, Li, ((A Maximum Likelihood Approach to Least Absolute Deviation Regression)), 2004 [3] Kutner, Neter Applied Linear Statistical Models, Fifth Edition, United States of America, 2004 [4] Leroy, Rousseeuw Robust Regression and Outlier detection, Second Edition, Wiley, United States of America, 1987 [5] Rawlings, Pantula Applied Regression Analysis, Second Edition, Springer, United States of America, 1998 [6] Weisberg Applied Linear Regression, Third Edition, Springer, United States of America, 2005 [7] Torres, ((Comparaci´ on de tres m´etodos de regresi´on lineal)), Revista Colombiana de Estad´ıtica, 2001. [8] Chen, Ying, ((Analysis of least absolute deviation)), 2000 [9] Powell, ((Least Absolute Deviations Estimation for the Censored Regression Model)), Massachusetts Institute of Technology, 1983 [10] Dielman, ((Least absolute value regression)), Neeley School of Business, 2003. [11] Pynnonen, Salmi, ((A report on least absolute deviation regression with ordinary linear programming)), University of Vassa, 2001. [12] Bloomfield, Steiger Least absolute deviations: theroy, applications, and algorithms, First Edition, Birkhaus, United States of America, 1983. [13] Farebrother, ((L1 norm and L2 Norm Estimation)), 2013 [14] van Zyl, ((A weighted least squares procceudre to approximate least absolute deviation)), University of the Free State, South Africa, 2002 [15] Huang, Ma, ((Least Absolute Deviations estimation)), University of Iowa, 2005 8

[16] Bai, ((Least Absolute Deviaton estimation of a shift)), University of Cambridge, 1995

9