Modelos no lineales de regresión

Traballo Fin de Grao Modelos no lineales de regresión Jose Antonio Cribeiro Ramallo 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Modelos no lineales de regresión Jose Antonio Cribeiro Ramallo 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA iii iv Trabajo propuesto Área de Coñecemento: Estítstica e Investigación Operativa Título: Modelos no lineales de regresión Breve descrición do contido Los modelos de regresión estudian el efecto de ciertas variables explicativas sobre una variable de interés, también llamada variable respuesta. Los primeros modelos de regresión son lineales, en el sentido de que los parámetros desconocidos son coecientes que multiplican a las variables explicativas o a funciones de ellas. Sin embargo, hay modelos muy conocidos de la Física, la Química o la Biología, como la ley de enfriamiento de Newton, los modelos de crecimiento o los modelos periódicos, en los cuales no es posible expresar los parámetros como coecientes. En este trabajo se revisarán los métodos de inferencia en modelos no lineales de regresión. Mediante simulaciones, se estudiarán las propiedades de los estimadores bajo condiciones diversas, como normalidad o no de la distribución del error, homocedasticidad o heterocedasticidad, u otras circunstancias. Asimismo, se ilustrará el trabajo aplicando modelos no lineales de regresión en datos reales. Recomendacións Outras observacións Índice general Resumen vii Introducción ix 1. Inferencia en modelos no lineales 1 1.1. Presentación del modelo no lineal . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Estimadores del modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1. Estimador de mínimos cuadrados . . . . . . . . . . . . . . . . . . . . 2 1.2.2. La varianza del error . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3. Comportamiento asintótico de los estimadores . . . . . . . . . . . . . . . . . 3 1.4. Intervalos de conanza y contrastes basados en la distribución asintótica . . 14 1.5. Intervalos de conanza basadas en el perl de RSS . . . . . . . . . . . . . . 14 2. Simulación 17 2.1. Motivación y Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2. Sesgo........................................ 18 2.3. Varianzaasintótica................................ 19 2.4. Intervalos de conanza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5. Contrastes..................................... 24 2.6. Distribuciones................................... 25 2.6.1. Efecto de la distribución del error . . . . . . . . . . . . . . . . . . . . 28 2.7. Conclusión..................................... 29 Código de R empleado 31 Bibliografía 47 v Resumen En este trabajo se ha hecho una revisión de los modelos no lineales. Se han estudiado las propiedades de consistencia y de distribución asintótica de los estimadores, y se han presentado los procedimientos para la inferencia sobre los parámetros del modelo. Además, se ha hecho un estudio por simulación en donde se han vericado varias propiedades teóricas de los estimadores y de los métodos inferenciales expuestos en la teoría, como intervalos de conanza y contrastes de hipótesis, acompañado de una breve conclusión de los resultados obtenidos Abstract In the following paper, we have revised the methods for the non linear regression models. We have studied the properties of consintency and the asymptotic distribution of the model estimator, on top of having presented a range of tools for the parameter inference. Furthermore, a simulation study has been taken where we have veried the multiple theorical properties of the estimators and the inferential methods exposed previosly, like the contrast statistics and condence intervals, followed by a brief conlcusion of the results that have been obtained. vii 2 CAPÍTULO 1. INFERENCIA EN MODELOS NO LINEALES más grande, el estimador converge a su valor verdadero y las tangentes se vuelven cada vez mejor aproximación local. Es por esto que el estimador tendrá un comportamiento cada vez más simétrico, y cada vez más parecido al de un estimador lineal dado por la tangente, tendiendo así a un comportamiento normal. El propósito será estimar el vector de parámetros θ∈Rp y la varianza del error σ2 a partir de una muestra, de las variables explicativas y la respuesta, de tamaño n . 1.2. Estimadores del modelo 1.2.1. Estimador de mínimos cuadrados Tal y como hacíamos en la regresión lineal, será conveniente tomar el estimador que nos garantice la menor suma residual de cuadrados: RSSn(θ) = n X i=1 (yi−m(xi, θ))2 (1.1) Vamos a denotar a este estimador por ˆ θ , tal que: ˆ θn=argminθ∈ΘRSSn(θ) es decir, ˆ θn sería el θ∈Θ tal que la derivada de (1.1) es cero: n X i=1 (yi−m(xi, θ))∂m ∂θ (xi, θ)=0. (1.2) Observación 1.2 . Si tomamos θ∈Rp , la derivada de (1.1) sería un vector gradiente. Al sistema de p ecuaciones y p incógnitas (1.2) lo vamos a conocer como las ecuaciones normales de regresión (al igual que en el caso lineal.) A diferencia que en el caso lineal, las ecuaciones normales de regresión son, en este caso, no lineales, por lo que no vamos a poder obtener una expresión explícita del estimador. Por lo tanto, la única forma de obtener el estimador de mínimos cuadrados será aproximándolo mediante algún algoritmo iterativo. Observación 1.3 . Recordemos que, en el modelo lineal, que la matriz X0X tuviera autovalores pequeños generaba problemas al invertir la matriz, suponiendo un incremento en la varianza de los estimadores (lo llamábamos problema de colinealidad). En el caso no lineal volvemos a tener este problema, y es que la inversión de la matriz hessiana ∂2RSS ∂2θ (simétrica y semidenida positiva) juega un papel crucial en algunos procesos 1.3. COMPORTAMIENTO ASINTÓTICO DE LOS ESTIMADORES 3 iterativos (como el método de Newton). Por lo que si ∂2RSS ∂2θ presenta autovalores pequeños, esta inversión se diculta; llegando a, incluso, hacer que el algoritmo no converja. El problema de la convergencia se puede resolver seleccionando correctamente los valores iniciales, utilizando métodos iterativos que tratan de paliar este problema (como el de Levenberg-Marquardt) o incluso tratar con otros métodos en donde no se emplean las derivadas. Las consecuencias sobre la varianza de los estimadores solo se pueden paliar modicando el modelo. 1.2.2. La varianza del error De forma idéntica al modelo lineal, nos encontramos con que ˆσ2= RSS (ˆ θn) n−p=Pn i=1(yi−m(xi,ˆ θn))2 n−p, es un estimador de la varianza del error que asintoticamente es una ji-cuadrado de n−p grados de libertad (lo veremos más adelante). 1.3. Comportamiento asintótico de los estimadores Es importante destacar que no hemos supuesto al principio que los errores del modelo tengan que ser normales para demostrar la normalidad asintótica de los estimadores en el modelo no lineal. Esto no es algo único del modelo no lineal, si no que también ocurre en la regresión lineal. Eicker (1963) y Grenander y Rosenblatt (1957) dieron condiciones sobre la función de regresión para la normalidad asintótica cuando los errores no son idénticamente distribuidos y condiciones para la eciencia asintótica. Nuestro objetivo va a ser demostrar la normalidad asintótica de los estimadores de mínimos cuadrados bajo ciertas hipótesis. La primera de ellas será: Hipótesis 1.4. yi=m(xi, θ0) + εi=fi(θ0) + εi, i = 1,2, . . . ; con fi funciones continuas conocidas en el subconjunto compacto Θ de un espacio Euclidiano y εi errores iid de media 0 y varianza σ2>0 ( θ0 y σ son desconocidos). Observación 1.5 . Durante el resto de esta sección vamos a considerar Θ como en la hipótesis 1.4, un subconjunto compacto de un espacio Euclidiano. 4 CAPÍTULO 1. INFERENCIA EN MODELOS NO LINEALES Observación 1.6 . Recordemos que estamos trabajando con una única variable explicativa y una única respuesta. Gallant[2] habla sobre el caso multivariante en su artículo, en donde establece una analogía entre el modelo no lineal multivariante asintótico con el lineal multivariante. Vamos a llamar a cualquier vector ˆ θ∈Θ que minimice Qn(θ) = 1 n n X i=1 (yi−fi(θ))2 como estimador de mínimos cuadrados de θ0 basado en los n primeros valores de y . Nuestro objetivo es obtener la distribución límite del estimador de los parámetros ˆ θ , para ello seguirémos la demostración dada por Jenrich [4]. Empezamos considerando una serie de conceptos que vamos a necesitar a lo largo de la demostración. Denición 1.7. Sean x= (xi) e y= (yi) dos sucesiones de números reales, y sea hx, yin=n−1 n X i=1 xiyi el producto interior n -ésimo de x e y . Si hx, yin converge cuando n−→ ∞ a un número hx, yi , vamos a conocer a este límite como el producto de colas ( tail product ) de x e y . Denición 1.8. Sean g= (gi), h = (hi) sucesiones de funciones de Θ , tal que g(θ) = (gi(θ)), h(θ) = (hi(θ)) (es decir, g y h las podemos considerar como funciones de Θ al espacio de sucesiones de elementos de Θ , ΘN ). Consideremos que hg(α), h(β)in−→ hg(α), h(β)i uniformemente ∀α, β ∈Θ . Llamaremos producto cruzado de colas ( tail cross product ) de g y h a la función: [g, h]:(α, β)∈Θ×Θ7→ hg(α), h(β)i ∈ R Nótese que si g y h son funciones continuas, entonces el producto cruzado de colas de ambas es una función continua. Observación 1.9 . Denotaremos a hx, xin como: kxk2 n=hx, xin y a su límite como: kxk2=hx, xi, 1.3. COMPORTAMIENTO ASINTÓTICO DE LOS ESTIMADORES 5 que lo llamaremos la norma de cola ( tail norm ) de x . Notese que si [g, g],[h, h], y [g, h] existen, entonces: kg(α)−h(β)kn−→ kg(α)−h(β)k uniformemente para todo α, β ∈Θ . Observación 1.10 . Se deduce facilmente que ka+bk2 n=kak2 n+ 2 ha, bin+kbk2 n Consideremos ahora dos teoremas que serán útiles a la hora de probar la existencia de [g, h] para dos funciones acotadas y continuas en X×Θ , con X un espacio Euclidiano. Teorema 1.11. Sean X un espaio Euclidiano y g una función continua y acotada en X×Θ . Sea (Fn) una sucesión de distribuciones en X que convergen a una distribución F Entonces: Zg(x, θ)dFn(x)−→ Zg(x, θ)dF (x). Teorema 1.12. Sea g una función en X×Θ con X un espacio Euclidiano. Sea g(x, θ) una función continua de θ para cada x y una función medible de x para cada θ , tal que: |g(x, θ)| ≤ h(x); para todo x, θ. con h integrable respecto a F , una distribución en X . Si x1, x2, . . . son una muestra aleatoria de F , entonces, para casi toda sucesión (xn) de muestras: n−1 n X i=1 g(xi, θ)−→ Zg(x, θ)dF (x) uniformemente, para todo θ∈Θ . Observación 1.13 . La demostración de estos teoremas se encuentra en Jenrich [4, p. 635- 636] De las deniciones anteriores, extraemos un resultado que vamos a necesitar más adelante: Sea ε= (εi) , con εi bajo la hipótesis 1.4. Teorema 1.14. Si existe el producto de cola de una sucesion x de números reales se verica que: hx, εin−→ 0, para casi todo ε . Este resultado se puede extender como: 6 CAPÍTULO 1. INFERENCIA EN MODELOS NO LINEALES Teorema 1.15. Bajo la hipótesis 1.4, si g es una sucesión de funciones de Θ y si [g, g] existe, tenemos que, para casi todo ε : hg(θ), εin−→ 0 uniformemente para todo θ∈Θ . Demostración. Como g es continua y kg(α)−g(β)kn converge uniformemente en α ; existe, para todo δ > 0 y para todo β∈Θ , un entorno N de β tal que: kg(α)−g(β)kn< δ, para todo α∈N y para todo n∈N sucientemente grande. Luego, usando el Teorema anterior en la inecuación: |hg(α), εin| ≤ kg(α)−g(β)knkεkn+|hg(β), εin|, obtenida a partir de la inecuación de Cauchy, deducimos que para casi todo ε , Θ es recubierto por entornos N tales que |hg(θ), εin|< δ para todo θ∈N , cuando n∈N es lo sucientemente grande. De esta forma, como Θ es compacto, hay una colección nita de estos entornos que recubren a N, por lo que para casi todo ε : |hg(θ), εin|< δ, para todo θ∈Θ con n lo sucientemente grande. Es decir, |hg(θ), εin−0|< δ =⇒ hg(θ), εin−→ 0. Vamos a establecer ahora una forma del Teorema Central del Límite que emplea las nuevas deniciones que hemos introducido y nos será de utilidad más adelante: Teorema 1.16. Sea ε según la hipótesis 1.4 y supongamos que la norma de cola de x , una sucesión de números reales, existe. Entonces: L√nhx, εin−→ N0,kxk2σ2 Este resultado puede ser extendido al caso multilineal como: 1.3. COMPORTAMIENTO ASINTÓTICO DE LOS ESTIMADORES 7 Corolario 1.17. Sea ε bajo las mismas hipótesis que en el Teorema anterior. Si existe cada producto de cola posible entre las sucesiones x1,...xp de números reales, entonces, la distribución conjunta de: √nhx1, εin,...,hxp, εin converge a N0, σ2X, siendo X la matriz tal que el elemento (i, j) -ésimo es hxi, xji . El siguiente resultado es de vital importancia, ya que nos garantiza que la hipótesis 1.4 es suciente para demostrar la existencia de estimadores de mínimos cuadrados medibles en el caso no lineal. Lema 1.18. Sea Q una función real en Θ×Y , con Y un espacio medible. Para cada θ∈Θ , sea Q(θ, y) una función medible de y , y para cada y∈Y una función continua de θ . Entonces, existe ˆ θ:Y−→ Θ medible, tal que, para todo y∈Y : Qˆ θ(y), y= ´ınf θ∈ΘQ(θ, y) Demostración. Sea (Θn) una sucesión creciente de subconjuntos nitos de Θ , cuyo límite es denso en Θ . Por el Th. de Weierstrass, para cada n hay una función medible, ˜ θn:Y−→ Θn tal que: Q˜ θn(y), y= ´ınf θ∈Θn Q(θ, y) Sea ˜ θn1 la primera componente de ˜ θn y ˆ θ1= l´ım n→∞ ˜ θn su límite (nótese que ˆ θ1 es medible por ser límite de medibles). Para cada y∈Y hay una subsucesión ˜ θni(y) de ˜ θn(y) , que converge a un punto ¯ θ=ˆ θ1(y),¯ θ2,...,¯ θp. 8 CAPÍTULO 1. INFERENCIA EN MODELOS NO LINEALES Ahora, como l´ım n→∞ Θn es denso en Θ : ´ınf (θ1,...,θn)∈ΘQˆ θ1(y), θ2, . . . , θp, y≤Q¯ θ, y (1.3) = l´ım i→∞ Q˜ θni(y), y (1.4) = l´ım i→∞ ´ınf θ∈Θni Q(θ, y) (1.5) = ´ınf θ∈ΘQ(θ, y). (1.6) Por lo tanto: ´ınf (θ1,...,θn)∈ΘQˆ θ1(y), θ2, . . . , θp, y= ´ınf θ∈ΘQ(θ, y) ; ∀y∈Y. Sea ahora, Q0((θ1, . . . , θp), y) = Qˆ θ1(y), θ2, . . . , θp, y, luego, Q0(θ, y) es una función continua de θ∀y∈Y , y una función medible de y∀θ∈Θ. Con el mismo razonamiento que hicimos para Q , encontramos que ´ınf (θ1,...,θn)∈ΘQˆ θ1(y),ˆ θ2(y), θ3, . . . , θp, y= ´ınf θ∈ΘQ(θ, y) ; ∀y∈Y. De esta forma, repitiendo esto encontramos las funciones medibles de Y a R : ˆ θ1(y),...,ˆ θp(y) ; tales que: Qˆ θ1(y),...,ˆ θp(y), y= ´ınf θ∈ΘQ(θ, y) ; ∀y∈Y. El Lema 1.18 junto con el Th. del valor medio nos otorga el siguiente resultado, que no es más que una versión del Th. del valor medio para funciones medibles. Proposición 1.19. Sea f:Y×Θ−→ R , con Y un espacio medible y Θ convexo. Para cada θ∈Θ , f(y, θ) es una función medible de y y para cada y∈Y una función continuamente diferenciable de θ . Sean θ1, θ2 funciones medibles de Y a Θ . Entonces, existe una función ¯ θ:Y−→ Θ medible, tal que: i)f(y, θ1(y)) −f(y, θ2(y)) = ∂ ∂θ fy, ¯ θ(y)(θ1(y)−θ2(y)) . ii)ˆ θ(y) está en el segmento que une θ1(y) con θ2(y) , para todo y∈Y . Demostración. Sea D(y, θ) la distania de θ al segmento que una θ1(y) y θ2(2) , y sea P(θ, y) = |f(y, θ1(y)) −f(y, θ2(y)) −∂ ∂θf(y, θ) (θ1(y)−θ2(y)) |+D(y, θ). 1.3. COMPORTAMIENTO ASINTÓTICO DE LOS ESTIMADORES 9 Si vemos a θ1 y θ2 como los límites puntuales de una sucesión de funciones simples y a ∂ ∂θ f(y, θ) como: ∂ ∂θf(y, θ) = l´ım h→0 f(y, θ +h)−f(y, θ) h, es claro que P es una función medible de y para todo θ∈Θ , y una función continua de θ para todo y∈Y . Por esto último, en virtud del Lema 1.18, encontramos que ¯ θ:Y→Θ medible, tal que ¯ θ(y) minimiza a P(θ, y) para todo y∈Y como función de θ . Además, por el Th de valor de medio, cuando θ=¯ θ(y) el primer y el segundo término de P(θ, y) son 0, luego: i)|f(y, θ1(y)) −f(y, θ2(y)) −∂ ∂θ fy, ¯ θ(y)(θ1(y)−θ2(y)) |= 0 =⇒f(y, θ1(y)) − f(y, θ2(y)) = ∂ ∂θ fy, ¯ θ(y)(θ1(y)−θ2(y)) ii)Dy, ¯ θ(y)=0=⇒¯ θ(y)∈L(θ1(y), θ2(y)) De ahora en adelante vamos a necesitar considerar una nueva condición sobre las fi de la hipótesis 1.4: Hipótesis 1.20. • ∃[f, f], con f= (fi) : Θ −→ ΘN • La función Q(θ) = kf(θ)−f(θ0)k2 tiene un único mínimo en θ=θ0 Observación 1.21 . Las hipótesis 1.4 y 1.20 son sucientes para probar la continuidad de Q(θ) Observación 1.22 . Notese que la Qn que denimos al principio coincide con Qn(θ) = kf(θ0) + ε−f(θ)k2 n Esta hipótesis es necesaria para demostrar una propiedad muy importante de los estimadores: la consistencia. Teorema 1.23. Sea ˆ θn una sucesión de EMC. Bajo las hipótesis 1.4 y 1.20, ˆ θn y ˆσ2 n=Qnˆ θn son estimadores fuertemente consistentes. Demostración. Consideremos primero que kεk=σ y que: hf(θ)−f(θ0), εin−→ 0 10 CAPÍTULO 1. INFERENCIA EN MODELOS NO LINEALES uniformemente, para todo θ∈Θ , gracias a la Ley fuerte de los grandes números y al Th. 1.15 respectivamente. Qn(θ) se puede escribir como: Qn(θ) = kf(θ0) + ε−f(θ)k2 n =kf(θ0)−f(θ)k2 n+ 2 hf(θ0)−f(θ), εin+kεk2 n, por lo tanto, encontramos que: Qn(θ)−→ kf(θ0)−f(θ)k2+θ2=Q(θ) + σ2 , uniformemente para todo θ∈Θ. Sean θn=ˆ θn(f(θ0) + ε) , θ0 el límite de la sucesión (θn) y (θni) una subsucesión cualquiera que converge a θ0 . Por la continuidad de Q y la convergencia uniforme de Qn a Q+σ2 , Qni(θni)−−−−−−−→ i→∞ Qθ0+σ2 Como σni es un EMC, Qni(θni)≤Qni(θ0) := kεk2 ni. Si i−→ ∞ , encontramos que: Qθ0+σ2≤σ2=⇒Qθ0= 0 ⇐⇒ θ0=θ0 (esto lo deducimos ya que el mínimo de Q es único y Q(θ)≥0 ) . En conclusión, lo que tenemos es que: θn−→ θ0, para casi todo ε=⇒ˆ θcs −−→ θ0 y ˆσ2 n=Qnˆσ2 ncs −−→ Q(θ0) + σ2=σ2. Observación 1.24 . Consideremos la notación: f0 ij(θ) = ∂ ∂θj fi(θ); f00 ijk (θ) = ∂2 ∂θj∂θ0 k fi(θ), de tal forma que: f0 j= (f0 ij); f00 jk = (f00 ijk ). Además: [an(θ)]j1,j2=f0 j1(θ), f0 j2(θ)n, 1.3. COMPORTAMIENTO ASINTÓTICO DE LOS ESTIMADORES 11 tal que [an(θ)]j1,j2 es la componente (j1, j2) -ésima de la matriz an(θ) : an(θ) =     [an(θ)]1,1··· [an(θ)]1,p . . .. . . [an(θ)]p,1··· [an(θ)]p,p     Consideremos también: [a(θ)]j1,j2=f0 j1(θ), f0 j2(θ) en donde [a(θ)]j1,j2 es la componente (j1, j2) -ésima de la matriz a(θ) : a(θ) =     [a(θ)]1,1··· [a(θ)]1,p . . .. . . [a(θ)]p,1··· [a(θ)]p,p     , y el producto de colas de f0 j1(θ) y f0 j2(θ) . Vamos a considerar ahora las dos últimas hipótesis que necesitamos para poder demostrar la normalidad asintótica de los estimadores. Sea f0 j1= (f0 ij1) y f00 j1j2= (f00 ij1j2) sucesiones de funciones, tales que: f0 j1, f00 j1j2: Θ −→ ΘN. Hipótesis 1.25. f0 j1, f00 j1,j2 existen, son continuas en Θ y todos los productos cruzados de colas posibles de f, f0 j1, f00 j1j2 existen. Hipótesis 1.26. El vector de los verdaderos parámetros, θ0 , es un punto interior de Θ y a(θ0) es no-singular. Teorema 1.27. Sea ˆ θn una sucesión de EMC de θ0 . Bajo las hipótesis 1.4, 1.20, 1.25 y 1.26; L√nˆ θn−θ0−→ N0, σ2a−1(θ). Además, anˆ θn es un estimador fuertemente consistente de a(θ0) . Demostración. Ya que ˆ θn−→ θ0 en casi todo punto, hay una sucesión (˜ θn) que es tail equivalent a (ˆ θn) (es decir, k(˜ θn)k2=k(ˆ θn)k2 ) tal que cada ˜ θn toma valores en un entorno convexo compacto de θ0 en ˚ Θ . 18 CAPÍTULO 2. SIMULACIÓN 2.2. Sesgo Hemos visto que ˆ θ es un estimador fuertemente consistente de θ , que además es asintóticamente insesgado, es decir, E(ˆ θn) converge a θ0 cuando n−→ ∞ . Con el objetivo de estudiar experimentalmente esta convergencia tenemos que considerar que el valor de E(ˆ θn) ha de ser aproximado, ya que carecemos de una expresión analítica. Por este motivo, vamos a simular S muestras del modelo (2.1) con σ= 5 , K= 20 , V= 120 ; de tal forma que vamos a tener S estimadores obtenidos cada uno de ellos con una de las S muestras:     (x1 1, y1 1). . . (x1 n, y1 n) . . .. . . (xs 1, ys 1). . . (xs n, ys n)     ←− ˆ θn1=ˆ K1,ˆ V1 . . . ←− ˆ θnS=ˆ KnS,ˆ VnS , en donde hemos tomado x1, . . . , xn de manera equiespaciada en [10,110] , yi=m(xi, θ)+εi para todo i , y con ε1, . . . , εn normales de media cero y desviación típica σ . Estos S estimadores simulados conforman una muestra de variables aleatorias independientes e idénticamente distribuidas ( ˆ θ1 d =. . . d =ˆ θs d =ˆ θ ) a la que llamaremos T= (K,V) . De esta muestra simulada conocemos que, dado n el tamaño muestral y ˆ θi el estimador obtenido en la i -ésima muestra simulada, en virtud del Teorema central del límite: 1 s s X i=1 ˆ θi=¯ T −−−−−−−→ N E(ˆ θn),V(ˆ θn) S!, Cuando S tiende a innito. Luego, tomando la media muestral de todos los ˆ θi simulados (llamemosle ¯ T=¯ K,¯ V ), para un S lo sucientemente grande como para minimizar los efectos de la varianza, obtendríamos una buena aproximación de E(ˆ θn) . En las grácas 2.1 podemos apreciar como el sesgo disminuye cada vez más a medida de que aumenta el tamaño muestral, llegando a estar cerca del cero. Observación 2.1 . En estas mismas grácas (las 2.1) podemos ver cómo a partir de n= 100 , las dos parecen sugerir pequeños aumentos y disminuciones en el sesgo para tamaños muestrales mayores. Esto se debe a la desviación típica de ¯ K y de ¯ V en relación al tamaño del sesgo. Por ejemplo, en n= 100 la desviación típica de ˆ K sería aproximadamente q1,35 S= 0,0082 , mientras que el sesgo vale 0,02 . Teniendo en cuenta que estamos trabajando con sesgos tan pequeños, es normal que tengamos pequeños errores de aproximación 2.3. VARIANZA ASINTÓTICA 19 (a) Función del sesgo de ˆ K (b) Función del sesgo de ˆ V Gráca 2.1: Sesgo de ˆ K y ˆ V obtenida a partir de 20.000 simulaciones 2.3. Varianza asintótica Tal y como hicimos con el sesgo, necesitamos encontrar aproximaciones de la varianza de ˆ Kn y ˆ Vn para poder comprobar mediante simulación que asintóticamente se aproximan a los valores teóricos de V(K) y V(V) respectivamente. Dados los estimadores simulados T , conocemos que Sc(K) = 1 S−1Ps i=1 ˆ Ki−¯ K2 , la cuasivarianza de los valores de ˆ Kn s, es una aproximación por Montecarlo de la varianza verdadera de ˆ Kn . Por lo tanto, con un S lo sucientemente grande como para minimizar el error de simulación, deberíamos obtener que la cuasivarianza se aproxima a la varianza verdadera de ˆ Kn , V ar(ˆ Kn) ). Por otro lado, tenemos la varianza de la distribución asintótica de ˆ Kn , que denotaremos AV ar(ˆ K) y que por el Teorema 1.27 vale AV ar(ˆ K) = V(K) = σ2 na−1(θ) . Ahora bien, como no conocemos la expresión de a(θ) , vamos a aproximarla por an(θ) , de tal forma que para calcular AV ar(ˆ Kn) , vamos a tomar: AV ar(ˆ K) = σ2 na−1 n(K), con σ y K conocidos. En concreto, queremos establecer que la razón de las varianzas converge a 1: Sc(K) AV ar(ˆ Kn)−→ 1, 20 CAPÍTULO 2. SIMULACIÓN cuando n−→ ∞ . De forma similar, establecemos que: Sc(V) AV ar(ˆ V)−→ 1, cuando n−→ ∞ . En las grácas 2.2 podemos observar un comportamiento tal como el esperado para una muestra simulada de tamaño S= 100000 . (a) ˆ Kn (b) ˆ Vn Gráca 2.2: Razón de varianzas muestrales y asintóticas de los estimadores Observación 2.2 . Nótese que, al igual que pasaba con el sesgo, las grácas presentan comportamientos no deseados a partir de n=25. Esto es debido, nuevamente, a que la variación de Sc(K) es relativamente grande en comparación a lo cerca que está de AV ar(ˆ Kn) . En el cuadro 2.1, podemos ver una comparación entre los valores de Sc(K) y AV ar(ˆ Kn) , y el intervalo de conanza para la varianza de K para reejar el error de aproximación. Vemos como, incluso, para valores pequeños de n AV ar(ˆ K) no está en los intervalos. 2.4. Intervalos de conanza A nivel teórico hemos visto dos tipos de intervalos para los parámetros del modelo: uno asintótico, asumiendo normalidad del estimador por tener un n lo sucientemente grande, y otro mediante el perl de RSS . Por lo tanto, sería de interés comprobar mediante simulación la cobertura de estos intervalos para diferentes tamaños muestrales. 2.4. INTERVALOS DE CONFIANZA 21 n Sc(K)AV ar(ˆ K) Intervalo de conanza 5 22.83 20.82 (22.62, 23.02) 10 12.34 11.91 (12.23, 12.44) 15 8.53 8.30 (8.45, 8.60) 20 6.51 6.37 (6.45, 6.56) 25 5.27 5.19 (5.22, 5.32) 30 4.38 4.35 (4.33, 4.41) 35 3.81 3.75 (3.77, 3.84) 40 3.34 3.30 (3.30, 3.36) 45 3.00 2.95 (2.97, 3.02) 50 2.69 2.66 (2.66, 2.71) Cuadro 2.1: Valores de Sc(K) , AV ar(ˆ K) , y el intervalo elemental basado en la X2 S−1 sobre la muestra K ; para diferentes tamaños muestrales Para cada par ˆ Kni,ˆ Vni , calcularemos los intervalos de conanza y de perl, y contaremos en cuántos de ellos están los verdaderos valores de K y V . Con un número de simulaciones S= 20000 y σ= 15 , para intervalos con nivel de conanza del 95 % , vemos que para tamaños muestrales pequeños los intervalos de perl tienen una cobertura más próxima al nivel que los asintóticos, mientras que los intervalos ansitóticos tiene coberturas por debajo del nivel. Para muestras mayores, los intervalos asintóticos comienzan a comportarse mejor (ver cuadro 2.2). n ICA pRSS K V K V 15 0.9390 0.9467 0.9471 0.9537 30 0.9426 0.9473 0.9485 0.9482 Cuadro 2.2: Cobertura de los intervalos de conanza asintóticos y de perl de RSS ( ICA y pRSS ) para K y V , con σ= 15 y 20000 simulaciones Si considerásemos σ= 5 , conseguimos un resultado diferente. Tal y como vemos en el cuadro 2.3, los intervalos asintóticos tienen más cobertura que el nivel nominal del 5 % para n= 15 Vamos a estudiar ahora la longitud de los intervalos. En el cuadro 2.4 observamos que la 22 CAPÍTULO 2. SIMULACIÓN n ICA pRSS K V K V 15 0.9647 0.9660 0.9470 0.949 30 0.9570 0.9576 0.9490 0.948 Cuadro 2.3: Cobertura de los intervalos de conanza asintóticos y de perl de RSS ( ICA y pRSS ) para K y V , con σ= 5 y 20000 simulacione longitud de los intervalos asintóticas es menor que los de perl. Notese que estos intervalos tienen menor cobertura, como veíamos en el cuadro 2.2. n ICA pRSS K V K V 15 12.45 ; 20.270 12.570 ; 20.296 30 8.544 ; 13.990 8.5650 ; 13.900 Cuadro 2.4: Longitud de los intervalos de conanza asintóticos (ICA) y de perl RSS (pRSS) para K y V , con σ= 5 y 20000 simulaciones Otra característica de los intervalos dados por el perl RSS es que son asimétricos. Así pues, dada una muestra única de datos, vamos a representar los intervalos sobre la función perl de los parámetros y veamos como se construyen. Consideremos, entonces: pRSS(K) = m´ın VRSS(K, V ), pRSS(V) = m´ın KRSS(K, V ); las funciones perl de K y de V . Podemos ver en la gráca 2.3 que, como era de esperar, no encontramos una asimetría muy grande en los intervalos para σ= 5 debido a que para desviaciones pequeñas los estimadores de nuestro modelo convergen más rápido a una normal. Si tomamos σ= 15 podemos apreciar mucho mejor la asimetría respecto al estimador, tal y como se ve en la gráca 2.4. Para una mejor visualización, en el cuadro 2.5 se han dado las razones de las asimetrías de los intervalos, es decir: Estimador − Extremo izquierdo del intervalo Extremo derecho del intervalo − Estimador 2.4. INTERVALOS DE CONFIANZA 23 ˆ Knˆ Vn σ= 15, n = 15 2.16 1.66 σ= 15, n = 30 1.54 1.33 σ= 5, n = 15 1.33 1.21 σ= 5, n = 30 1.17 1.11 Cuadro 2.5: Razón de asimetrías de los intervalos de perl. (a) K, n = 15 (b) K, n = 30 (c) V, n = 15 (d) V, n = 30 Gráca 2.3: Función perl de ambos parámetros para n= 15 y n= 30 , con σ= 5 . En rojo el valor del estimador y en azul los intervalos. 24 CAPÍTULO 2. SIMULACIÓN (a) K, n = 15 (b) K, n = 30 (c) V, n = 15 (d) V, n = 30 Gráca 2.4: Función perl de ambos parámetros para n= 15 y n= 30 , con σ= 15 . En rojo el valor del estimador y en azul los intervalos. 2.5. Contrastes Planteamos los siguientes contrastes unilaterales sobre el valor de los parámetros: (H0:K= 20; H1:K6= 20. (2.2) (H0:V= 120; H1:V6= 120. (2.3) Hemos visto con anterioridad que usaríamos para estos test los estadísticos: TKni =ˆ Kni−K0 √AV ar(ˆ Kn), TVni =ˆ Vni−V0 √AV ar(ˆ Vn)∼Tn−2; asumiendo normalidad de los estimadores, y en donde hemos considerado esta vez σ , y K y V desconocidos, de tal forma que hemos calculado las varianzas de sus distribuciones 2.6. DISTRIBUCIONES 25 asintóticas como: AV ar(ˆ Kn) = ˆσ2 n nan(ˆ Kn) y AV ar(ˆ Vn) = ˆσ2 n nan(ˆ Vn) . Obviamente, los estimadores no siempre son normales, por lo que cabría preguntarse como evoluciona el nivel de estos tests en función del tamaño muestral. Dada la muestra simulada de S estimadores, vamos a calcular cada TKni y TVni para así poder averiguar con cuánta frecuencia los estadísticos no verican que: −tα≤ˆ T≤tα, es decir, cuantos estadísticos caen en la región de rechazo, siendo tα un valor de una T de Student de n−2 grados de libertad, que deja a su derecha una región de probabilidad α . Podemos ver en los cuadros 2.6 y 2.7 los diferentes porcentajes de rechazo de los contrastes para distintos valores de σ . n H0:K= 20 H0:V= 129 15 6.5 5.1 30 5.6 4.9 Cuadro 2.6: Porcentajes de rechazo de los dos contrastes según el tamaño muestral con nivel de signicación del 5%, σ= 15 y 100000 simulaciones n H0:K= 20 H0:V= 129 15 5.1 4.9 30 5.1 5.0 Cuadro 2.7: Porcentaje de rechazo de los dos contrastes según el tamaño muestral con nivel de signicación del 5%, σ= 5 y 100000 simulaciones Estos resultados coinciden con lo obtenido en la sección 2.4 : al reducir la varianza, las aproximaciones asintóticas mejoran. 2.6. Distribuciones Durante todo momento hemos considerado que la distribución de ˆ θn se aproxima por una normal de media θ y varianza AV ar(ˆ θn) . Hemos, incluso, construido intervalos y tests de contraste basados en esto, cuyos niveles mejoraban al aumentar el tamaño muestral n . Consideremos T= (K,V) la muestra de estimadores simulados. Para estudiar la densidad de éstos vamos a tomar la estandarización de las muestras K y V , de tal forma que de seguir K y V una distribución normal, la muestra estandarizada tendrá que seguir una normal de 26 CAPÍTULO 2. SIMULACIÓN media 0 y varianza 1 ; y vamos a utilizar un estimador de la función de densidad tipo núcleo. (a) Densidad aproximada de ˆ Kn (b) Densidad aproximada de ˆ Vn Gráca 2.5: Funciones de densidad de ˆ Kn y ˆ Vn , con σ= 15 y 100000 simulaciones. En la gráca 2.5 podemos ver cierta tendencia a una distribución normal al aumentar el tamaño muestral, algo coherente con los resultados obtenidos en las secciones 2.4 y 2.5. Ahora bien, considerando una varianza menor, los resultados obtenidos en apartados anteriores sugerirían una mayor normalidad de los estimadores, por lo que sería interesante estudiar estas densidades, pero para un σ menor. Los resultados vistos en la gráca 2.6 (a) Densidad aproximada de ˆ Kn (b) Distribución aproximada de ˆ Vn Gráca 2.6: Funciones de densidad de ˆ Kn y ˆ Vn , con σ= 5 y 100000 simulaciones. sugieren que una disminución de la varianza de los errores conlleva una convergencia más 2.6. DISTRIBUCIONES 27 (a) ˆ K, n = 15 (b) ˆ K, n = 30 (c) ˆ V , n = 15 (d) ˆ V , n = 30 Gráca 2.7: Q-Q Plot de ambos estimadores para n= 15 y n= 30 , con σ= 15 . rápida a una normal. Podemos ver estos resultados con más claridad mirando los grácos cuantil-cuantil de los diferentes estimadores en la gráca 2.7 (a mas alejados de la diagonal, mayor discrepancia respecto de la normal). Cabe destacar que para σ= 5 , tal y como se ve en la gráca 2.8, los cuantiles más grandes y los más pequeños se acercan mucho más a su valor teórico, indicándonos que los estimadores son cada vez más simétricos. Observación 2.3 . En las grácas Cuantil-Cuantil de una muestra de tamaño n (también llamadas Q-Q Plots) representamos n+ 1 puntos, dados por los (n+ 1) -cuantiles de la muestra con los correspondientes cuantiles de la distribución de comparación (en este caso la normal de media 0 y varianza 1). De esta manera, si la distribución de la muestra fuese normal, sus cuantiles tendrían que coincidir con los teóricos, acercándose a la diagonal. 34 CÓDIGO DE R EMPLEADO Sigma=cbind(c(sum(Mk^2),sum(Mv*Mk)),c(sum(Mv*Mk),sum(Mv^2))) Sigma=solve(Sigma) Sigma=Sigma*sigma^2 VAK[i,]=Sigma[1,1] VAV[i,]=Sigma[2,2] VAC[i,]=Sigma[1,2] } #Calculo de la varianza aproximada por simulacion v=matrix(0,ncol=2,nrow=m) int=matrix(0,ncol=2,nrow=m) for (i in 1:m){ set.seed(123456) n=5*i x=seq(10,110,length=n) theta=matrix(0,nrow=ns,ncol=2) for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) } int[i,]=sigma.test(theta[,1])$conf.int[1:2] v[i,]=apply(theta, 2,var) } #Representacion de las graficas v/cbind(VAK,VAV) par(mfrow=c(1,2)) plot(v[,1]/VAK~Index[1:10],xlab='Tamano de la muestra',ylab='Razon de varianzas de K (Var. muestral sobre asintotitca)',main='100.000 simulaciones ') lines(v[,1]/VAK~Index[1:10]) plot(v[,2]/VAV~Index[1:10],xlab='Tamano de la muestra',ylab='Razon de varianzas de V (Var. muestral sobre asintotitca)',main='100.000 simulaciones ') lines(v[,2]/VAV~Index[1:10]) #Intervalos de confianza para la varianza de theta[,1] int Intervalos de Conanza 35 #####COBERTURA DE LOS INTERVALOS pRRS Y ASINTOTICO en 2.4#### set.seed(123456) ns=20000 V=120 K=20 m=2 theta=matrix(0,nrow=ns,ncol=2) pRSS_IC=matrix(0,nrow=m,ncol=2) ICA=matrix(0,nrow=m,ncol=2) contador=matrix(0,nrow=ns,ncol=2) contador2=matrix(0,nrow=ns,ncol=2) difeICA=matrix(0,ncol=2,nrow=m) difeRSS=matrix(0,ncol=2,nrow=m) #para sigma=5 sigma=5 for (i in 1:m){ n=15*i# Tamano muestral x=seq(10,110,length=n) pRSS_ICs=matrix(0,nrow=ns,ncol=4)#Intervalos mediante el pRSS ICAs=matrix(0,nrow=ns,ncol=4)#Intervalos asintoticos for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) pRSS_ICs[is,]=c(confint(m)[1,],confint(m)[2,]) #pRSS ICAs[is ,]=c(theta[is,1]+qt(0.05/2,n-2)*summary(m)$coefficients[1,2], theta[is,1]-qt(0.05/2,n-2)*summary(m)$coefficients[1,2], theta[is,2]+qt(0.05/2,n-2)*summary(m)$coefficients[2,2],theta[is,2]- qt(0.05/2,n-2)*summary(m)$coefficients[2,2]) if (20 >= pRSS_ICs[is,1] & 20 <= pRSS_ICs[is,2]){ contador[is,1]=1 36 CÓDIGO DE R EMPLEADO } if (120 >= pRSS_ICs[is,3] & 120 <= pRSS_ICs[is,4]){ contador[is,2]=1 } if (20 >= ICAs[is,1] & 20 <= ICAs[is,2]){ contador2[is,1]=1 } if (120 >= ICAs[is,3] & 120 <= ICAs[is,4]){ contador2[is,2]=1 } } pRSS_IC[i,]=colSums(contador) ICA[i,]=colSums(contador2) difeRSS[i,1]=(pRSS_ICs[,2]-pRSS_ICs[,1])/ns difeICA[i,1]=(ICAs[,2]-ICAs[,1])/ns difeRSS[i,2]=(pRSS_ICs[,4]-pRSS_ICs[,3])/ns difeICA[i,2]=(ICAs[,4]-ICAs[,3])/ns } pRSS_IC ICA difeRSS difeICA #para sigma=15 sigma=15 for (i in 1:m){ n=15*i# Tamano muestral x=seq(10,110,length=n) pRSS_ICs=matrix(0,nrow=ns,ncol=4)#Intervalos mediante el pRSS ICAs=matrix(0,nrow=ns,ncol=4)#Intervalos asintoticos 37 for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) pRSS_ICs[is,]=c(confint(m)[1,],confint(m)[2,]) #pRSS ICAs[is ,]=c(theta[is,1]+qt(0.05/2,n-2)*summary(m)$coefficients[1,2], theta[is,1]-qt(0.05/2,n-2)*summary(m)$coefficients[1,2],theta[is,2]+ qt(0.05/2,n-2)*summary(m)$coefficients[2,2],theta[is,2]-qt(0.05/2,n-2)* summary(m)$coefficients[2,2]) if (20 >= pRSS_ICs[is,1] & 20 <= pRSS_ICs[is,2]){ contador[is,1]=1 } if (120 >= pRSS_ICs[is,3] & 120 <= pRSS_ICs[is,4]){ contador[is,2]=1 } if (20 >= ICAs[is,1] & 20 <= ICAs[is,2]){ contador2[is,1]=1 } if (120 >= ICAs[is,3] & 120 <= ICAs[is,4]){ contador2[is,2]=1 } } pRSS_IC[i,]=colSums(contador) ICA[i,]=colSums(contador2) difeRSS[i,1]=(pRSS_ICs[,2]-pRSS_ICs[,1])/ns difeICA[i,1]=(ICAs[,2]-ICAs[,1])/ns difeRSS[i,2]=(pRSS_ICs[,4]-pRSS_ICs[,3])/ns difeICA[i,2]=(ICAs[,4]-ICAs[,3])/ns } pRSS_IC ICA 38 CÓDIGO DE R EMPLEADO difeRSS difeICA Representación de los perles y las razones de asimetría ######pRSS en 2.4########## Asim=matrix(0,ncol=2,nrow=4) vn=c(15,30) vs=c(15,5) j=0 for (sigma in vs){ for (n in vn){ set.seed(123456) j=j+1 p=50 q=250 x=seq(10,110,length=n) K=20 V=120 y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta=matrix(0,nrow=1,ncol=2) theta=coef(m) RSSK=matrix(0,ncol=1,nrow=q-59) RSSV=matrix(0,ncol=1,nrow=p+1) for (i in 60:q){ V=i RSSK[i-59]=deviance(nls(y~V*x/(K+x),start=list(K=20))) } V=120 for (i in 0:p){ K=i RSSV[i+1]=deviance(nls(y~V*x/(K+x),start=list(V=120))) } 39 K=20 #Representemos la funcion IndexK=60:q IndexV=0:p plot(RSSK~IndexK ,xlab='Valor de V',ylab=' Residuo ',main='perfil RSS para K=20',xlim=c(75,175),ylim=c(0,15000)) #Valor del estimador de V V=coef(m)[2] points(c(coef(m)[2]),c(deviance(nls(y~V*x/(K+x),start=list(K=20)))), col="red",pch=9) #Intervalo de V sobre la funcion V=confint(m)[2,] y1=deviance(nls(y~V[1]*x/(K+x),start=list(K=20))) y2=deviance(nls(y~V[2]*x/(K+x),start=list(K=20))) points(confint(m)[2,],c(y1,y2),col="blue",pch=19) lines(confint(m)[2,],c(y1,y2),col="blue",type='o') #Recta points(c(coef(m)[2]),c(y1),col="blue",pch=3) plot(RSSV~IndexV ,xlab='Valor de K',ylab='Residuo',main='perfil RSS para V=120',xlim=c(0,50),ylim=c(0,15000)) #Valor del estimador de K K=coef(m)[1] points(c(coef(m)[1]),c(deviance(nls(y~V*x/(K+x),start=list(V=120)))), col="red",pch=9) #Intervalo de K sobre la funcion K=confint(m)[1,] y1=deviance(nls(y~V*x/(K[1]+x),start=list(V=120))) y2=deviance(nls(y~V*x/(K[2]+x),start=list(V=120))) points(confint(m)[1,],c(y1,y2),col="blue",pch=19) 40 CÓDIGO DE R EMPLEADO lines(confint(m)[1,],c(y1,y2),col="blue",type='o') #Asimetrias Asim[i,]=c((coef(m)[1]-K[1])/(K[2]-coef(m)[1]),(coef(m)[2]-confint(m)[2,1])/(confint(m)[2,2]-coef(m)[2])) #Recta points(c(coef(m)[1]),c(y1),col="blue",pch=3) } } Contrastes ######CONTRASTES DE 2.5###### set.seed(123456) par(mfrow=c(1,1)) pval=matrix(0,ncol=2,nrow=ns) ns=100000 theta=matrix(0,nrow=ns,ncol=2) K=20 V=120 nivelsd5=matrix(0,ncol=2,nrow=2) nivelsd15=matrix(0,ncol=2,nrow=2) contador=matrix(0,nrow=ns,ncol=2) ThK=matrix(0,ncol=1,nrow=ns) ThV=matrix(0,ncol=1,nrow=ns) sigma=5 i=0 for (n in vn){ i=i+1 set.seed(123456) x=seq(10,110,length=n) for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) ThK=(summary(m)$coef[,1][1] - 20)/(summary(m)$coef[,2][1]) ThV=(summary(m)$coef[,1][2] - 120)/(summary(m)$coef[,2][2]) 41 pval[is ,]=c(2*pt(-abs(ThK),n-2),2*pt(-abs(ThV),n-2)) } nivelsd5[i,1]=length(pval[1:is,1][pval[1:is,1]<=0.05])/(is) nivelsd5[i,2]=length(pval[1:is,2][pval[1:is,2]<=0.05])/(is) } i=0 sigma=15 for (n in vn){ i=i+1 set.seed(123456) x=seq(10,110,length=n) for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) ThK=(summary(m)$coef[,1][1] - 20)/(summary(m)$coef[,2][1]) ThV=(summary(m)$coef[,1][2] - 120)/(summary(m)$coef[,2][2]) pval[is ,]=c(2*pt(-abs(ThK),n-2),2*pt(-abs(ThV),n-2)) } nivelsd15[i,1]=length(pval[1:is,1][pval[1:is,1]<=0.05])/(is) nivelsd15[i,2]=length(pval[1:is,2][pval[1:is,2]<=0.05])/(is) } Distribuciones para errores normales y sus Q-Q Plot ######DISTRIBUCIONES DE 2.6######### set.seed(123456) ns=100000 vs=c(5,15) for (sigma in vs){ n=15 # Tamano muestral x=seq(10,110,length=n) theta=matrix(0,nrow=ns,ncol=2) theta15=matrix(0,nrow=ns,ncol=2) 42 CÓDIGO DE R EMPLEADO for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) } theta15=theta aux=theta15[,1] auxe=(aux-mean(aux))/sd(aux) plot(density(auxe),col='dark green',type='l',lwd=1,xlim=c(-4,4), xlab="Valor de x", ylab="Densidad estandarizada", main="Comparacion de las distribuciones al aumentar n") x <- seq(-4,4, length=100) hx <- dnorm(x) degf <- c(1,3,8,30) colors <- c("red","blue","darkgreen","gold","black") labels <- c("df=1","df=3","df=8","df=30","normal") lines(x, hx, type="l", lty=2, xlab="Valor de x", ylab="Densidad estandarizada", main="Comparaci on de las distribuciones al aumentar n") set.seed(123456) n=30 # Tama no muestral x=seq(10,110,length=n) theta=matrix(0,nrow=ns,ncol=2) theta30=matrix(0,nrow=ns,ncol=2) for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) } aux=theta[,1] auxe=(aux-mean(aux))/sd(aux) lines(density(auxe)) theta30[,1]=theta[,1] theta30[,2]=theta[,2] 43 set.seed(123456) n=50 # Tama no muestral x=seq(10,110,length=n) theta=matrix(0,nrow=ns,ncol=2) for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) } aux=theta[,1] auxe=(aux-mean(aux))/sd(aux) lines(density(auxe),col=4,type='l') set.seed(123456) n=100 # Tama no muestral x=seq(10,110,length=n) theta=matrix(0,nrow=ns,ncol=2) for (is in 1:ns){ y=V*x/(K+x)+rnorm(n,sd=sigma) m=nls(y~V*x/(K+x),start=list(K=20,V=120)) theta[is,]=coef(m) } aux=theta[,1] auxe=(aux-mean(aux))/sd(aux) lines(density(auxe),col=6) colors=c('dark green','black',4,6) labels=c('n=15','n=30','n=50','n=100') legend("topright", inset=.05,title="Distribuciones", labels, lwd=2, lty=c(1,1,1,1,2), col=colors) #######QQ PLOTS######### aux=theta15[,2] auxe=(aux-mean(aux))/sd(aux) qqnorm(auxe ,main='Q-Q Plot del estimador de V con n=15',ylim=c(-4,10)) abline(a=0,b=1)