Inferencia estadística con datos sesgados por longitud

Traballo Fin de Grao Inferencia estadística con datos sesgados por longitud Irene Sanmartín Dopico 2018/2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Inferencia estadística con datos sesgados por longitud Irene Sanmartín Dopico Xullo 2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA iii iv Trabajo propuesto Área de Coñecemento: Estatística e Investigación Operativa Título: Inferencia estadística con datos sesgados por longitud Breve descrición do contido Los datos sesgados por longitud surgen cuando los individuos no tienen la misma probabilidad de forma parte de la muestra. Esto es frecuente en múltiples situaciones de muestreo, por ejemplo si se toman turistas al azar de entre los que se encuentran en el lugar de destino, pues los turistas con estancias más largas tendrán más probabilidad de ser encuestados. Para que los estimadores (por ejemplo, de la duración media de la estancia) no estén sesgados, se han diseñado estimadores especícos con este tipo de datos. Cabe destacar que estos estimadores están relacionados con la media armónica. En este trabajo se revisarán estos estimadores y se ilustrarán con datos simulados y datos reales. Recomendacións Outras observacións Índice general Resumen viii Introducción xi 1. Estimación con datos sesgados por longitud 1 1.1. Estimación con una muestra no sesgada . . . . . . . . . . . . . . . . . . . . 2 1.1.1. Estimación de la función de distribución . . . . . . . . . . . . . . . . 2 1.1.2. Estimación de la media y la varianza . . . . . . . . . . . . . . . . . . 5 1.2. Estimación con una muestra sesgada . . . . . . . . . . . . . . . . . . . . . . 9 1.2.1. Distribución de la variable sesgada . . . . . . . . . . . . . . . . . . . 9 1.2.2. Estimación de la media con una muestra sesgada: media armónica . 13 1.2.3. Estimación de la función de distribución a partir de una muestra sesgadaporlongitud........................... 17 2. Estudio de simulación sobre la media armónica 23 2.1. Modelos considerados en la simulación y sus versiones sesgadas . . . . . . . 23 2.2. Procedimiento general de simulación . . . . . . . . . . . . . . . . . . . . . . 30 2.3. Distribuciones continuas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.4. Distribuciones discretas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3. Regresión con datos sesgados por longitud 45 3.1. Modelo lineal general con datos no sesgados . . . . . . . . . . . . . . . . . . 46 3.1.1. Estimación de los parámetros del modelo : β y σ2. .......... 48 3.1.2. Propiedades de los estimadores . . . . . . . . . . . . . . . . . . . . . 50 3.2. Modelo lineal general con datos sesgados por longitud . . . . . . . . . . . . 52 3.2.1. Estimación de los parámetros del modelo . . . . . . . . . . . . . . . 53 3.2.2. Propiedades de los estimadores . . . . . . . . . . . . . . . . . . . . . 54 3.3. Simulaciones ................................... 55 v vi ÍNDICE GENERAL A. Código R para las simulaciones 57 A.1. Distribuciones continuas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 A.2. Distribuciones discretas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 A.3. Regresión con datos sesgados por longitud . . . . . . . . . . . . . . . . . . . 65 Bibliografía 67 xiv INTRODUCCIÓN Capítulo 1 Estimación con datos sesgados por longitud Los valores que obtenemos de una muestra nos permiten conocer los valores que encontraríamos en una determinada población. Luego cuestiones como la forma en la que se seleccionen los sujetos o el tamaño de la muestra, tienen una gran importancia en el momento de poder determinar en qué medida es representativa de la población a la cual se reere. Una muestra aleatoria o probabilística es aquella en la que todos los sujetos de la población tienen la misma probabilidad de ser escogidos. Las muestras aleatorias aseguran o garantizan mejor el poder extrapolar los resultados, y en ellas tenemos más seguridad de que se encuentren representadas las características importantes de la población, en la proporción que les corresponde. Si el 20% de la población tiene la característica A, (una determinada edad, una determinada situación económica, etc.) podemos esperar que en la muestra también habrá en torno a un 20% con esa característica. Si la muestra no es aleatoria (no probabilística) puede suceder que esté sesgada y que por lo tanto no sea representativa de la población general, porque predominan más unos determinados tipos de sujetos que otros. Por ejemplo, si hacemos una pregunta a los conductores que se paran ante un semáforo, prescindimos de los que utilizan otro medio de transporte, y si hacemos la pregunta a la salida de una estación de metro, prescindimos de los que tienen y utilizan coche particular, etc. En la primera sección de este capítulo veremos cómo estimar la función de distribución 1 2 CAPÍTULO 1. ESTIMACIÓN CON DATOS SESGADOS POR LONGITUD con una muestra no sesgada, así como la media y la varianza poblacionales. También mencionaremos algunos resultados relacionados con propiedades de la función de distribución empírica, como el Teorema de Glivenko-Cantelli. En la segunda sección, vamos a volver a estimar la función de distribución, pero en este caso, considerando muestras sesgadas por longitud. Que una variable esté sesgada por longitud, quiere decir que en lugar de observarla a ella, observaremos otra que toma los mismos valores que la primera pero con probabilidades diferentes. Por lo tanto, veremos que para hacer esta estimación, tendremos que relacionar la función de distribución de una variable sesgada, con la de otra variable no sesgada. También veremos que la media armónica es un estimador de la media poblacional de la variable sesgada por longitud y aproximaremos su distribución haciendo uso del Teorema central del límite y del Método delta. 1.1. Estimación con una muestra no sesgada Vamos a comenzar trabajando con una muestra no sesgada. En esta sección veremos quiénes son los estimadores de la función de distribución, de la media y de la varianza poblacional, así como algunas propiedades. 1.1.1. Estimación de la función de distribución Denición 1.1. Sea X1, ..., Xn una muestra aleatoria de X con función de distribución F . Se dene la función de distribución empírica como la función Fn(x) = 1 n n X i=1 I(Xi≤x), x ∈R, (1.1) que a cada número real x le asigna la proporción de valores observados que son menores o iguales que x . Es inmediato comprobar que Fn así denida es una función de distribución que cumple: 1. Fn(x)∈[0,1] para todo x∈R. 2. Fn es continua por la derecha. 3. Fn es no decreciente. 4. l´ım x→−∞ Fn(x)=0 . 1.1. ESTIMACIÓN CON UNA MUESTRA NO SESGADA 3 5. l´ım x→∞ Fn(x) = 1. En general, Fn se puede considerar como una función de distribución de una variable aleatoria discreta (que podemos llamar Xe ), que asigna probabilidad 1 n a cada uno de los n valores Xi , con i= 1,2, ..., n . Además Fn es un estimador de F . xix1x2 ... xn pi=P(Xe=xi) 1/n 1/n ... 1/n De este modo, Fn asigna a un conjunto A del espacio muestral de X la probabilidad empírica 1 nXI(Xi∈A), que es la proporción de datos Xi que pertenecen a A . Si se ja el valor de x entonces la variable aleatoria I(Xi≤x) toma valores 0,1 , con probabilidad 1−p, p respectivamente. Es decir, es una Bernoulli de parámetro p=F(x) . De ahí se deduce que Fn es una variable aleatoria y que nFn(x) tiene distribución binomial con parámetros n y p=F(x) . Es decir, con media nF(x) y varianza nF(x)(1−F(x)) . Así, E(nFn(x)) = nF(x) =⇒E(Fn(x)) = F(x) (1.2) V ar(nFn(x)) = nF(x)(1 −F(x)) =⇒V ar(Fn(x)) = 1 nF(x)(1 −F(x)). (1.3) Como hemos visto, la esperanza de la función de distribución empírica coincide con la función de distribución teórica, E(Fn(x)) = F(x) . Por lo tanto Fn(x) es un estimador insesgado de F(x) . De lo anterior se sigue que la función de distribución empírica es un proceso estocástico: si consideramos un espacio probabilístico (Ω, A, P) donde están denidas las sucesiones de variables aleatorias {Xn}n≥1 a partir de las cuales deniremos la función de distribución empírica, tenemos que Fn : (Ω, A, P)×(R, B)−→ [0,1] (ω, x)−→ Fn(x)(ω) = 1 nPn i=1 I(Xi≤x)(Xi(ω)) . 4 CAPÍTULO 1. ESTIMACIÓN CON DATOS SESGADOS POR LONGITUD Fijado x , Fn(x)(·) : (Ω, A, P)−→ [0,1] es una variable aleatoria. Fijado ω , Fn(·)(ω) : R−→ [0,1] es una función de distribución. Por lo tanto, la función de distribución empírica es una función de distribución aleatoria. Podemos englobar lo dicho anteriormente en el siguiente teorema. Teorema 1.2. Sea {Xn}n≥1 , sucesión de variables aleatorias independientes e idénticamente distribuidas denidas en el espacio de probabilidad (Ω, A, P) con función de distribución común F . Se denota por Fn la función de distribución empírica obtenida de las n primeras variables aleatorias X1, ..., Xn . Sea x∈R . Se verica lo siguiente: 1. P(Fn(x) = j/n) = n j!F(x)j(1 −F(x))n−j , j= 0, ..., n. 2. E(Fn(x)) = F(x) , V ar(Fn(x)) = 1 nF(x) (1 −F(x)). 3. Fn(x)→F(x) casi seguro. 4. Se tiene √n(Fn(x)−F(x)) pF(x)(1 −F(x)) d −→ Z (1.4) donde Z es una variable aleatoria con distribución normal estándar, y la convergencia es en distribución. Demostración. Los apartados (1) y (2) son consecuencia inmediata del hecho de que nFn(x)∼B(n, p =F(x)) , como hemos visto. Por otro lado, si denimos Yi=I(Xi≤x) , se tiene que Fn(x) es igual a la media aritmética de las variables aleatorias Y1, ..., Yn . Así, el apartado (3) es una aplicación inmediata de la ley fuerte de los grandes números y el apartado (4) es consecuencia del teorema central de límite. Como hemos visto en la asignatura de Probabilidad y Estatística , la convergencia casi segura, implica convergencia en probabilidad, luego del teorema anterior deducimos que Fn(x)p −→ F(x), x ∈R. Pero se tiene un resultado más potente, que es el llamado Teorema de Glivenko-Cantelli. Éste nos dice que Fn(x) converge a F(x) de forma uniforme para todos los valores de x . 1.1. ESTIMACIÓN CON UNA MUESTRA NO SESGADA 5 Teorema 1.3 (Glivenko-Cantelli) . Sea {Xn}n≥1 , sucesión de variables aleatorias independientes e idénticamente distribuidas denidas en el espacio de probabilidad (Ω, A, P) con función de distribución común F . Se denota por Fn la función de distribución empírica obtenida de las n primeras variables aleatorias X1, ..., Xn . Entonces: Supx∈R|Fn(x)−F(x)|c.s. −→ 0. La demostración del teorema anterior puede encontrarse en Vélez y García (1993), p.36. Es importante resaltar que según el apartado (3) del Teorema 1.2, las distribuciones empíricas asociadas a muestras de tamaño n convergen débilmente a la distribución de probabilidad teórica identicada por F , para casi todas las muestras de tamaño innito que se extraigan de F . Esta es una de las consecuencias más importantes de dicho teorema: la distribución empírica converge débilmente con probabilidad 1 a la poblacional cuando el tamaño de la muestra tiende a innito. Esto garantiza la posibilidad de realizar inferencia estadística. Los aspectos probabilísticos de una característica X , medida en una población, se resumen de forma estilizada en una distribución de probabilidad F , la cual puede ser aproximada mediante las distribuciones empíricas Fn obtenidas por muestreo de la población en estudio. El Teorema de Glivenko-Cantelli arma que esas aproximaciones son uniformes en x . Por esta razón el Teorema de Glivenko-Cantelli se llama a veces Teorema Fundamental de la Estadística Matemática: da una fundamentación de la inferencia estadística, cuyo objetivo principal consiste en extraer información sobre F a partir de las observaciones muestrales. 1.1.2. Estimación de la media y la varianza Denición 1.4. Dada una muestra aleatoria simple X1, X2, ..., Xn de X , con E(X) = µ y V ar(X) = σ2. La media muestral es el estadístico obtenido tomando la media aritmética de los elementos de la muestra. La denotaremos mediante ¯ X : ¯ X= n X i=1 1 nXi La media aritmética muestral es un estimador insesgado de la media poblacional, pues: E(¯ X) = E"n X i=1 1 nXi#=1 nE"n X i=1 Xi#=1 n(E[X1] + ... +E[Xn]) = 1 nn µ =µ 6 CAPÍTULO 1. ESTIMACIÓN CON DATOS SESGADOS POR LONGITUD Por otra parte, la varianza de la media muestral sería: V ar(¯ X) = V ar "n X i=1 1 nXi#=1 n2 n X i=1 V ar(Xi) = 1 n2n σ2=σ2 n Como podemos ver, la media muestral en general es diferente a la media poblacional. La magnitud del error que estamos cometiendo dependerá de la distribución de la media muestral, y en particular de su variabilidad. A medida que aumenta el tamaño muestral n , el valor de la varianza de la media muestral decrece, reduciéndose así el error. El requisito mínimo deseable para un estimador es que a medida que el tamaño de la muestra crece, el valor del estimador tienda a ser el valor del parámetro poblacional, propiedad que se denomina consistencia. Denición 1.5. Se dene el error cuadrático medio de un estimador como ECM(ˆ θ) = Eh(ˆ θ−θ)2i. Y diremos que un estimador converge en media cuadrática al verdadero valor del parámetro que está estimando, si su ECM tiende a cero cuando el tamaño de la muestra se va a innito. Además también se cumple que ECM(ˆ θ) = V ar(ˆ θ)+(Sesgo(ˆ θ))2. Denición 1.6. Se dice que un estimador es consistente si converge en probabilidad al valor verdadero del parámetro que pretende estimar, cuando el número de datos de la muestra tiende a innito. La convergencia en media cuadrática implica la convergencia en probabilidad, tal y como vimos en la asignatura de Probabilidad y Estadística . Luego si un estimador converge en media cuadrática es consistente. Además, como el error cuadrático medio es la suma de la varianza y del cuadrado del sesgo y ambos sumandos son no negativos, es inmediato el siguiente resultado: 1.1. ESTIMACIÓN CON UNA MUESTRA NO SESGADA 7 Proposición 1.7. La sucesión {ˆ θn}∞ n=1 de estimadores de un parámetro θ , es consistente en media cuadrática, si y sólo si se cumplen las dos condiciones siguientes: 1. Es asintóticamente insesgado. Esto es, l´ım n→∞ E(ˆ θn) = θ . 2. l´ım n→∞ V ar(ˆ θn)=0 . Luego si el estimador es insesgado, la consistencia en media cuadrática equivale a l´ım n→∞ V ar(ˆ θn)=0 (ya que el sesgo es igual a cero). La media muestral es un estimador consistente de la media poblacional, pues: l´ım n→∞ E(ˆ θ) = l´ım n→∞ E(¯ X) = µ l´ım n→∞ V ar(ˆ θ) = l´ım n→∞ V ar(¯ X) = l´ım n→∞ σ2 n= 0 Denición 1.8. Dada una muestra aleatoria simple X1, X2, ..., Xn de X , con E(X) = µ y V ar(X) = σ2. La varianza muestral se dene como : S2=1 n n X i=1 (Xi−¯ X)2 Vamos a calcular la esperanza de la varianza muestral y comprobar así que es un estimador sesgado de la varianza poblacional. Para ello vamos a realizar algunos cálculos previos sumando y restando la esperanza de la variable aleatoria poblacional. S2=1 n n X i=1 (Xi−¯ X)2=1 n n X i=1 (Xi−¯ X+µ−µ)2=1 n n X i=1 (Xi−µ)−(¯ X−µ)2. Desarrollando el cuadrado: S2=1 n n X i=1 (Xi−µ)2+ ( ¯ X−µ)2−2(Xi−µ)( ¯ X−µ)= =1 n"n X i=0 (Xi−µ)2+n(¯ X−µ)2−2( ¯ X−µ) n X i=1 (Xi−µ)#= 8 CAPÍTULO 1. ESTIMACIÓN CON DATOS SESGADOS POR LONGITUD =1 n"n X i=1 (Xi−µ)2+n(¯ X−µ)2−2( ¯ X−µ)(n¯ X−nµ)#= =1 n"n X i=1 (Xi−µ)2−(¯ X−µ)2# Entonces: E(S2) = 1 nE"n X i=1 (Xi−µ)2−n(¯ X−µ)2#=1 n n X i=1 E(Xi−µ)2−E(¯ X−µ)2= =σ2−E[( ¯ X−µ)2]. Pero por denición y empleando resultados anteriores, se tiene que E[( ¯ X−µ)2] = ECM(¯ X) = V ar(¯ X)+(Sesgo(¯ X))2=σ2 n. Luego, E(S2) = n−1 nσ2. Por lo tanto concluimos que la varianza muestral es un estimador sesgado de la varianza poblacional, como ya habíamos dicho. Ahora vamos a denir un estimador insesgado de la varianza poblacional y después probaremos que lo es. Denición 1.9. Dada una muestra aleatoria simple X1, X2, ..., Xn de X , con E(X) = µ y V ar(X) = σ2. La cuasi-varianza muestral se dene como : Sc2=1 n−1 n X i=1 (Xi−¯ X)2 La cuasi-varianza muestral es un estimador insesgado de la varianza poblacional. Como S2 c=n n−1·S2 , podemos obtener de lo anterior que ES2 c=En n−1S2=n n−1ES2=n n−1n−1 nσ2=σ2. Por lo tanto la cuasi-varianza muestral es un estimador insesgado de la varianza. 1.2. ESTIMACIÓN CON UNA MUESTRA SESGADA 9 1.2. Estimación con una muestra sesgada El sesgo muestral, a veces también llamado efecto de selección o error muestral es una distorsión de un análisis estadístico que surge debido al método de recolección de muestras. Es importante tenerlo en cuenta, pues sino algunas de las conclusiones podrían ser erróneas. En esta sección vamos a estimar la función de distribución de variables sesgadas por longitud, tanto en el caso discreto, como en el continuo. Ahora bien, que una variable esté sesgada por longitud, quiere decir que en lugar de observarla a ella, observaremos otra que toma los mismos valores que la primera, pero con probabilidades diferentes. Por lo tanto, para hacer esta estimación, tendremos que relacionar la función de distribución de la variable sesgada con la de otra variable no sesgada. Veremos que un estimador de la media poblacional de la variable sesgada por longitud, es la media armónica, y aproximaremos su distribución haciendo uso del Teorema central del límite y del Método delta. Aplicando propiedades llegaremos a una relación entre las esperanzas de las variables, sesgada y no sesgada, que vamos a denotar por X e Y respectivamente. Concretamente E(X)≤E(Y) , lo cual resulta intuitivo, pues en una muestra sesgada por longitud los datos con mayor longitud tienen más probabilidad de ser seleccionados que el resto. Finalmente, en el último apartado nos centraremos en la estimación de la función de distribución no ponderada a partir de una muestra sesgada por longitud. También calcularemos su distribución asintótica, haciendo uso nuevamente del Teorema central del límite y del Método delta, pero en el caso multivariante. 1.2.1. Distribución de la variable sesgada Comencemos exponiendo un ejemplo sencillo para ver el signicado de muestra sesgada por longitud. Vamos a considerar una urna con bolas de diferentes tipos y sus frecuencias, y veremos cómo varían dichas frecuencias, en el caso de suponer que la muestra estuviera sesgada. Supongamos que tenemos una urna con 100 bolas: 25 marcadas con un uno, 50 marcadas con un dos, y otras 25 marcadas con un tres. Valores 1 2 3 Frecuencias 25 50 25 16 CAPÍTULO 1. ESTIMACIÓN CON DATOS SESGADOS POR LONGITUD Por lo tanto, √n1 ˆµ−1 µd −→ N(0, V ar(1/Y )). (1.6) Si denimos µ−1=E(1/X) , tenemos que: E1 Y2=Z1 y2·g(y)dy =Z1 y2·y·f(y) µdy =1 µZ1 y·f(y)dy =1 µE1 X=µ−1 µ E1 Y2 =1 µ2 Así, la varianza es: V ar 1 Y=E1 Y2−E1 Y2 =µ−1 µ−1 µ2 =µµ−1−1 µ2. Y sustituyendo en (1.6): √n1 ˆµ−1 µd −→ N0,µµ−1−1 µ2. Pero como lo que queremos ver es la distribución de ˆµ , tenemos que aplicar el método delta considerando h(x)=1/x : √nh1 ˆµ−h1 µ d −→ N 0, h01 µ2 σ2!. (1.7) h1 ˆµ= ˆµ h1 µ=µ 1.2. ESTIMACIÓN CON UNA MUESTRA SESGADA 17 h01 µ2 =µ4=⇒h01 µ2 σ2= (µµ−1−1)µ2 Por lo tanto, sustituyendo en (1.7): √n(ˆµ−µ)d −→ N(0,(µµ−1−1)µ2). Es decir, la distribución límite de ˆµ es una normal con media µ y varianza (µµ−1−1)µ2 n . Así, si queremos calcular un intervalo de conanza asintótico de nivel 1−α para la media: P −zα/2≤ˆµ−µ p(µµ−1−1)µ2≤zα/2!= 1 −α Y haciéndo cálculos llegamos a: Pˆµ−zα/2p(µµ−1−1)µ2≤µ≤ˆµ+zα/2p(µµ−1−1)µ2= 1 −α Luego un intervalo de conanza asintótico para la media, sería: ˆµ−zα/2p(µµ−1−1)µ2,ˆµ+zα/2p(µµ−1−1)µ2 donde zα/2 es el valor de una distribución normal estándar que deja a su derecha una probabilidad de α/2 para un intervalo de conanza de (1 −α) . 1.2.3. Estimación de la función de distribución a partir de una muestra sesgada por longitud En este apartado vamos a estimar la función de distribución no ponderada (cdf) evaluando en z . Posteriormente aproximaremos su distribución, empleando de nuevo el Teorema central del límite y el Método delta. Se puede ver fácilmente que: 18 CAPÍTULO 1. ESTIMACIÓN CON DATOS SESGADOS POR LONGITUD E1 YI{Y≤z} E1 Y=µZz 0 1 yg(y)dy =µZz 0 1 y yf(y) µdy =Zz 0 f(y)dy =F(z) Es decir, F(z) = E1 YI{Y≤z} E1 Y , lo cual se puede estimar mediante la media muestral: ˆ F(z) = 1 nPn i=1 1 Yi P{1 Yi≤z} 1 nPn i=1 1 Yi . Para calcular la distribución límite de ˆ F(z) , vamos a hacer uso del Teorema central del límite y del Método delta en el caso multidimensional. Comencemos recordando estos resultados. Teorema 1.14 (Teorema de Levy-Lindeberg multivariante) . Sea {Xn} una sucesión de vectores aleatorios mutuamente independientes y con la misma distribución que cierto vector X , del que suponemos que tiene vector de medias y matriz de covarianzas, que denotaremos µ=E(X) y Σ = Cov(X, X) . Entonces se tiene que 1 √n n X i=1 (Xi−µ)d −→ N(0,Σ). Teorema 1.15 (Método delta multivariante) . Si {Tn} es una sucesión de vectores aleatorios de dimensión k tales que √n(Tn−θ)d −→ Nk(0,Σ) y g:Rk−→ Rm es una función diferenciable en θ , siendo Dg(θ) la matriz jacobiana de g en θ , entonces √n(g(Tn)−g(θ)) d −→ N(0, Dg(θ) Σ Dg(θ)t) 1.2. ESTIMACIÓN CON UNA MUESTRA SESGADA 19 Denamos ˆ F(z) = Pn i=1 Ui(z) Pn i=1 Vi , donde Ui(z) = (1/Yi si Yi≤z 0 si Yi> z Vi=1 Yi . La distribución asintótica de ˆ F(z) = Pn i=1 Ui(z) Pn i=1 Vi , se encuentra al examinar la distribución conjunta de numerador y denominador. Fijado un i , vamos a calcular la esperanza, la varianza y la covarianza de Ui y Vi . En primer lugar, denotemos como µr(z) = Zz 0 xrf(x)dx . Luego, en particular se tiene que µ0(z) = F(z) . Entonces: Eg{Ui(z)}=µ0(z) µ Eg(Vi) = E1 Yi=1 E(Xi)=1 µ V arg{Ui(z)}=E(U2 i)−(E(Ui))2=µ−1(z) µ−µ2 0(z) µ2=µµ−1(z)−(µ0(z))2 µ2 V arg(Vi) = E(V2 i)−(E(Vi))2=E1 Y2 i−E1 Yi2 =µ−1(z) µ µ2 0(z) µ2=µµ−1(z)−1 µ2 Covg{Ui(z), Vi}=E(Ui·Vi)−E(Ui)·E(Vi) = µ−1(z) µ µ0(z) µ=µ−1 µ−1 µ·µ0 µ= =µµ−1(z)−µ0(z) µ2 Consideremos ahora el vector (An, Bn) , siendo: An=1 n n X i=1 Ui(z) = 1 n n X i=1 1 Yi P{1 Yi≤z} Bn=1 n n X i=1 Vi(z) = 1 n n X i=1 1 Yi 20 CAPÍTULO 1. ESTIMACIÓN CON DATOS SESGADOS POR LONGITUD Tenemos (U1(z), V1), ...(Un(z), Vn), ... una sucesión de vectores aleatorios, mutuamente independientes y con la misma distribución, con vector de medias µ0(z)/µ 1/µ ! , y matriz de covarianzas Σ =     µµ−1(z)−(µ0(z))2 µ2 µµ−1(z)−µ0(z) µ2 µµ−1(z)−µ0(z) µ2 µµ−1(z)−1 µ2     . Aplicando el Teorema central del límite para el caso multivariante, se tiene: √n (An, Bn)− µ0(z)/µ 1/µ !! d −→ N2(0,Σ) Pero como queríamos la distribución asintótica de ˆ F(z) = Pn i=1 Ui(z) Pn i=1 Vi , tenemos que aplicar el Método delta para el caso multivariante. Consideremos la función diferenciable g: (0,+∞)×(0,+∞)−→ R. (a, b)−→ a/b Considerando dicho teorema, se tiene que: √n g(An, Bn)−g µ0(z)/µ 1/µ !! d −→ N2 0, Dg µ0(z)/µ 1/µ !ΣDg µ0(z)/µ 1/µ !t , (1.8) donde Dg(θ) es la matriz jacobiana de g en µ0(z)/µ 1/µ ! . Hacemos los cálculos necesarios: g(An, Bn) = An Bn g µ0(z)/µ 1/µ !=µ0(z) Dg(a, b) =   1 b−a b2  1.2. ESTIMACIÓN CON UNA MUESTRA SESGADA 21 Dg µ0(z)/µ 1/µ !ΣDg µ0(z)/µ 1/µ !t = =µ−µ0(z)µ    µµ−1(z)−(µ0(z))2 µ2 µµ−1(z)−µ0(z) µ2 µµ−1(z)−µ0(z) µ2 µµ−1(z)−1 µ2     µ −µ0(z)µ!= =µ−1(z)(1 −µ0(z)) µ−1(z)(1 −µ0(z)) µ −µ0(z)µ!= =µµ−1(z)−2µµ−1(z)µ0(z)+(µ0(z))2µµ−1(z). Y sustituimos en (1.8). Así llegamos a: √nAn Bn−µ0(z)d −→ N20, µµ−1(z)−2µµ−1(z)µ0(z)+(µ0(z))2µµ−1(z) Luego siempre que todos los términos sean nitos, el numerador y el denominador de ˆ F(z) tienen asintóticamente una distribución normal bivariada. Por lo tanto ˆ F(z) es asintóticamente normalmente distribuido con media µ0(z) = F(z) (1.9) y varianza µµ−1(z)−2µµ−1(z)µ0(z) + µµ−1(µ0(z))2 n. (1.10) 22 CAPÍTULO 1. ESTIMACIÓN CON DATOS SESGADOS POR LONGITUD Capítulo 2 Estudio de simulación sobre la media armónica El presente estudio de simulación tiene por objetivo realizar un análisis de las propiedades de la media armónica. En él consideramos distribuciones continuas y discretas, para cada una de las cuales se simulan mil muestras sesgadas. Se calcula la media armónica de cada muestra sesgada, y luego se hace una comparación de diferentes propiedades de ésta, como son el sesgo, la varianza, el ECM y la varianza asintótica. 2.1. Modelos considerados en la simulación y sus versiones sesgadas En nuestro estudio vamos a emplear los modelos recogidos en el Cuadro 2.1 y sus formas modicadas bajo sesgo por longitud. Es sencillo ver cómo se obtienen algunas de estas formas sesgadas,si empleamos la densidad ponderada de la que hemos hablado en el capítulo anterior: g(y) = yf(y) µ , donde f es la densidad no sesgada y µ su media. Para simular tenemos dos opciones: - En una muestra establecer un sesgo de muestreo. Es decir, un sesgo en la que se recoge una muestra de tal manera que algunos miembros de la población destinada tienen menos probabilidades de ser incluidos que otros. El resultado es una muestra sesgada, una muestra no aleatoria de una población en la que todas las personas o instancias, no tienen las mismas probabilidades de ser seleccionados. - Emplear distribuciones sesgadas por longitud. Las más básicas se encuentran en el 23 24 CAPÍTULO 2. ESTUDIO DE SIMULACIÓN SOBRE LA MEDIA ARMÓNICA Cuadro 2.1. Este es el método que emplearemos en nuestro estudio. Variable aleatoria (V.a.) Función de densidad ó probabilidad V.a. sesgada por longitud Uniforme(0,1) 1 2x Gamma, G(p, a)1 Γ(p)apxp−1e−ax G(p+ 1, a) Beta, B(p, q)1 β(p, q)xp−1(1 −x)q−1B(p+ 1, q) Binomial, B(n, p) n x!px(1 −p)n−x1 + B(n−1, p) Binomial negativa, BN(k, p) x−1 k−1!qk(1 −p)x−k1 + BN(k+ 1, p) Poisson, Po(λ)e−λλx x!1 + Po(λ) Hipergeométrica, H(n, M, N) n x!M(x)(N−M)(n−x)N(n)1 + H(n−1, M −1, N −1) Equiprobable 1 k,∀j∈ {1, ..., k}j k(k+1) 2 ,∀j∈ {1, ..., k} Cuadro 2.1: Algunas distribuciones básicas y sus formas sesgadas por longitud. En primer lugar consideremos las distribuciones continuas Gamma y Beta y las distribuciones discretas Equiprobable y Poisson. Veamos cómo se obtienen sus formas sesgadas. Gamma Este modelo depende de dos parámetros positivos: a y p . Su función de densidad es de la forma f(x) = 1 Γ(p)apxp−1e−ax, x > 0, y su esperanza es p/a . La Γ(p) es la denominada función Gamma de Euler que representa la siguiente integral: Γ(p) = R∞ 0xp−1e−xdx . Además es fácil comprobar que verica que Γ(p) = (p−1) Γ(p−1) . En efecto, si hacemos integración por partes considerando u=xp−1du = (p−1)xp−2 dv =−e−xdx v =−e−x 2.1. MODELOS CONSIDERADOS EN LA SIMULACIÓN Y SUS VERSIONES SESGADAS 25 llegamos a Γ(p) = Z∞ 0 xp−1exdy = (p−1) Z∞ 0 e−xxp−2= (p−1)Γ(p−1). Repitiendo el proceso sucesivamente, se llega a que Γ(p)=(p−1)(p−2)...Γ(1) , para p entero positivo. Por otra parte, por integración directa vemos que Γ(1) = 1 . Entonces, deducimos que Γ(p)=(p−1)! para p entero positivo. Comprobemos que su forma sesgada por longitud coincide con la escrita en el Cuadro 2.1. Sea f(x) la función de densidad de una Gamma(p, a) . Empleando la distribución no ponderada calculada previamente se tiene: g(x) = xf(x) µ=x apxp−1e−ax Γ(p)µ (a) =ap+1 xpe−ax Γ(p)p (b) =ap+1 xpe−ax Γ(p+ 1) . En la igualdad (a) se sustituye µ=p/a y en la igualdad (b) se aplica que Γ(p+ 1) = pΓ(p) . Luego, como queríamos probar, g(x) es la función de densidad de una Gamma(p+1, a) . Beta Su función de densidad para valores 0<x<1 está representada por la expresión f(x) = 1 β(p, q)xp−1(1 −x)q−1, donde β(p, q) es la función beta, denida por β(p, q) = Z1 0 xp−1(1 −x)q−1dy =Γ(p) Γ(q) Γ(p+q). Veamos cuál es la forma sesgada por longitud de esta distribución: 32 CAPÍTULO 2. ESTUDIO DE SIMULACIÓN SOBRE LA MEDIA ARMÓNICA Por último, como hemos dicho, en las tablas también va a aparecer la varianza de la media aritmética simple de la variable no sesgada X , que no es ninguna característica de la media armónica. El motivo de añadir esta columna es hacer una comparación de las muestras anteriores con otras no sesgadas. Nótese que la varianza de la media aritmética simple no debería diferir en gran cantidad de la varianza de la media armónica. Luego tampoco lo debería hacer de la varianza asintótica, ni del ECM. 2.3. Distribuciones continuas Vamos a comenzar simulando muestras sesgadas por longitud de distribuciones continuas. Concretamente una Uniforme, una Gamma y una Beta. Para cada una de ellas consideraremos varios casos, y haremos comparaciones mediante representaciones grácas y tablas. En cada tabla vamos a aproximar las propiedades de la media armónica explicadas en la sección anterior, así como su varianza asintótica y la varianza de la media aritmética simple. Una vez que hallemos estas aproximaciones, haremos una comparación de cómo varían los resultados obtenidos dependiendo del tamaño muestral. También compararemos las diferentes tablas de las distribuciones. Por otra parte, como se ha explicado en el capítulo anterior, en una muestra sesgada por longitud los datos con mayor longitud tienen mayor probabilidad de ser seleccionados. Y claramente el sesgo por longitud afecta más cuando hay valores de x próximos al cero con probabilidades altas. Este es el motivo por el que algunas funciones de densidad, se deforman más que otras al pasar a su forma sesgada, tal y como vamos a ver en esta sección. Uniforme(a,b) La distribución Uniforme es el modelo (absolutamente) continuo más simple. Corresponde al caso de una variable aleatoria que sólo puede tomar valores comprendidos entre dos extremos a y b , de manera que todos los intervalos de una misma longitud (dentro de 2.3. DISTRIBUCIONES CONTINUAS 33 (a, b) ) tienen la misma probabilidad. En nuestro caso trabajaremos con una Uniforme(0,1) y con una Uniforme(9,10) . La primera es un caso particular de la distribución beta, que se corresponde con una beta de parámetros a= 1 y b= 1. Luego su forma sesgada por longitud la podemos deducir del Cuadro 2.1. Sin embargo, en dicha tabla no tenemos una forma sesgada para la Uniforme(a, b) , con a6= 0 y b6= 1 . Por lo tanto para hacer simulaciones con esta distribución, tendremos que hallar la función de distribución sesgada G(y) , y calcular Y=G−1(U) , siendo U una variable con distribución Uniforme(a, b) . Como conocemos la función de densidad no ponderada de una Uniforme, a partir de ella podemos calcular la ponderada o sesgada, g(y) = yf(y) µy =y a+b 2(b−a) =2y b2−a2. Por lo tanto, la función de distribución es G(y) = Zy a 2u b2−a2du =1 b2−a2u2y a=y2−a2 b2−a2, y ∈[a, b]. Sea ahora U∈Uniforme(a, b) , pongamos Y=G−1(U) , entonces: y2−a2 b2−a2=u=⇒y2=a2+u(b2−a2) =⇒Y=pa2+U(b2−a2) Para calcular la varianza asintótica, vamos a obtener µ−1 para las dos Uniformes que hemos considerado. En el caso de Uniforme(0,1) : µ−1=E1 X=Z1 0 1 xf(x)dx =Z1 0 1 x 1 1−0dx = [ln(x)]1 0= +∞, 34 CAPÍTULO 2. ESTUDIO DE SIMULACIÓN SOBRE LA MEDIA ARMÓNICA y en el caso de la Uniforme(9,10) : µ−1=E1 X=Z10 9 1 xf(x)dx =Z10 9 1 x 1 10 −9dx = [ln(x)]10 9=ln(10) −ln(9). En el Cuadro 2.2 se presentan para una distribución Uniforme(0,1) los valores de sesgo, varianza, ECM y varianza asintótica de la media armónica como estimador de la media de X . Se añade una columna con σ2/n , que es la varianza de ¯ X , que sería el estimador natural de la media con datos no sesgados. Así podremos considerar a la media armónica como el estimador óptimo de la media en condiciones ideales. Se muestran los resultados para tamaños muestrales crecientes n= 20,50,100,500 . En el Cuadro 2.3 se presentan los resultados para la Uniforme(9,10). n Sesgo Varianza ECM Varianza asintótica σ2/n 20 2,24 1,11 1,16 Inf 0,03 50 1,19 0,53 0,53 Inf 0,013 100 0,85 0,36 0,37 Inf 0,000694 500 0,40 0,09 0,09 Inf 0,0000139 Cuadro 2.2: Valores de sesgo, varianza, ECM y varianza asintótica de la media armónica en función del tamaño muestral, n , para una Uniforme(0,1). (Los resultados vienen en centésimas ). Se añade el valor σ2/n , también en centésimas. n Sesgo Varianza ECM Varianza asintótica σ2/n 20 1,742 3,932 3,935 4,173 4,16 50 0,120 1,72 1,72 1,66 1,66 100 −0,85 0,83 0,83 0,83 0,83 500 −0,12 0,17 0,17 0,16 0,16 Cuadro 2.3: Valores de sesgo, varianza, ECM y varianza asintótica de la media armónica en función del tamaño muestral, n , para una Uniforme(9,10). (Los resultados vienen en milésimas). Se añade el valor σ2/n , también en milésimas. 2.3. DISTRIBUCIONES CONTINUAS 35 En el Cuadro 2.2 vemos que tanto el sesgo como la Varianza y el ECM convergen a cero cuando el tamaño muestral tiende a innito. Esto muestra que la media armónica es un estimador consistente en media cuadrática. También presentamos la varianza asintótica que hemos obtenido en el capítulo anterior, y vemos que su valor es +∞ . Esto se debe a que µ−1= +∞. Aunque en los demás casos que consideraremos la varianza asintótica se parecerá mucho a la aproximación de la varianza. Si hacemos una comparación de los Cuadros 2.2 y 2.3, vemos que la media armónica para la Uniforme(0,1) tiene mayor sesgo que para la Uniforme(9,10) . El hecho de que la Uniforme(9,10) esté poco sesgada, se ve reejado en que la V ar(¯ X) y la varianza asintótica son muy parecidas. Por otra parte, que la varianza asintótica de la Uniforme(0,1) sea innito es consecuencia de que µ−1=∞ , como vimos anteriormente. Un sesgo pequeño implica que la varianza y el ECM de la media armónica son muy parecidos. Vemos que en la Uniforme(9,10) el sesgo aproximado es más pequeño que en la Uniforme(0,1), y en consecuencia la varianza y el ECM son mucho más parecidos que en el primer caso. En la Figura 2.1 se muestra una representación de la función de densidad de las dos Uniformes de las que hemos hablado y sus formas sesgadas por longitud. Además también añadimos la Uniforme(100,101) . Las funciones de densidad no sesgadas están representadas en rojo, y las sesgadas en verde. Por lo tanto, en la primera columna de la gura podemos ver una Uniforme(0,1) , una Uniforme(9,10) y una Uniforme(100,101) , y en la segunda sus respectivas formas sesgadas. La función de densidad de la Uniforme(0,1) se deforma mucho más al pasar a su forma sesgada que la de la Uniforme(9,10) . Y a su vez, las dos anteriores se deforman más que la función de densidad de la Uniforme(100,101) , en la que a penas vemos cambio al pasar de una forma a otra. Como las funciones de densidad de las Uniformes son rectas, podemos ver de forma numérica lo que se deforman al pasar de la forma no sesgada a la sesgada. Simplemente tenemos que hallar y comparar las pendientes de dichas rectas. 36 CAPÍTULO 2. ESTUDIO DE SIMULACIÓN SOBRE LA MEDIA ARMÓNICA Figura 2.1: Uniforme no sesgada y sesgada Es conocido que cualquier función de densidad de una Uniforme no sesgada es una recta horizontal con pendiente cero. Por otra parte, es sencillo ver que la pendiente de la recta asociada a la forma sesgada de la Uniforme(0,1) es 2, mientras que las de la Uniforme(9,10) y Uniforme(100,101) son 0,1053 y 0,0099 respectivamente. Por lo tanto, la función de densidad de la Uniforme(0,1) es la que más se deforma al pasar a su forma sesgada, seguida de la Uniforme(9,10) , y por último la Uniforme(100,101) . Esto es debido a que en la función de densidad de la Uniforme(0,1) hay valores de x próximos al cero que toman probabilidades altas, mientras que en las Uniforme(9,10) y Uniforme(100,101) las x que toman probabilidades altas están muy alejadas del cero. Estamos hablando de sesgo por longitud, luego a medida que consideramos valores de x próximos al cero que toman probabilidades altas, afecta en mayor medida. Es más, vimos que la función de densidad sesgada de una Uniforme(a, b) , es g(y) = 2y b2−a2 , luego si a y b toman valores muy grandes, la función de densidad sesgada es muy parecida a la no sesgada, pues: g(a) = 2a b2−a2=2a (a+b)(b−a) (a) ≃1 b−a 2.3. DISTRIBUCIONES CONTINUAS 37 g(b) = 2b b2−a2=2b (a+b)(b−a) (a) ≃1 b−a Donde en (a) estamos usando que como a y b toman valores muy grandes, 2a y 2b se parecen mucho a (a+b) . Así, a medida que x toma valores más grandes, la forma sesgada por longitud diere menos de la no ponderada o no segada. Gamma(p,a) En una distribución Gamma, según los valores que tome el parámetro de forma, p , la función de densidad presenta perles muy diversos. Con valores de p menores o iguales que 1, la función de densidad muestra un perl decreciente; en cambio, si p es mayor que la unidad, la función de densidad crece hasta el valor x=p−1 a y decrece a partir de este valor. Vamos a considerar dos casos, uno con p=2 y a=1 y otro con p=8 y a=1. Como se maniesta en el Cuadro 2.1, la forma sesgada por longitud de una Gamma(2,1) es una Gamma(3,1) , y la de una Gamma(8,1) , es una Gamma(9,1). Figura 2.2: Distribuciones Gamma no sesgadas y sesgadas 38 CAPÍTULO 2. ESTUDIO DE SIMULACIÓN SOBRE LA MEDIA ARMÓNICA En la Figura 2.2, podemos ver una Gamma(2,1) y una Gamma(8,1) con sus respectivas formas sesgadas por longitud. Tal y como ocurría con las Uniformes, la Gamma(2,1) se deforma más al pasar a su forma sesgada que la Gamma(8,1) . Esto es debido a que en el primer caso hay valores de x próximos al cero que toman probabilidades altas. A continuación vamos a ver quién es µ−1=E(1/X) , necesaria para el cálculo de la varianza asintótica. µ−1=E1 X=Z∞ 0 1 xf(x)dx (a) =Z∞ 0 1 x 1 Γ(p)apxp−1e−ax (b)dx = =Z∞ 0 (ax)p−2 Γ(p)a2e−ax dx (c) =a2 p−1Z∞ 0 (p−1) (ax)p−2e−ax Γ(p)dx (d) =a2 p−1. En la igualdad (a) sustituimos la función de densidad de una Gamma(p, a) , que es f(x) = 1 Γ(p)apxp−1e−ax . En (b) simplemente hacemos operaciones. En (c) multiplicamos y dividimos por (p−1) para obtener dentro una Γ(p) . Y nalmente en la igualdad (d) simplicamos las dos Γ(p) empleando lo siguiente: Γ(p) = Z∞ 0 xp−1exdx, luego si llamamos u=xp−1du = (p−1)xp−2dx dv =e−xdx v =−e−x Integrando por partes, se obtiene lo siguiente: Γ(p) = −e−xxp−1∞ 0+Z∞ 0 e−x(p−1)xp−2dx = (p−1) Z∞ 0 xp−2e−xdx. En los Cuadros 2.4 y 2.5 se presentan las aproximaciones de las características de la media armónica y la varianza de la media aritmética simple en el caso de considerar una Gamma(2,1) y una Gamma(8,1) . 2.3. DISTRIBUCIONES CONTINUAS 39 n Sesgo Varianza ECM Varianza asintótica σ2/n 20 0.0996 0.1451 0.1550 0.2 0.1 50 0.0332 0.0715 0.0726 0.08 0.04 100 0.0154 0.0375 0.0377 0.04 0.02 500 0.0006 0.0076 0.0076 0.008 0.004 Cuadro 2.4: Valores de sesgo, varianza, ECM y varianza asintótica de la media armónica en función del tamaño muestral, n , para una Gamma(2,1). Se añade el valor σ2/n . n Sesgo Varianza ECM Varianza asintótica σ2/n 20 0.0606 0.482 0.4864 0.4571 0.4 50 0.0157 0.1845 0.1847 0.1828 0.16 100 -0.0014 0.0932 0.0932 0.0914 0.08 500 -0.0085 0.0188 0.0188 0.0183 0.016 Cuadro 2.5: Valores de sesgo, varianza, ECM y varianza asintótica de la media armónica en función del tamaño muestral, n , para una Gamma(8,1). Se añade el valor σ2/n . La aproximación de la varianza de la media armónica a medida que aumenta el tamaño muestral se parece más al ECM en ambos casos. Esto es debido a que el sesgo tiende a cero. Así mismo la varianza asintótica también toma un valor similar a ambos. Si nos jamos en la columna del sesgo de las dos tablas, vemos que la Gamma(2,1) está más sesgada que la Gamma(8,1) . En consecuencia, la varianza asintótica de la Gamma(8,1) se parece más a V ar(¯ X) = σ2/n , que en el primer caso. Es decir, en una Gamma(p, a) un incremento del valor de p produce una disminución del sesgo aproximado de la media armónica. Y en consecuencia la varianza asintótica de la Gamma(p, a) se parece más a V ar(¯ X) = σ2/n que en el caso de otra Gamma(p, a) que tuviera un valor de p inferior. Beta(p,q) La distribución beta es adecuada para variables aleatorias continuas que toman valores en el intervalo (0,1) , lo que la hace muy apropiada para modelar proporciones. En nuestro caso vamos a considerar una Beta(2,1) y una Beta(9,1) . 40 CAPÍTULO 2. ESTUDIO DE SIMULACIÓN SOBRE LA MEDIA ARMÓNICA Comencemos viendo cómo obtener µ−1 , necesaria para el cálculo de la varianza asintótica. E1 X=Z+∞ −∞ 1 xf(x)(a) =1 β(p, q)Z+∞ −∞ 1 xxp−1(1 −x)q−1dx (b) = =1 β(p, q)Z1 0 x(p−1)−1(1 −x)q−1=β(p−1, q) β(p, q) (c) = Γ(p−1) Γ(p) Γ(p−1 + q) Γ(p) Γ(q) Γ(p+q) (d) =p+q−1 p−1. En la igualdad (a) sustituimos la función de densidad de una Beta(p, q) , que es, f(x) = 1 β(p, q)xp−1(1 −x)q−1. En (b) agrupamos términos para conseguir una función β(p−1, q) , sabiendo que β(p, q) = Z1 0 xp−1(1 −x)q−1dy. En (c) empleamos la igualdad demostrada anteriormente: β(p, q) = Γ(p) Γ(q) Γ(p+q). Y nalmente en (d) simplicamos empleando que Γ(p−1) = Γ(p) p−1 . Consideremos los Cuadros 2.6 y 2.7. Nuevamente vemos que en ambos casos el sesgo aproximado se va a cero a medida que aumenta el tamaño muestral, y que la varianza de la media armónica y el ECM se parecen cada vez más. Como en el caso de la Gamma, volvemos a tener que la Beta(2,1) está más sesgada que la Beta(9,1) . Luego la varianza asintótica de la segunda se parece más a V ar(¯ X) que en el primer caso. 2.3. DISTRIBUCIONES CONTINUAS 41 n Sesgo Varianza ECM Varianza asintótica σ2/n 20 0.78 0.52 0.53 0.74 0.015 50 0.13 0.25 0.25 0.29 0.00617 100 0.19 0.13 0.13 0.14 0.00308 500 0.07 0.027 0.027 0.029 0.000617 Cuadro 2.6: Valores de sesgo, varianza, ECM y varianza asintótica de la media armónica en función del tamaño muestral, n , para una Beta(2,1). (Los resultados vienen en centésimas). Se añade el valor σ2/n , también en centésimas. n Sesgo Varianza ECM Varianza asintótica σ2/n 20 -0.74 0.48 0.48 0.5 0.00334 50 -0.7 2 2 0.2 0.00133 100 -0.4 0.1 0.1 0.1 0.000669 500 -0.0396 0.0209 0.209 0.02025 0.000133 Cuadro 2.7: Valores de sesgo, varianza, ECM y varianza asintótica de la media armónica en función del tamaño muestral, n , para una Beta(9,1). (Los resultados vienen en milésimas). Se añade el valor σ2/n , también en milésimas. En la Figura 2.3 se hace una representación de las funciones de densidad correspondientes a ambas distribuciones. En la primera la vemos una Beta(2,1) y su forma sesgada por longitud, y en la segunda una Beta(9,1) y su forma sesgada por longitud. La Beta(2,1) se deforma más al pasar a su forma sesgada por longitud que la Beta(9,1) . Esto se debe a que en la Beta(9,1) los valores de x tienen probabilidad 0 hasta que x es mayor que 0,6 ; mientras que en el otro caso x siempre toma probabilidades positivas, incluso en valores próximos a cero. 48 CAPÍTULO 3. REGRESIÓN CON DATOS SESGADOS POR LONGITUD Sustituyendo cada vector o matriz por un símbolo, llegaríamos a la siguiente expresión abreviada: Y=Xβ +ε, donde Y es el vector de respuestas, la matriz X es una matriz n×p, (cada la representa a un individuo y cada columna a cierta característica) β es el vector de parámetros y ε es un vector que contiene los errores y verica ε∈Nn(0, σ2). La representación del modelo de regresión múltiple permite considerar un contexto más general, denominado modelo lineal general, en el cual caben modelos de regresión con variable explicativa discreta:     Y1 ... Yn     =    x1,1. . . x1,p ... . . . ... xn,1. . . xn,p    ·    β1 ... βp     +    ε1 ... εn     . Bajo el modelo lineal general, X es una matriz no aleatoria, β es un vector de parámetros que hay que estimar y ε∈Nn(0, σ2In) , siendo σ2 la varianza del error, que también hay que estimar, e In la matriz identidad de orden n . Esta última expresión aglutina las suposiciones de homocedasticidad, normalidad e independencia de los errores. 3.1.1. Estimación de los parámetros del modelo : β y σ2. Nuestro problema ahora se va a centrar en la estimación del vector de parámetros β y de la varianza del error σ2. Comenzaremos con la estimación de β y plantearemos este procedimiento desde el método de mínimos cuadrados. Escogeremos como estimador aquel ˆ β donde se alcance m´ın β n X i=1 (Yi−xiβ)2, siendo xi la la i-esima de la matriz del diseño X . 3.1. MODELO LINEAL GENERAL CON DATOS NO SESGADOS 49 En notación matricial, el problema de minimización se puede expresar de manera equivalente así: m´ın β(Y−Xβ)0(Y−Xβ) = m´ın βφ(β), donde φ es la función objetivo. Si derivamos la función objetivo respecto de β , e igualamos a cero, se obtiene lo que se conoce como las ecuaciones normales de regresión, X0Xβ =X0Y, cuya solución es el estimador de β por mínimos cuadrados, dado por: ˆ β= (X0X)−1X0Y. Una vez que se han obtenido los estimadores de los parámetros, ˆ β , se pueden calcular los ajustes o predicciones para los individuos de la muestra, de la siguiente manera: ˆ Yi=xiˆ β, i ∈ {1, ..., n}, o lo que es lo mismo ˆ Y=Xˆ β. Antes de seguir con la estimación de la varianza, vamos a dar una interpretación geométrica a las predicciones obtenidas con ˆ β . Para ello partimos de la expresión matricial ˆ Y=Xˆ β=X(X0X)−1X0Y=HY, donde H=X(X0X)−1X0Y se conoce como matriz hat. Así, las predicciones ˆ Y se obtienen aplicando la matriz hat a las observaciones Y . Además ˆ Y es una combinación lineal de las columnas de X , y dentro de todas las posibles combinaciones, es la que tiene menor distancia a Y . Por ello, podemos decir que la predicción ˆ Y es la proyección de Y sobre el espacio formado por todas las combinaciones lineales de las columnas de X . Además, también podemos decir que H es una matriz de proyección sobre dicho espacio. Como corresponde a una matriz de proyección, es una matriz n×n simétrica, idempotente y de rango p. Vamos cerrar esta sección con la estimación de la varianza. 50 CAPÍTULO 3. REGRESIÓN CON DATOS SESGADOS POR LONGITUD Es conocido que los residuos se denen como la diferencia entre las observaciones y las predicciones, esto es, ˆεi=Yi−ˆ Yi=Yi−xiˆ β, i ∈ {1, ..., n}. Por lo tanto, se puede formar un vector de residuos: Y−ˆ Y= (In−H)Y=MY, donde M=In−H se conoce como la matriz generadora de residuos. Es fácil ver que M es una matriz simétrica, idempotente, de rango (n−p) y ortogonal a la matriz hat, es decir, MH = 0. Como los errores no se observan, para estimar su varianza σ2 , emplearemos los residuos que acabamos de denir. Así, el estimador de la varianza del error sería ˆσ2=1 n−p n X i=1 ˆεi2=1 n−p n X i=1 (Yi−xiˆ β)2=RSS n−p. Hemos extendido la notación RSS para la suma residual de cuadrados, y hemos empleado como denominador (n−p) para conseguir un estimador insesgado de la varianza del error. 3.1.2. Propiedades de los estimadores En este apartado vamos a deducir las propiedades de los estimadores de los parámetros asociados al modelo lineal genera, es decir de ˆ β y de ˆσ2. Estas propiedades serán su comportamiento en media (veremos que son insesgados), su varianza y su distribución e independencia. En primer lugar vamos a ver que ˆ β es un estimador insesgado de β . Para ello es suciente emplear que E(ε)=0. E(ˆ β)=(X0X)−1X0E(Y)=(X0X)−1X0Xβ =β Vamos a obtener la matriz de covarianzas del vector de aleatorio ˆ β . Basta tener en cuenta la hipótesis de homocedasticidad, es decir , V ar(Yi) = σ2 . Por lo que la matriz de covarianzas del vector de respuestas se puede escribir como Cov(Y, Y ) = σ2In. Cov(ˆ β, ˆ β) = Cov (X0X)−1X0Y, (X0X)−1X0Y= = (X0X)−1X0σ2InX(X0X)−1=σ2(X0X)−1 3.1. MODELO LINEAL GENERAL CON DATOS NO SESGADOS 51 Para la estimación de cada coeciente del modelo, que es una componente del vector β , se toma la componente correspondiente de ˆ β , que es un estimador insesgado. Y su varianza se obtiene multiplicando σ2 por el elemento correspondiente de la diagonal de (X0X)−1. Vamos a mencionar un par de propiedades que utilizaremos en la demostración del siguiente teorema. Propiedad 1. Si X∈Nm(µ, Σ) y C es una matriz p×m de rango p , entonces CX ∈Np(Cµ, CΣC0). Propiedad 2. Sea X1, ..., Xm∈N(0, σ2) una muestra aleatoria simple formada por m observaciones independientes de una misma distribución normal de media cero y varianza σ2, y X= (X1, ..., Xm)0∈Nm(0, σ2Im) el vector aleatorio construido con las observaciones. (i) Si A es una matriz simétrica de orden m×m , idempotente (A2=A) y de rango r≤m , entonces X0AX ∈σ2χ2 r. (ii) Si A es una matriz en las condiciones anteriores y b es un vector de dimensión m , tal que Ab = 0 , entonces X0AX y b0X son independientes. (iii) Si A y B son dos matrices en las condiciones anteriores y AB = 0, entonces X0AX y X0BX son independientes. Teorema 3.1. Supongamos que los errores ε1, ..., εn son independientes y tienen distribución común N(0, σ2) , y X es una matriz n×p de rango p . Entonces (i)ˆ β∈Np(β, σ2(X0X)−1) 52 CAPÍTULO 3. REGRESIÓN CON DATOS SESGADOS POR LONGITUD (ii)RSS σ2=(n−p)ˆσ2 σ2∈χ2 n−p (iii)ˆ β y RSS (ó ˆσ2 ) son independientes. Demostración. Para obtener (i) simplemente tenemos que aplicar la Propiedad 1 (sobre transformaciones lineales de vectores aleatorios normales), teniendo en cuenta que ˆ β= (X0X)−1X0Y , e Y∈Nn(Xβ, σ2I). Por lo tanto ˆ β tendrá distribución normal, cuyo vector de medias y matriz de covarianzas han sido calculados anteriormente. Para obtener (ii) tenemos que emplear el apartado (i) de la Propiedad 2 (sobre transformaciones cuadráticas aplicadas a una muestra de variables normales), pues según hemos visto RSS =ε0(In−H)ε , y como ε= (ε1, ..., εn)0 está en las condiciones de la Propiedad 2, y la matriz (In−H) es simétrica, idempotente y de rango (n−p) , entonces RSS ∈σ2χ2 n−p. Finalmente el apartado (iii) se obtiene también de la Propiedad 2, en este caso de su apartado (ii) , para lo cual basta con observar que ˆ β= (X0X)−1X0Y= (X0X)−1X0Xβ + (X0X)−1X0ε=β+ (X0X)−1X0ε, mientras que RSS =ε0(In−H)ε . Como (X0X)−1X0(In−H) = 0, entonces ε0(In−H)ε y (X0X)−1X0ε son independientes, y por tanto, también lo son ˆ β y RSS . 3.2. Modelo lineal general con datos sesgados por longitud En esta sección vamos a seguir tratando con el modelo lineal general, con la diferencia de que ahora consideraremos que la variable respuesta Y está sesgada por longitud. Daremos una estimación de los parámetros del modelo, que van a ser diferentes a los vistos en la sección anterior debido a la presencia del sesgo, y mencionaremos algunas de sus propiedades. Hemos visto que un modelo de regresión lineal, homocedástico, con errores normales e independientes, del que extraemos una muestra bajo diseño aleatorio nos proporciona datos del tipo (X1, Y1), ..., (Xn, Yn) . Seguimos suponiendo el modelo lineal para las variables originales, pero ahora suponemos que el proceso de observación está sometido a cierto sesgo. De este modo, la muestra resultante será del tipo (Xw 1, Y w 1), ..., (Xw n, Y w n), siendo w una función que indica la probabilidad de observar cada dato bajo el mecanismo de sesgo. 3.2. MODELO LINEAL GENERAL CON DATOS SESGADOS POR LONGITUD 53 Por ejemplo, en el caso de sesgo por longitud, w es proporcional al valor de Y . Por lo tanto, bajo sesgo por longitud, nuestra muestra (Xw 1, Y w 1), ..., (Xw n, Y w n), es una muestra i.i.d. de una variable aleatoria con distribución Fw y cuya densidad viene dada por: dFw(x, y) = y µY dF(x, y). siendo F(x, y) la distribución bivariada de (X, Y ) y µY=E(Y) . 3.2.1. Estimación de los parámetros del modelo Uno de los procedimientos más utilizados en la literatura para obtener una estimación de β , es el método de mínimos cuadrados ponderados, que consiste en: m´ın β n X i=1 1 wi (Yw i−xw iβ)2, (3.1) siendo wi=ω(xω i, Y ω i) . La introducción de los recíprocos de ωi , como pesos en los mínimos cuadrados ponderados, es la manera de corregir el sesgo de los datos para obtener un estimador consistente de los coecientes del modelo lineal en la distribución original de (X, Y ) . En concreto, en el caso de sesgo por longitud, basta con tomar el recíproco de las respuestas (ver Cristóbal y Alcalá (2000)), obteniendo así el siguiente problema de optimización: m´ın β n X i=1 1 Yw i (Yw i−xw iβ)2. La solución al problema anterior es ˆ β=(Xw)0WXw−1(Xw)0WY w, donde Yw es un vector columna con las observaciones Yw i , Xw es una matriz n×p dada por Xw=    xw 1 . . . xw n     =    xw 1,1··· xw 1,p . . . ··· . . . xw n,1··· xw n,p     y W es la siguiente matriz diagonal W=    (Yw 1)−1··· 0 . . ..... . . 0··· (Yw n)−1     . 54 CAPÍTULO 3. REGRESIÓN CON DATOS SESGADOS POR LONGITUD 3.2.2. Propiedades de los estimadores Aunque el estimador por mínimos cuadrados ponderados, ˆ β , es fácil de calcular, sus propiedades de sesgo y varianza no son sencillas de obtener, pues involucra unos pesos aleatorios, (Yw i)−1 en el caso de sesgo por longitud en Y . Aquí vamos a proporcionar las propiedades asintóticas que obtuvieron Ojeda Cabrera y Van Keilegom (2009) para el estimador ˆ β con función de sesgo arbitraria ω(x, y) y en un contexto de regresión paramétrica no necesariamente lineal. Así, se considera un modelo de regresión del tipo Y=m(X, β) + σ(X)ε, donde Y es la variable respuesta, X la variable explicativa, ε el error, que se supone independiente de X (cumple que E(ε)=0 , V ar(ε)=1 ), m(x, β) es la función de regresión , siendo β un parámetro desconocido que debemos estimar, y σ2(x) = V ar(Y/X =x) es la función de varianza condicional. Claramente, el modelo de regresión puede no ser lineal. Cabría, por ejemplo, un modelo periódico como éste m(X, β) = β0+β1·sen 2πX β2, donde la función de regresión no es lineal respecto del vector de parámetros β= (β0, β1, β2) . A continuación enunciamos el resultado sobre la distribución límite del estimador β , que obtuvieron Ojeda Cabrera y Van Keilegom (2009), página 2839. En su Lemma 3.1 se considera una función de peso ω(x, y) arbitraria. Lema 3.2. Consideremos el modelo de regresión Y=m(X, β) + σ(X)ε, del cual se observa una muestra bajo sesgo de selección. Entonces, bajo ciertas condiciones de regularidad en m y σ , el estimador por mínimos cuadrados ponderados, β , presenta la siguiente distribución límite cuando n tiende a innito: √n(ˆ β−β)d −→ N(0,Σ), siendo Ω = E  ∂mβ(X) ∂β ∂mβ(X)T ∂β   y Σ = µwΩ−1E  (Y−∂mβ(X))2 w(X, Y ) ∂mβ(X) ∂β ∂mβ(X)T ∂β  Ω−1. 3.3. SIMULACIONES 55 En el Lema 3.2 se denota µw=E(ω(X, Y )) . Observamos que la matriz Ω juega el papel de la matriz X0X de los modelos lineales. De hecho, en la matriz de covarianzas aparece Ω junto con una matriz ponderada en función del cuadrado del erro y la función de ponderación por sesgo de selección, ω . 3.3. Simulaciones En esta sección se presenta un estudio de simulación en el que se han estudiado las propiedades de sesgo y varianza del estimador por mínimos cuadrados ponderados ˆ β , y se han comparado con el estimador por mínimos cuadrados ordinarios (sin ponderación). Hemos considerado el modelo de regresión lineal simple Y=β0+β1X+ε, donde β0=β1= 3 , X tiene distribución uniforme en el intervalo [0,1] y ε∈χ2 2−2, esto es, el error tiene distribución ji-cuadrado con dos grados de libertad, a la cual le hemos restado 2 para que tenga media cero. Generar muestras de este modelo, añadiendo sesgo por longitud en Y , no es tan sencillo como en el tema anterior, en el cual no había regresión, y había muchos casos de distribuciones cuya versión sesgada era conocida o fácilmente deducible. En el contexto de regresión que nos ocupa, la versión sesgada del vector (X, Y ) no tiene distribución conocida. En estas circunstancias, hemos optado por generar una cantidad ingente de datos del modelo original (un millón en este caso) y tomar muestras de esta pseudo-población con probabilidades proporcionales a los valores de Y . De este modo, resulta una muestra (Xw 1, Y w 1),...,(Xw n, Y w n) obtenida bajo sesgo por longitud de Y . Realmente no hemos tomado una muestra sino mil, para obtener así mil réplicas del estimador por mínimos cuadrados ponderados, ˆ β , cada una calculada sobre una muestra. A partir de estas mil réplicas podremos aproximar el sesgo y la varianza del estimador. En el Cuadro 3.1 se muestran el sesgo y la varianza aproximados en mil simulaciones del modelo, para los estimadores de los coecientes de regresión, por mínimos cuadrados ponderados (bajo el título Con ponderación en el cuadro) y por mínimos cuadrados ordinarios (bajo el título Sin ponderación en el cuadro). Se han considerado dos tamaños de muestra, n= 50 y n= 100 . 56 CAPÍTULO 3. REGRESIÓN CON DATOS SESGADOS POR LONGITUD Respecto del sesgo, observamos que el estimador sin ponderación presenta un sesgo considerable (tanto para la ordenada en el origen como para la pendiente), bastante mayor que el sesgo del estimador con ponderación, y lo que es más grave, ese sesgo no se reduce al aumentar el tamaño muestral. Esto es coherente con la previsible falta de consistencia del estimador sin ponderación, pues no corrige el sesgo de observación, y por tanto, no converge a los verdaderos valores de los coecientes. Por el contrario, el estimador con ponderación reduce su sesgo con el tamaño muestral, y como era previsible, converge a los valores verdaderos de los coecientes. Sin embargo, debemos destacar que es un estimador sesgado, particularmente para tamaños de muestra pequeños o moderados, pues la corrección del sesgo de observación no es perfecta. En relación con la varianza, tanto el estimador con ponderación como el estimador ordinario presentan varianzas decrecientes con el tamaño muestral. Aun así, el estimador con ponderación también muestra varianzas menores que el estimador ordinario. En consecuencia, hemos obtenido las propiedades que cabía esperar de los estimadores, con ponderación por sesgo y sin ella. Podemos concluir que la ponderación por sesgo es necesaria para obtener estimaciones correctas de los coecientes de regresión, como ya ocurría con la media armónica respecto de la media usual en capítulos anteriores. Método Coecientes Sesgo Varianza n=50 n=100 n=50 n=100 Sin ponderación Ordenada en el origen 1.2135 1.2131 0.6697 0.3462 Pendiente -0.5839 -0.5909 1.7435 0.8760 Con ponderación Ordenada en el origen 0.0779 0.02024 0.3309 0.1713 Pendiente -0.0767 -0.0102 0.8806 0.4463 Cuadro 3.1: Sesgo y varianza de los estimadores de los coecientes de regresión, con ponderación y sin ponderación. Apéndice A Código R para las simulaciones A.1. Distribuciones continuas Modelo: F es uniforme en [a,b] #- - - Generación de los datos set.seed(123456) a=9 b=10 media=(a+b)/2 varianza=(b-a)2/12 n=20 # Tamaño muestral ns=1000 # Mil muestras simuladas vm=c() # Vector que recogerá las mil medias aritmética vma=c() # Vector que recogerá las mil medias armónica #- - - Generamos muestras for (is in 1:ns){ u=runif(n,0,1) x=sqrt(a2+u*(b2-a2)); x # Muestra sesgada # Media aritmética simple vm[is]=mean(x) # Media armónica vma[is]=1/mean(1/x) } 57 64 APÉNDICE A. CÓDIGO R PARA LAS SIMULACIONES } #- - - Propiedades de la media aritmética simple sesgo=mean(vm)-lambda et_hat=sd(vm); et_hat # Error típico aproximado por simulación #- - - Propiedades de la media armónica sesgo=mean(vma)-lambda et_hat=sd(vma); et_hat # Error típico aproximado por simulación var=(et_hat)2 ECM=var+sesgo2 sigma2/n new=(pi 2/6)*lambda ((media2)*(media*new-1))/n # Varianza asintótica #- - - Representación gráca par(mfrow=c(1,2)) #- - Poisson(2) no sesgada lambda=2 .x <- 0:30 plot(.x, dpois(.x, lambda=lambda,log = FALSE),col=2,xlim=c(0,10),ylim=c(0,1.5), xlab="x", ylab="Probabilidad", main="Poisson(2)", type="h") points(.x, dpois(.x, lambda=lambda), col=2,pch=16) abline(h=0, col="gray") #- - Poisson(2) sesgada plot(.x, 1+dpois(.x, lambda=lambda,log = FALSE),col=3,xlim=c(0,10),ylim=c(0,1.5), xlab= 2 ", ylab="Probabilidad", main="1+Poisson(2)", type="h") points(.x, 1+dpois(.x, lambda=lambda), col=3,pch=16) abline(h=0, col="gray") A.3. REGRESIÓN CON DATOS SESGADOS POR LONGITUD 65 A.3. Regresión con datos sesgados por longitud set.seed(123456) nn=1000000 # Tamaño poblacional de referencia beta=c(3,3) # Coecientes de regresión xg=runif(nn) eps=rchisq(nn,df=2)-2 yg=beta[1]+beta[2]*xg+eps #plot(xg,yg,ylim=c(0,max(y))) #abline(a=2,b=3) p=yg/sum(yg) n=50 # Generamos muestras ns=1000 # Mil muestras simuladas m_sin=matrix(0,nrow=ns,ncol=2) # Matriz que recogerá los mil pares de coecientes (sin ponderación) m_con=matrix(0,nrow=ns,ncol=2) # Matriz que recogerá los mil pares de coecientes (con ponderación) for (is in 1:ns){ #- - - Generación de la muestra sesgada ind=sample(nn,n,prob=p) x=xg[ind] y=yg[ind] # plot(x,y,ylim=c(0,max(y))) #- - - Ajuste lineal sin ponderación beta_sin=coef(lm(yx)) m_sin[is,]=beta_sin #- - - Ajuste lineal con ponderación beta_con=coef(lm(yx,weights=1/y)) 66 APÉNDICE A. CÓDIGO R PARA LAS SIMULACIONES m_con[is,]=beta_con if (is==oor(is/5)*5){ cat(Sample,is,Coef_sin",beta_sin,Coef_con,beta_con,\n)} } # Propiedades de los coecientes sin ponderación sesgo_sin=colMeans(m_sin)-beta; sesgo_sin et_sin=c(sd(m_sin[,1]),sd(m_sin[,2])); et_sin # Error típico aproximado por simulación # Propiedades de los coecientes con ponderación sesgo_con=colMeans(m_con)-beta; sesgo_con et_con=c(sd(m_con[,1]),sd(m_con[,2])); et_con # Error típico aproximado por simulación Bibliografía [1] Feller,W. (1971), An Introduction to Probability Theory and Its Applications. Vol II, John Wiley & Sons Inc, New York. [2] Patil, G. P. (1984). Studies in statistical ecology involving weighted distributions. Statistics: Applications and New Directions , Indian Stat. Inst., Calcutta, 478-503. [3] Cox, D.R.. (1969). Some Sampling Problems in Technology. In: New Developments in Survey Sampling , John Wiley, New York, 506-527. [4] Vardi, Y. (1982). Nonparametric estimation of a regression function from recurrence times. Ann. Stat. , 10 , 616-620. [5] Patil, G. P. (1978). Weighted Distributions and Size-Biased Sampling with Applications to Wildlife Populations and Human Families. Biometrics , 34 , 179-189. [6] Vélez, R. y A. García (1993) Principios de Inferencia Estadística. , UNED. [7] Jorge L. Ojeda Cabrera and Ingrid Van Keilegom (2009).Goodness-of-t tests for parametric regression with selection biased data Journal of Statistical Planning and Inference Vol 139, , 8 , 2836 - 2850. [8] Cristóbal, J. A., Ojeda, J. L., Alcalá, J. T (2004).. Condence bands in nonparametric regression with length biased data. Annals of the Institute of Statistical Mathematics , 56 (3) , 175 - 196. 67