Transcript
UNIVERSIDAD AUTÓNOMA CHAPINGO DIVISIÓN DE CIENCIAS FORESTALES “DISCRIMINACIÓN ENTRE LAS DISTRIBUCIONES GAUSIANA INVERSA Y WEIBULL Y UNA APLICACIÓN EN BIOMETRÍA FORESTAL.” TESIS PROFESIONAL Que como requisito parcial para obtener el título de: LICENCIADO EN ESTADÍSTICA PRESENTA: Acosta Percastegui Alan
Chapingo, México. Julio del 2009.
La presente tesis titulada "Discriminación entre las distribuciones Gausiana Inversa y Weibull y una aplicación en Biometría Forestal" fue realizada por el C. Alan Acosta Percastegui, bajo la dirección del Dr. Carlos L. Cíntora González. Ha sido revisada y aprobada por el Comité Revisor y Jurado examinador que se indica para obtener el título de Licenciado en Estadística.
Presidente
Dr. Carlos L. Cíntora González
Secretario
Dr. Gerardo H. Terrazas González
Vocal
Lc. Margarito Soriano Montero
Suplente
M. en C. Alejandro Corona Ambriz
Suplente
M. en C. Ángel Leyva Ovalle
AGRADECIMIENTOS
A Dios, por permitirme concluir este primer paso en mi formación profesional, y por poder compartir con mi familia y amigos este logro.
A mis padres, Joel Acosta Martínez y Yolanda Percastegui González, por su apoyo incondicional, sin importar que tan difíciles fueran los tiempos, por su amor de padres para con sus hijos, MUCHAS GRACIAS.
A mis hermanas, Arianna y Jeannette, por todo su apoyo y cariño.
A toda mi familia, que por cuestiones de espacio no los nombro, pero sé que cuento con su apoyo.
A la Universidad Autónoma Chapingo, mi alma mater, por la invaluable oportunidad que me brindó para formarme profesionalmente.
Al Dr. Carlos L. Cíntora González, por su asesoría constante, tiempo y esfuerzo prestados para la realización de la presente.
Al Dr. Gerardo H. Terrazas González, Lic. Margarito Soriano Montero, M. en C. Ángel Leyva Ovalle y al M. en C. Alejandro Corona Ambriz, por la revisión, y puntuales observaciones para la realización de la presente.
A mis amigos: Arisel, Kaliman, Fer, Lore, Chumi, Eli, Francisco, Roy, Ramón, José, Mari, Micah, a los miembros de la generación 2004-2008 de la Lic. en Estadística, una disculpa a los que omite por cuestiones de espacio, pero sé que cuento con su inapreciable amistad.
“Cada uno de nosotros lo sabe todo. Sólo necesitamos abrir nuestras mentes para escuchar nuestra propia sabiduría”
´ Indice general
´ Indice de Cuadros
IV
´ Indice de Figuras
VI
Resumen
VII
Abstract
IX
1. Introduccion ´
1
2. Objetivos
3
3. Revision ´ de Literatura
4
3.1. Distribucion ´ Gausiana Inversa . . . . . . . . . . . . . . . . . . . . . . . . .
4
3.1.1. Propiedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
3.1.2. Estimacion ´ por M´axima Verosimilitud . . . . . . . . . . . . . . . .
7
3.2. Distribucion ´ Weibull . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3.2.1. Propiedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
3.2.2. Estimacion ´ por M´axima Verosimilitud . . . . . . . . . . . . . . . .
11
3.3. Informacion ´ de Kullback-Leibler o Distancia Entre dos Modelos . . . .
13
I
3.3.1. La realidad f , considerada como constante . . . . . . . . . . . . . .
14
3.4. Criterio de Informacion ´ de Akaike . . . . . . . . . . . . . . . . . . . . . .
16
3.4.1. Criterio de Informacion ´ de segundo orden
. . . . . . . . . . . . .
19
3.5. Prueba de Vuong . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
3.6. M´etodo de Simulacion ´ de Monte Carlo . . . . . . . . . . . . . . . . . . . .
22
4. Metodolog´ıa
24
5. Aplicacion, ´ Simulacion ´ y Discusion ´
27
5.1. Aplicacion ´ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
5.1.1. Estimacion ´ Distribucion ´ Weibull . . . . . . . . . . . . . . . . . . .
28
5.1.2. Estimacion ´ Distribucion ´ Gausiana Inversa . . . . . . . . . . . . .
31
5.1.3.
Criterio de Informacion ´ AICc . . . . . . . . . . . . . . . . . . . . .
33
5.1.4. Estad´ıstico de Vuong . . . . . . . . . . . . . . . . . . . . . . . . . .
33
5.1.5. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
5.2. Simulacion ´ y Evaluacion ´ del Comportamiento de los Estad´ısticos . . . .
36
5.2.1. Simulacion ´ del criterio AICc . . . . . . . . . . . . . . . . . . . . . .
36
5.2.2. Simulacion ´ del criterio de Vuong . . . . . . . . . . . . . . . . . . .
41
5.3. Discusion ´ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
45
6. Conclusiones y Recomendaciones
48
´ APENDICES
51
A. Codigo ´ R usado para calcular el criterio AICc
52
II
B. Codigo ´ R usado para la simulacion ´ de Monte Carlo
56
Referencias
66
III
´ Indice de cuadros ´ forestal ex4.1. Di´ametros normales de n = 67 a´ rboles de pino en la estacion perimental Zoquiapan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
5.1. Valores estimados de AICc para cada una de las densidades . . . . . . . .
33
5.2. Fracciones de rechazo observadas. Modelo verdadero IG (µ = 0.5, λ = ´ basada en el criterio AICc con nivel de significancia 1), y discriminacion nominal de α = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
5.3. Fracciones de rechazo observadas. Modelo verdadero IG (µ = 0.5, λ = ´ basada en el criterio AICc con nivel de significancia 0.8), y discriminacion α = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
5.4. Fracciones de rechazo observadas. Modelo verdadero IG (µ = 0.5, λ = ´ basada en el criterio AICc con nivel de significancia 2.5), y discriminacion α = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
5.5. Fracciones de rechazo observadas. Modelo verdadero Weibull (c = 1.5, α = ´ basada en el criterio AICc con nivel de significancia 0.5), y discriminacion de α = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
5.6. Fracciones de rechazo observadas. Modelo verdadero Weibull (c = 1.5, α = ´ basada en el criterio AICc con nivel de significancia 0.95), y discriminacion α = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
IV
39
5.7. Fracciones de rechazo observadas. Modelo verdadero Weibull (c = 1.5, α = ´ basada en el criterio AICc con nivel de significancia 2.5) y discriminacion α = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
´ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8. Fracciones de Seleccion
40
5.9. Fracciones de rechazo, modelo verdadero IG (µ = 0.8, λ = 1.5), criterio de ´ de Vuong con nivel de significancia α = 0.05 . . . . . . . . discriminacion
42
5.10. Fracciones de rechazo, modelo verdadero Weibull (c = 4, α = 1.5), criterio ´ de Vuong con nivel de significancia α = 0.05 . . . . . . de discriminacion
43
´ 5.11. Fracciones de rechazo, modelo verdadero es χ215 , criterio de discriminacion de Vuong con nivel de significancia α = 0.05 . . . . . . . . . . . . . . . . .
43
5.12. Fracciones de rechazo, modelo verdadero IG (µ = 0.8, λ = 1.5), criterio de ´ de Vuong con nivel de significancia α = 0.05 . . . . . . . . discriminacion
44
5.13. Fracciones de rechazo, modelo verdadero Weibull (c = 4, α = 1.5), criterio ´ Vuong . . . . . . . . . . . . . . . . . . . . . . . . . . . . . de discriminacion
44
´ 5.14. Fracciones de rechazo, modelo verdadero χ215 , criterio de discriminacion Vuong . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
V
45
´ Indice de figuras ´ IG (µ = 1, λ = 41 ) y IG (µ = 1,λ = 1) . . . . . . . . . . . . . . . . . 3.1. Distribucion
5
´ Weibull (c = 1.5, α = 1.25) y Weibull (c = 2,α = 1.25) respectiva3.2. Distribucion mente (con ξ 0 = 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
5.1. Histograma de n = 67 observaciones de di´ametros normales . . . . . . . . . . .
27
´ estimada Weibull (cˆ = 1.8291474, αˆ = 0.292557) . . . . . . . . . . . 5.2. Distribucion
30
´ estimada Weibull (cˆ = 1.8291474, αˆ = 0.292557) vs histograma muestral 31 5.3. Distribucion ´ estimada IG (µˆ = 0.2576866, λˆ = 0.9933674) . . . . . . . . . . . . . 5.4. Distribucion
32
´ estimada IG (µˆ = 0.2576866, λˆ = 0.9933674) vs histograma muestral 35 5.5. Distribucion
VI
Resumen El presente trabajo tuvo como finalidad estudiar el comportamiento de dos criterios para ´ discriminar entre las densidades Weibull y Gausiana Inversa. Los criterios de seleccion de modelos son los estad´ısticos AICc y Vuong, y se aplicaron a una muestra de 67 ob´ de ambas servaciones correspondientes a di´ametros normales. Se realizo´ la estimacion distribuciones por el m´etodo de m´axima verosimilitud, y el c´alculo de ambos criterios de ´ fue programado en el software R, resultando mejor modelo la distribudiscriminacion ´ Gausiana Inversa con par´ametros estimados µˆ = 0.2576866 y λˆ = 0.9933674 con cion un valor de AICc = −104.0795. La prueba de Vuong rechazo´ la igualdad de modelos en ´ Gausiana Inversa con un nivel de significancia de α = 0.05. Se anafavor de distribucion lizo´ el comportamiento de ambos criterios por medio de simulaciones de Monte Carlo, para lo cual se fijaron pseudo poblaciones determinadas por modelos Weibull y Gausiano ˜ 10,20,30,50 y 60 respecInverso, para los cuales se simularon 10,000 muestras de tamano tivamente. Para cada una de las pseudo poblaciones se estimaron ambos modelos y se ´ en la que la pseudo poblacion ´ evaluaron los criterios AICc y Vuong. Para la situacion ´ de rechazo del criterio estaba determinada por un modelo Gausiano Inverso la fraccion ˜ (cercana a 0), es decir, el nivel de significancia observado fue AICc y Vuong fue pequena menor o igual que el nivel nominal de ambas pruebas, es decir, ambos criterios eleg´ıan ´ Weibull, cuando al modelo Gausiano Inverso como mejor modelo sobre la distribucion ´ en la que la pseuen realidad e´ ste es el que genero´ la muestra aleatoria. Para la situacion ´ estaba determinada por alguna distribucion ´ Weibull el comportamiento fue do poblacion ´ para todas las situaciones en las que la pseudo poblacion ´ estaba detersimilar. M´as aun, ´ modelo Gausiano Inverso, se encontro´ que con un tamano ˜ de muestra minada por algun de n = 30 y mayor, el nivel de significancia observado es igual o menor que el nivel nominal en ambas pruebas (AICc y Vuong). Cada una de las rutinas fueron escritas en el software R, y se presentan en los ap´endices A y B. VII
Palabras Clave: AICc , Vuong, Monte Carlo, densidad Gausiana Inversa, densidad Weibull, R.
VIII
Abstract The behavior of two criteria for discriminating between an Inverse Gaussian and a Weibull density was studied. The criteria are based on the AICc statistic and Vuong’s test, and were applied to a random sample of 67 diameters at breast height obtained from pines ´ Forestal Experimental Zoquiapan, M´exico. Maximum likelihood esgrowing at Estacion timates of the parameters were obtained and the computations were implemented using the R language; numerical optimization was used to estimate the Weibull density, whereas analytic expressions for the estimators were available for the Inverse Gaussian density. The best model, according to the discrimination criteria, was the Inverse Gaussian denb = 0.2576866, b sity with estimates µ λ = 0.9933674 and AICc = −104.0795. Vuong’s test rejected the null hypothesis of equal models assumption, selecting the Inverse Gaussian over the Weibull density with a nominal significance level of α = 0.05. To study the behavior of the above discrimination criteria, as a function of the sample size, Monte Carlo simulations were performed; pseudo-populations, determined by the Inverse Gaussian or Weibull density, were considered and 10,000 random samples of sizes 10, 20, 30, 50 and 60 respectively were simulated from each of the above densities. Using each simulated sample the maximum likelihood estimates and the AICc and Vuong’s statistic were estimated. When the pseudo-population corresponded to an Inverse Gaussian model, the fraction of times that the null hypothesis was rejected, when it was indeed true, was less or equal than the nominal significance level, for both tests, i.e. the Inverse Gaussian was selected as the best model when the sample was generated from this model. The same behavior was observed when the true model was the Weibull density. It was generally found that, if the sample was generated from the Inverse Gaussian model, then a sample size of 30 was enough to attain an empirical significance level close to the nominal one for both, the AICc and Vuong test. Key words: AICc , Vuong, Monte Carlo, Inverse Gausian, Weibull, R. IX
Cap´ıtulo 1 Introduccion ´ ´ Weibull, as´ı como la Gausiana Inversa, o distribucion ´ de Wald, son La distribucion modelos de probabilidad ampliamente utilizados en situaciones en las cuales el soporte ´ de la variable aleatoria de inter´es son los numeros reales positivos. Dentro de esta gama de situaciones se encuentran las de los datos relacionados con tiempos de vida y an´alisis de supervivencia. ´ Weibull, obtenida por el ingeniero Sueco Walodi Weibull, ha sido La distribucion ampliamente usada en an´alisis de tiempos de vida, ver por ejemplo Lawless (1982), Kalbfleisch and Prentice (2002),o Klein and Moeschberger (2003), y ha sido estudiada extensivamente por diversos autores, entre los cuales pueden citarse a Smith and Weissman (1985), Bain and Antle (1969), y otros m´as. En su forma simple, el modelo Weibull consta de dos par´ametros, uno de escala y otro ´ pero en este caso ocurre a de forma, y puede considerarse uno adicional, de localizacion, ´ de verosimilitud falla en proporcionar los estimadores de m´axima menudo que la funcion ˜ y el par´ametro de forma es menor verosimilitud cuando las observaciones son pequenas que la unidad, Smith and Weissman (1985). La densidad Gausiana Inversa, por su parte, es una densidad asociada a una variable aleatoria continua con soporte en los reales positivos, la cual pertenece a la familia expo1
nencial multi-param´etrica, ver Fahrmeir and Tutz (2001), o Dobson (2001). La densidad ha sido usada en diversas aplicaciones, Jorgensen (1982), Seshadri (1993), pero parece no haber sido usada en aplicaciones forestales. El problema de discriminar entre modelos estad´ısticos competitivos ocurre cuando se tiene un conjunto de observaciones, correspondientes a realizaciones de variables aleatorias observables, y se dispone de dos o m´as modelos candidatos que pueden usarse para ´ consiste modelar las respuestas observadas. Dadas las observaciones, la discriminacion en escoger uno de los modelos competitivos, usando como criterio una medida de la ´ entre las medidas respectivas obtenidas. calidad del ajuste del modelo y la comparacion Aunque el problema de discriminar entre modelos competitivos ha sido estudiado en ´ relativa al problema de algunos casos, Cox (1961), Vuong (1989), no se tiene informacion discriminar entre los modelos Gausiano inverso y Weibull y, por ende, se considera de ´ al respecto. inter´es efectuar un estudio que proporcione informacion
2
Cap´ıtulo 2 Objetivos En el presente trabajo se tienen los siguientes objetivos: ´ que permitan escoger entre los 1. Establecer criterios estad´ısticos de discriminacion modelos Weibull y Gausiano inverso. ´ Monte Carlo el desempeno ˜ de los criterios en dis2. Estudiar mediante simulacion ´ cusion. ´ de los criterios para la discriminacion ´ entre dos modelos de 3. Mostrar una aplicacion ´ de distribuciones diam´etricas de a´ rboles. utilidad para la descripcion
3
Cap´ıtulo 3 Revision ´ de Literatura 3.1.
Distribucion ´ Gausiana Inversa
La variable aleatoria X sigue la densidad Gausiana Inversa con par´ametros µ y λ si ¾ λ 2 exp − 2 ( x − µ) f X ( x | µ, λ) = 2µ x ¸1/2 ½ µ ¶¾ · λ x µ λ exp − −2+ = 2µ µ x 2πx3 ·
λ 2πx3
¸1/2
½
(3.1) (3.2)
donde µ > 0, λ > 0 y x > 0. Cuando X sigue la densidad anterior se anota X ∼ IG (µ, λ). En la Figura (3.1) se muestran dos casos espec´ıficos de tal densidad.
4
2.5
IG(µ = 1 λ = 1 4)
0.0
0.5
1.0
fX(x)
1.5
2.0
IG(µ = 1 λ = 1)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
x
´ IG (µ = 1, λ = 41 ) y IG (µ = 1,λ = 1) Figura 3.1: Distribucion
La densidad (3.2) puede escribirse en cualquiera de las siguientes tres formas: ·
µφ f x ( x | µ, φ) = 2πx3
¸1/2
½
1 e exp − φ 2 φ
µ
x µ + µ x
¶¾ ,
(3.3)
¸ ½ µ ¶¾ 1 φ2 x λ λ 1/2 φ e exp − + f x ( x | φ, µ) = , 2 λ x 2πx3 ½ ¾¸ ¸ · · 1 λ λ 1/2 1/2 λx − (2λ) + exp − f x ( x | α, λ) = , 2 2x 2πx3 ·
(3.4) (3.5) 2
a las cuales se les representar´a respectivamente por IG (µ, φµ), IG ( λφ , λ) e IG ( α2 , λ) rep spectivamente. En las ecuaciones anteriores φ = λ/µ, α = (2µ), mientras que λ,µ,φ, y α son todas positivas. ´ de distribucion ´ acumulada correspondiente a (3.2) es La funcion ( r µ (r µ ¶) ¶) x λ x λ −1 + e2λ/µ Φ − +1 FX ( x | µ, λ) = Φ x µ x µ 5
(3.6)
´ de distribucion ´ normal est´andar. donde Φ(·) representa la funcion ´ generadora de momentos, es decir De (3.2) se obtiene la funcion ( µ ¶1/2 ) 2µ2 t λ 1− 1− , Ψ X (t; µ, λ) = log m X (t) = µ λ
(3.7)
´ caracter´ıstica donde m X (t) = E(etX ), t = 1, 2, .., n, es decir, t-´esimo momento; la funcion correspondiente es:
"
λ exp µ
(
µ
2iµ2 t 1− 1− λ
¶1/2 )# .
(3.8)
Las funciones generadoras acumuladas para (3.3),(3.4) y (3.5) respectivamente son: ( ¶ ) µ 2µt 1/2 (3.9) Ψ X (t; µ, φ) = φ 1 − 1 − φ ( ¶ ) µ 2λt 1/2 (3.10) Ψ X (t; φ, λ) = φ 1 − 1 − 2 φ ( ¶1/2 ) µ t , Ψ X (t; α, λ) = 21/2 λ α1/2 − α − (3.11) λ y los primeros cuatro momentos correspondientes a (3.2) son k1 = µ µ3 k2 = λ k3 = k4 =
3µ5 λ2
5µ7 λ3
En general, para r ≥ 2 µ2r−1
kr = 1 × 3 × 5 × ... × (2r − 3) λr−1
donde kr (.) representa el r-´esimo momento de la densidad Gausiana Inversa, ver Johnson et al. (1994). 6
3.1.1.
Propiedades
Algunas propiedades de la densidad Gausiana Inversa son las siguientes: 1) Si X ∼ IG (µ, λ) y a > 0, entonces aX ∼ IG ( aµ, aλ). ´ IG (µi , λi ) 2) Considere las variables X1 , ..., Xn mutuamente independientes, con distribucion para i = 1, ..., n, entonces n
∑ (µi−2 λi Xi ) ∼ IG(µ, λ)
i =1
donde
µ
n
µ=
∑
i =1
λi µi
¶ y
λ = µ2 .
3) Si µi = µ y λi = λ para i = 1, ..., n, se cumple que λ µ2
n
λ
∑ Xi = µ2 nX¯ ∼ IG(nλ/µ, n2 λ2 /µ2 )
i =1
donde X¯ ∼ IG (µ, nλ). 4) Si las variables X1 , X2 , ..., Xn son mutuamente independientes e id´enticamente distribuidas IG (µ, λ), entonces Q1 = y
nλ( X¯ − µ)2 µ2 X¯ (
Q − Q1 = λ
n
)
∑ Xi−1 − nX¯ −1
i =1
son mutuamente independientes y Q − Q1 se distribuye como una variable aleatoria con densidad χ2n−1 . Ver Johnson et al. (1994).
3.1.2.
Estimacion ´ por M´axima Verosimilitud
Los estimadores de m´axima verosimilitud (EMV) se obtienen en la forma usual; sean X1 , ..., Xn una serie de observaciones provenientes de distribuciones Gausianas Inversas 7
IG (µ, λi ) para i = 1, ..., n y sea λi = wi λ0 , donde λ0 es desconocida, pero wi es un valor positivo conocido. Los EMV de µ y λ satisfacen las ecuaciones " #" # −1 n
µˆ =
∑ w i Xi
i =1
n
∑ wi
(3.12)
i =1
1 n 1 = ∑ wi ( Xi−1 − X¯ −1 ) n i =1 λˆ
(3.13)
Cuando las wi son todas iguales a la unidad, entonces (3.12) y (3.13) se convierten respectivamente en ∑n X µˆ = X¯ = i=1 i n
(3.14)
n 1 = V = n−1 ∑ ( Xi−1 − X¯ −1 ) λˆ0
(3.15)
i =1
´ es IG (µ, nλ0 ), Puede verificarse que X¯ es un estad´ıstico suficiente para µ, cuya distribucion y Tweedie (1957) demostro´ que X¯ y λˆ 0 son independientes. Puede mostrarse tambi´en que si V = λˆ 0−1 entonces
(λ0 n)V ∼ χ2n−1 . ´ (3.15) se obtiene la siguiente relacion ´ Combinando el resultado anterior y la expresion n
∑ wi (Xi−1 − X¯ −1 ) ∼ λ0−1 × (χ2 )
(3.16)
i =1
donde χ2 denota la densidad Ji-cuadrada con n − 1 grados de libertad. A partir de la ´ (3.16) es posible construir intervalos de confianza para λ0 , [Tweedie (1957)]. expresion Tambi´en note que n
(n − 1)−1 ∑ wi ( Xi−1 − X¯ −1 ) i =1
es un estimador insesgado de 1/λ0 , y es de hecho el estimador insesgado de m´ınima varianza, Roy and Wasan (1968). ´ para el estimador insesgado de 1/λ0 es S2 / X¯ 3 , donde Una aproximacion S2 =
∑in=1 ( Xi − X¯ −1 ) . n−1 8
3.2.
Distribucion ´ Weibull
´ Weibull, si para c > 0, α > 0, y ξ 0 se Se dir´a que la variable X sigue una distribucion cumple que:
µ Y=
X − ξ0 α
¶
∼ exp(1),
X>0
(3.17)
´ exponencial est´andar con funcion ´ de densidad donde exp(1) representa la distribucion de probabilidad f Y (y) = e−y , y > 0.
(3.18)
´ de densidad de probabilidad de la variable aleatoria X es De acuerdo con esto, la funcion c f X (x) = α
µ
x − ξ0 α
¶ c −1
c
e{( x−ξ 0 )/α} ,
x > ξ0.
(3.19)
La densidad (3.19) se denomina densidad Weibull con par´ametros c > 0, α > 0, y ξ 0 y se anotar´a X ∼ Weibull (c, α, ξ 0 ), Johnson et al. (1994). Si X ∼ Weibull (c, α, ξ 0 ), entonces la ´ de distribucion ´ acumulada es: funcion c
FX ( x ) = 1 − e−{( x−ξ 0 )/α} , x > ξ 0 .
(3.20)
´ de densidad Weibull en (3.19) tiende a 0 conforme x → ξ 0 , y hay Para c > 1 la funcion un solo modelo
µ x=α
c−1 c
¶1/c
+ ξ0.
(3.21)
Este valor tiende a α + ξ 0 r´apidamente cuando c → ∞. Para 0 < c ≤ 1 el modelo es ξ 0 y ´ decreciente en x para toda x > ξ 0 . la densidad es una funcion ´ es De (3.20) puede ver que la mediana de la distribucion α(log 2)1/c Adem´as note que para cualquier valor de c . FX (ξ 0 + α) = 1 − e−1 = 0.63 9
(3.22)
´ considera ξ 0 = 0 y α = 1, y as´ı la funcion ´ de densiLa forma est´andar de la distribucion dad est´andar es c
f X ( x ) = cx c−1 e− x , x > 0, c > 0,
(3.23)
´ de distribucion ´ acumulada es mientras que la correspondiente funcion c
FX ( x ) = 1 − e− x , c > 0, x > 0.
3.2.1.
(3.24)
Propiedades
´ Weibull de tres par´ametros (3.19) Los momentos correspondientes a la distribucion 0
´ X = ξ 0 + αX. pueden obtenerse a partir de la densidad (3.23)usando la transformacion Haciendo esto se obtiene
¶ 1 +1 E[ X ] = Γ c ¶ ½ µ ¶¾2 µ 1 2 +1 − Γ +1 Var ( X ) = Γ . c c µ
Cuando c es grande, las expresiones anteriores se pueden aproximar mediante µ 2 ¶ π 1 0.57722 0.98905 γ 2 +γ ' 1− + E[ X ] = 1 − + 2 c 6 c 2c c2 Var ( X ) '
0.64493 π2 = 2 6c c2
donde γ es la constante de Euler definida como " γ = l´ımn→∞
# Z µ ¶ ∞ 1 1 1 ∑ k − log(n) = 1 b xc − x ≈ 0.57721566... k =1 n
para la cual b x c se define como
b x c = m´ax {k ∈ Z | k ≤ x } donde Z representa el conjunto de los enteros positivos. En las expresiones anteriores, ´ Gama definida como Γ(.) representa la funcion Γ(t) =
Z ∞ 0
x t−1 e− x dx 10
t > 0.
Sean X1 , X2 , ..., Xn variables aleatorias mutuamente independientes e id´enticamente ´ de distribucion ´ acumulada dada distribuidas tales que Xi ∼ Weibull (c, α, ξ 0 ) con funcion por (3.20), y sean X(1) ≤ X(2) ≤ ... ≤ X(n) las estad´ısticas de orden obtenidas de las n ´ de densidad del estad´ıstico m´as pequeno ˜ X(1) es variables anteriores. La funcion f X(1) ( x ) = n {1 − FX ( x )}n−1 f X ( x ) µ ¶ nc x − ξ 0 c−1 −n{( x−ξ 0 )/α}c e , x > ξ0. = α α
(3.25) (3.26)
´ anterior se sigue f´acilmente que X(1) tiene tambi´en distribucion ´ Weibull, De la ecuacion ´ de densidad de X(r) , para excepto que ahora α es reemplazada por αn−1/c . La funcion 1 ≤ r ≤ n, es f X(r ) ( x ) =
c c n! (1 − e− x )r−1 e− x (n−r+1) cx c−1 , x > 0. (r − 1) ! ( n − r ) !
(3.27)
De (3.27), se obtiene el k-´esimo momento de X(r) h E ( X(r ) )
k
i
n! = (r − 1) ! ( n − r ) ! n! = (r − 1) ! ( n − r ) !
Z ∞ 0 r −1
x
k
n µ
∑ (−1)
i
i =0
µ
=
1−e
k n! Γ 1+ (r − 1) ! ( n − r ) ! c
xc
o r −1
r−1 i
¶ r −1
∑
i =0
¶Z
e− x ∞
0
c ( n −r +1)
e− x
cx c−1 dx
c ( n −r + i +1)
x k cx c−1 dx
(−1)i (r−i 1) (n − r + i + 1)1+(k/c)
ver por ejemplo Johnson et al. (1994).
3.2.2.
Estimacion ´ por M´axima Verosimilitud
´ m´as comun ´ es cuando el umbral ξ 0 se considera conocido pero c y α La situacion ´ se desconocidas, por lo cual deben de ser estimados; dos ejemplos de esta distribucion ˜ n de una distribucion ´ muestran la Figura 3.2. Dada una muestra aleatoria de tamano ´ de distribucion ´ Weibull de dos par´ametros con funcion f X (x) =
c ³ x ´c−1 −( x/α)c e , x > 0, α α 11
(3.28)
0.8
Weibull(c = 1.5 α = 1.25)
0.0
0.2
0.4
fX(x)
0.6
Weibull(c = 2 α = 1.25)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
x
´ Weibull (c = 1.5, α = 1.25) y Weibull (c = 2,α = 1.25) respectivamente Figura 3.2: Distribucion (con ξ 0 = 0)
los estimadores de m´axima verosimilitud cˆ y αˆ satisfacen las ecuaciones ( αˆ = y
( cˆ =
n
∑ Xicˆ log Xi
i =1
1 n cˆ Xi n i∑ =1
)(
n
∑ Xicˆ
i =1
)1/cˆ (3.29)
) −1
n
−
−1
1 log Xi n i∑ =1
.
(3.30)
Si ξ 0 no es igual a cero, entonces cada Xi se remplaza por Xi − ξ 0 en las ecuaciones anteriores. El valor cˆ se obtiene al resolver (3.30) y al usarlo en (3.29) se obtiene αˆ ; estas ´ anal´ıtica. soluciones se obtienen mediante m´etodos num´ericos ya que no existe solucion ˆ entonces αˆ en (3.29) debe ser el estimador Debe notarse que si c fuese conocido igual a c, de m´axima verosimilitud de α. Si el par´ametro ξ 0 es tambi´en desconocido, entonces los estimadores de m´axima verosimil12
ˆ αˆ y ξˆ0 satisfacen las ecuaciones itud c, ( αˆ = ( cˆ =
n
1 ( Xi − ξˆ0 )cˆ ∑ n i =1
)(
n
∑ (Xi − ξˆ0 )cˆ log(Xi − ξˆ0 )
)1/cˆ
) −1
n
∑ (Xi − ξˆ0 )cˆ
i =1
(3.31)
n
−
i =1
−1
1 log( Xi − ξˆ0 ) n i∑ =1
(3.32)
y n
n
i =1
i =1
(cˆ − 1) ∑ ( Xi − ξˆ0 )−1 = cˆαˆ −cˆ ∑ ( Xi − ξˆ0 )cˆ−1
(3.33)
Si el valor ξˆ0 satisface (3.31)-(3.33) y es m´as grande que X(1) , entonces e´ ste es el estimador de m´axima verosimilitud de ξ 0 . Adem´as si el EMV ξˆ0 es tal que ξˆ0 = X(1) , entonces ˆ Johnson et al. (1994) las expresiones en (3.31) y (3.32) son satisfechas al usar αˆ y c.
3.3.
Informacion ´ de Kullback-Leibler o Distancia Entre dos Modelos
Suponga que f y g son dos modelos de probabilidad completamente conocidos, de los cuales f representa la ”realidad o la verdad” y se denotar´a a g como un modelo de aprox´ en t´erminos de una distribucion ´ de probabilidades. De acuerdo con lo anterior, imacion ´ (K-L) entre los modelos f y g como Kullback y Leibler definieron la informacion µ
Z
I ( f , g) =
f ( x ) log
f (x) g( x | θ )
¶ dx
´ la notacion ´ I ( f , g) donde log denota el logaritmo natural. De acuerdo con su definicion, ´ perdida cuando se usa g para aproximar a f ”. Como una interrepresenta ¨ınformacion ´ heur´ıstica, ”I ( f , g) es la distancia entre los modelos g y f ”. pretacion Para escoger entre dos modelos candidatos, un principio consiste en buscar el modelo ´ posible; esto es equivalente a minimizar que pierda la menor cantidad de informacion
13
´ se perI ( f , g) con respecto a g. La densidad f se asume completamente conocida y solo mite variar a g sobre todo el espacio de modelos posibles, indizado por θ. La distancia ´ es g Kullback-Leibler (K-L) es una medida de ineficiencia al suponer que la distribucion ´ es f . cuando la verdadera distribucion La distancia de Kullback-Leibler puede ser conceptualizada como una distancia dirigida entre dos modelos, por decir f y g. Estrictamente hablando es una medida de ”discrepancia”, ya que no es una distancia puesto que la medida de f a g no es la misma que la medida de g a f . ´ A un nivel heur´ıstico, ¨ınformacion”se define como (loge ( f ( x )) para funciones de ´ de Kullback-Leibler es un tipo de densidad de probabilidad continua. La informacion ¨ınformacion ´ cruzada”. El lado derecho es el valor esperado del logaritmo de la razon ´ de dos distribuciones ( f y g) y puede visualizarse como el promedio, con respecto a f , de log( f /g). La distancia de K-L (I ( f , g)) siempre es positiva, excepto cuando las dos distribuciones f y g son id´enticas, es decir, I ( f , g ) = 0 ⇔ f ( x ) = g ( x ), Burnham and Anderson (2002).
3.3.1.
La realidad f , considerada como constante
Todo lo anterior hace pensar que se debe conocer en su totalidad las distribuciones ´ f y g para poder calcular la distancia K-L entre ambos modelos. Sin embargo, si solo se considera la distancia relativa entre ambos modelos, conocer en su totalidad ambos ´ modelos es inutil, de donde I ( f , g) puede escribirse equivalentemente como Z
I ( f , g) =
Z
f ( x ) log( f ( x ))dx −
14
f ( x ) log( g( x | θ ))dx.
´ anterior es una Note que cada uno de los dos t´erminos del lado derecho de la ecuacion esperanza estad´ıstica con respecto a f . Por lo tanto la distancia K-L puede ser expresada como la diferencia de dos esperanzas. I ( f , g) = E f [log( f ( X ))] − E f [log( g( X | θ ))] ´ ´ facilita la derivacion ´ del Criterio cada una de ellas con respecto a f . Esta ultima expresion ´ de Akaike (AIC). de Informacion ´ de la distribucion ´ La esperanza E f [log( f ( X ))] es una constante que depende solo verdadera desconocida, es decir, se desconoce en el an´alisis a f , por lo tanto, tratando este t´ermino desconocido como una constante, es posible calcular una medida dirigida relativa, Burnham and Anderson (2002). Claramente, si se calcula la esperanza E f [log( g( X | θ ))], se puede estimar I ( f , g), excepto por una constante C, concretamente E f [log( f ( X ))],
I ( f , g) = C − E f [log( g( X | θ ))] o I ( f , g) − C = − E f [log( g( X | θ ))]. El t´ermino ( I ( f , g) − C ) es una distancia relativa entre los modelos f y g; entonces, ´ de un E f [log( g( X | θ ))] se convierte en la medida de calidad de inter´es para la seleccion mejor modelo. Para dos modelos g1 y g2 , si I ( f , g1 ) < I ( f , g2 ) entonces g1 es mejor, luego I ( f , g1 ) − C < I ( f , g2 ) − C y por lo tanto
− E f [log( g1 ( X | θ ))] < − E f [log( g2 ( X | θ )). 15
´ Mas aun I ( f , g2 ) − I ( f , g1 ) ≡ − E f [log( g2 ( X | θ ))] + E f [log( g1 ( X | θ ))] y entonces se sabe que g1 es mejor modelo que g2 . Sin conocer C no es posible saber absolutamente qu´e tan buen modelo es g1 , pero es posible identificar el hecho de que es mejor modelo que g2 . T´ıpicamente, se postulan modelos a-priori gi ( x | θ ) y se desea seleccionar el mejor entre ellos para tomarlo como base para el an´alisis e inferencia. Es decir se selecciona el modelo que tenga la menor distancia. Alternativamente, se selecciona al modelo que ´ con respecto a la realidad. El concepto de realpierda la menor cantidad de informacion ´ es idad f , se convierte en una constante y no es necesario conocerla. En la pr´actica, solo posible obtener un estimador relativo de la distancia K-L para cada uno de los modelos ´ gi ( x | θˆ), Burnham and Anderson (2002). de aproximacion
3.4.
Criterio de Informacion ´ de Akaike
´ K-L basado en el logaAkaike (1973) encontro´ una manera de calcular la informacion ´ de verosimilitud en su punto m´aximo. Dado un modelo estructurado ritmo de la funcion ´ para el cual existe un unico valor de θ que minimiza la distancia K-L (I ( f , g)), el minimizador depende de f , de la estructura del modelo g, del espacio de par´ametros, y del espacio muestral (es decir, de la estructura y naturaleza de los datos observados). En este ´ de MV, sea e´ ste valor igual a θ0 . sentido, hay un valor ”verdadero”de θ bajo la estimacion ´ Entonces θ0 es el mejor valor de θ para el modelo g; en efecto, la p´erdida de informacion de K-L es minimizada por θ0 . En el an´alisis de datos los par´ametros del modelo deben estimarse, y hay generalmente incertidumbre sobre estas estimaciones. Los modelos basados en estimaciones de ´ del caso en el que los los par´ametros, es decir θˆ y no θ, representan una gran distincion ´ afecta la forma en que se usa la distancia de par´ametros son conocidos. Esta distincion 16
´ del modelo. La diferencia entre tener θ o θ0 y tener el K-L como base para la seleccion ´ del modelo estimador θˆ es muy importante, y b´asicamente cambia el criterio de seleccion a minimizar la distancia K-L esperada, en vez de minimizar la distancia conocida K-L sobre el conjunto R conformado por los modelos considerados. £ ¤ ´ calcular EY EX log( g( X | θˆ(y))) por medio de la minimizacion ´ Aunque es tentador solo de log(L(θˆ) | datos) para cada gi . Akaike (1973) mostro´ que maximizar el logaritmo de la ´ de verosimilitud genera un estimador sesgado positivamente del modelo objetifuncion vo. Tambi´en encontro´ que bajo ciertas condiciones este sesgo es aproximadamente igual ´ a K, el numero de par´ametros estimables en el modelo aproximado. Este es un resultado ´ asintotico de fundamental importancia [Burnham and Anderson (2002)]: Un estimador aproximadamente insesgado de £ ¤ EY EX log( g( X | θˆ))) para muestras grandes y ”buenos”modelos, es decir, modelos adecuados o que propicien ´ es una buena estimacion log(L(θˆ | datos)) − K. Este resultado es equivalente a log(L(θˆ | datos)) − K = constante − Eˆ θˆ [ I ( f , gˆ )] ´ del sesgo (K) ver donde gˆ = g(· | θ ). Para m´as detalles sobre el t´ermino de correccion Takeuchi (1976). ´ llamado AIC; multiplicando log(L(θˆ | Akaike (1973) definio´ un criterio de informacion y)) − K por -2 obtuvo AIC = −2 log(L(θˆ | y)) + 2k. Por lo tanto, m´as que tener una medida simple de la distancia directa entre dos modelos (es decir la distancia K-L), se tiene en lugar un estimador de la esperanza de la distancia
17
´ relativa entre el modelo candidato y el mecanismo desconocido (quiz´as de dimension infinita) que en realidad genero´ los datos. ´ de verosimilEl t´ermino log(L(θˆ | y)) es un valor num´erico del logaritmo de la funcion itud en su punto m´aximo y este punto m´aximo corresponde a los valores de los esti´ madores de m´axima verosimilitud. El numero de par´ametros estimados es denotado por K. En algunas clases de modelos existen excepcionalmente param´etros que no son estimables a partir de los datos, y no son contados en K. Por ejemplo la no estimabilidad ocurre en el an´alisis de datos generados por conteos cuando una celda no tiene observaciones, por lo tanto un par´ametro no es estimable. En la pr´actica, se calcula el AIC para cada uno de los modelos candidatos y se selecciona el que tenga el menor valor de AIC. Este es el modelo ”m´as cercano” al modelo desconocido que en realidad genero´ los datos, ´ dentro de todos los modelos candidatos considerados en la discriminacion. ´ el AIC Por supuesto, dentro de todos los modelos considerados en la discriminacion seleccionar´a el mejor de ellos dentro de este conjunto, y si todos los modelos son muy ´ al modelo verdadero, el criterio AIC seleccionar´a el pobres respecto a la aproximacion ´ mejor de ellos, pero incluso e´ ste mejor modelo relativamente ser´a pobre; es por esta razon ´ de los modelos a discriminar, con el unico ´ que se debe se hacer una buena seleccion ´ proposito de seleccionar un modelo adecuado. ˜ anadiendo ˜ El criterio I ( f , g) puede hacerse m´as pequeno m´as par´ametros conocidos ´ de en el modelo aproximado g, por lo tanto, para un conjunto de datos fijos, la adicion par´ametros en el modelo gi le permitir´a ser muy cercano a f . Sin embargo, cuando estos par´ametros deben estimarse se agrega incertidumbre al c´alculo de la distancia relativa ´ momento, la adicion ´ de un par´ametro m´as tendr´a el efecto contrario de K-L. En algun ´ la estimacion ´ de la distancia relativa de al deseado (reducir Eθˆ[ I ( f ,g)] ). En esa situacion, K-L se incrementa debido al ¨ruido” en los par´ametros estimados que en realidad no son ´ necesarios para conseguir un buen modelo. Este fenomeno puede ser observado mini-
18
´ mizando el criterio de informacion AIC = −2 log(L(θˆ | y)) + 2K donde el primer t´ermino del lado derecho decrece conforme m´as par´ametros se adicionan al modelo aproximado; mientras que el segundo t´ermino (2K) se incrementa conforme ´ se incrementa el numero de par´ametros al modelo aproximado. Algunos investigadores han considerado que K sea una medida de ´´complejidad´´, pero esto no es necesario, ´ aun cuando no es irracional. Considere a K principalmente como una simple expresion ´ ´ de verosimilitud, como un estimador del sesgo asintotico en el logaritmo de la funcion £ ¤ de EY EX log( g( x | θˆ)(y))) Shibata (1983), Linhart and Zuchini (1986), Bozdogan (1987), y Sakamoto (1991).
3.4.1.
Criterio de Informacion ´ de segundo orden
´ ´ con El AIC puede operar mal si el numero de par´ametros es muy grande en relacion ˜ de la muestra. Sugiura (1978), Sakamoto (1991) derivaron una variante de seel tamano gundo orden del AIC llamado AICc . ˜ que reHurvich y Tsai (1989) estudiaron el comportamiento para muestras pequenas sulto´ en un criterio llamado AICc µ AICc = −2 log(L(θˆ | y)) + 2K
n n−K−1
¶
´ donde el t´ermino 2K es multiplicado por n/(n − K − 1) y se llama factor de correccion. Este puede ser escrito como: 2K (K + 1) AICc = −2 log(L(θˆ | y)) + 2K + n−K−1 o equivalentemente AICc = AIC +
19
2K (K + 1) n−K−1
˜ de muestra. A menos que el tamano ˜ de muestra sea grande con donde n es el tamano ´ respecto al numero de par´ametros estimados, se recomienda el uso del AICc . ´ Si n es grande El criterio AICc simplemente tiene adicionado un t´ermino de correccion. ´ de segundo orden es despreciable y el criterio AIC con respecto a K, entonces la correccion funciona de forma correcta. Bedrick and Tsai (1994) proveen una mejora adicional pero es complicada de calcular (para m´as detalles verHurvich and Tsai (1991) y Hurvich and Tsai (995a) y Hurvich and Tsai (995b), y Hurvich et al. (1990)). En general, se recomienda ´ n/K es pequena ˜ (¡40). Llegar a una decision ´ sobre el uso el uso del AICc cuando la razon ´ del tamano ˜ de la razon ´ n/K. Si la razon ´ n/K es suficientede AIC y AICc depender´a solo mente grande, entonces AIC y AICc son similares y seleccionaran el mismo modelo. Debe usarse AIC o AICc consistentemente en un an´alisis, en vez de mezclar ambos criterios. Pocos paquetes de software proveen valores de AICc , pero e´ stos pueden ser calculados f´acilmente.
3.5.
Prueba de Vuong
´ entre dos modelos no anidados, es decir, La prueba de Vuong considera la eleccion modelos tales que uno no puede ser expresado en t´erminos del otro y viceversa. Sean Fθ y Gγ modelos con densidades f (y | xi , θ ) y g(y | xi , γ). El estad´ıstico LR para discriminar entre el modelo Fθ contra Gγ es ˆ γˆ ) ≡ L f (θˆ) − L g (γˆ ) = LR(θ,
n
f (y | x , θˆ)
∑ log g(y | xii, γˆ )
i =1
donde θˆ y γˆ son los estimadores de m´axima verosimilitud correspondientes a las densidades Fθ y Gγ respectivamente. En el caso especial donde los modelos son anidados, Fθ ⊂ Gγ , se obtiene el resultado ˆ γ) ˆ se distribuye como Ji-cuadrada general que establece que la variables aleatoria 2 LR(θ, ´ bajo la hipotesis nula Gy = Fθ . Considere la prueba para modelos no anidados, Fθ * Gy 20
´ ji-cuadrada no es apropiada, Cameron and Trivedi y Gy * Fθ , entonces la distribucion (1998). ´ al problema mediante la aplicacion ´ del Teorema Cox (1961,1962) propuso la solucion ´ de que Fθ es el modelo verdadero. Este enfoque es Central de L´ımite, bajo la suposicion dif´ıcil de implementar para el an´alisis, ya que se requiere obtener E f [log [ f (y | xi , θ )/g(y | xi , γ)]]; donde, E f denota la esperanza con respecto a la densidad f (y | xi , θ ). Adem´as, si un estad´ıstico similar se obtiene con los papeles de Fθ y Gγ invertidos, es posible encontrar ambos casos, que el modelo Fθ es rechazado en favor de Gγ , y que el modelo Gγ es rechazado en favor de Fθ . Voung (1989) en lugar de discriminar entre modelos usando sus distancias con respecto al proceso real que genero´ los datos, con densidad h0 (y | xi ), propuso el estad´ıstico A B 1 ˆ γˆ )/wˆ 2 = √ LR(θ, n
TLR,NN =
donde 1 n f (y | xi , θˆ) A = √ ∑ ln n i=1 g(y | xi , γˆ ) Ã !2 Ã !2 1 n n 1 f (y | xi , θˆ) f (y | xi , θˆ) ln − ln B= n i∑ g(y | xi , γˆ ) n i∑ g(y | xi , γˆ ) =1 =1 y
à !2 à !2 n n ˆ ˆ 1 1 f ( y | x , θ ) f ( y | x , θ ) i i − ln , wˆ 2 = ∑ ln ˆ) n i =1 g(y | xi , γˆ ) n i∑ g ( y | x i, γ =1
´ de la varianza de es una estimacion
ˆ γˆ ).Para √1 LR ( θ, n
modelos estrictamente no anidados
´ a una N (0, 1) Cameron and Trivedi (1998), es decir, TLR,NN converge en distribucion TLR,NN d N (0, 1) − → bajo
· H0 : Eh
¸ f ( y | xi , θ ) =0 log g ( y | xi , γ ) 21
donde Eh denota la esperanza con respecto a h0 (y | xi ). Por lo tanto, se rechaza la hipotesis nula de equivalencia de modelos en favor de Fθ con un nivel de significancia α, si ´ TLR,NN > zα (o si TLR,NN < −zα ). La hipotesis nula no se rechaza si | TLR,NN |≤ z α2 , Cameron and Trivedi (1998). Donde zα es el cuantil superior α de la densidad N (0, 1), es decir, el valor zα , tal que P(z ≥ zα ) = α.
3.6.
M´etodo de Simulacion ´ de Monte Carlo
El m´etodo de Monte Carlo es un m´etodo num´erico que permite resolver, en forma ´ de variables aleatorias, problemas matem´aticos. El aproximada y mediante la simulacion m´etodo fue bautizado as´ı por su analog´ıa con los juegos de ruleta de los casinos, el m´as ´ fue propuesta en 1856 por c´elebre de los cuales es el de Monte Carlo, cuya construccion ´ el pr´ıncipe Carlos III de Monaco. La importancia actual del m´etodo Monte Carlo se basa en la existencia de problemas ´ por m´etodos exclusivamente anal´ıticos o num´ericos, pero que que tienen dif´ıcil solucion pueden asociarse a un modelo probabil´ıstico artificial, el cual puede estudiarse en forma sencilla simulando variables aleatorias correspondientes al mismo. La idea fundamental ´ de Monte Carlo para realizar inferencia, es que las caracter´ısticas detr´as de la simulacion ´ de un estad´ıstico pueden ser observadas generando muestras aleatorias de la poblacion de inter´es repetidamente y observando el comportamiento del estad´ıstico sobre todas las ´ del estad´ıstico tomando muestras simuladas. En otras palabras, calcular la distribucion ´ y guardando los valores del estad´ıstico para cada una muestras al azar de la poblacion de las muestras. Los valores observados del estad´ıstico para cada muestra se usan para ´ de e´ ste. calcular la distribucion ´ de Monte Carlo consiste en definir una El primer paso para realizar la simulacion ´ que debe suponerse representa a la poblacion ´ leg´ıtima en todos sus pseudo poblacion aspectos relevantes. El uso de la palabra pseudo enfatiza el hecho de que se obtienen 22
´ muestras mediante el uso de una computadora y de numeros pseudo aleatorios. ´ debe ser tal que sea posible tomar muestras usando la comLa pseudo poblacion putadora. El procedimiento b´asico de Monte Carlo para estudiar propiedades de un estimador, consta de los siguientes pasos ´ o modelo que representa la verdadera poblacion ´ 1. Determinar la pseudo poblacion de inter´es. ´ 2. Usar un m´etodo de muestreo para muestrear la pseudo poblacion. 3. Calcular el valor del estad´ıstico y almacenarlo. 4. Repetir los pasos 2 y 3 M ocasiones. 5. Usar los M valores obtenidos en el paso 4 para estudiar emp´ıricamente la distribu´ del estad´ıstico. cion ´ el analista Es importante tener en mente que cuando se muestrea la pseudo poblacion, ´ estad´ıstica. debe de asegurar que todas las caracter´ıstica relevantes reflejen la situacion ˜ de muestra y estrategia de muestreo se deben de usar Por ejemplo, el mismo tamano cuando se trata de comprender el rendimiento de una estad´ıstica. Esto quiere decir que la ´ de la estad´ıstica obtenida v´ıa simulacion ´ de Monte Carlo es eficaz solamente distribucion para las condiciones del procedimiento de muestreo y las suposiciones sobre la pseudo ´ Wendy and Angel (2002). poblacion,
23
Cap´ıtulo 4 Metodolog´ıa ´ entre las distribuSe definieron los objetivos para el desarrollo de la discriminacion ´ de literatura para las ciones Weibull y Gausiana Inversa, y se llevo´ a cabo una revision distribuciones Weibull y Gausiana Inversa consideradas en el presente estudio, detallando las funciones de densidad, momentos, estimadores y propiedades sobre variables aleatorias para cada una de las distribuciones. ´ bibliogr´afica sobre criterios de discriminacion ´ enAs´ı mismo, se realizo´ una revision tre distribuciones candidatas a modelar cierto proceso que genero´ una serie de obser´ vaciones, as´ı como sus diferentes aplicaciones; de igual manera se realizo´ una revision ´ de estad´ısticas, como lo es la simulacion ´ de bibliogr´afica de las t´ecnicas de validacion ´ se procedio´ a la aplicacion ´ de los criterios Monte Carlo. Recabada toda esta informacion ´ bibliogr´afica a una serie de observaciones reales provenientes analizados en la revision de una muestra aleatoria formada por mediciones de di´ametros normales de a´ rboles de pino. ´ bibliogr´afica se procedio´ a un estudio de simulacion ´ de Monte Despu´es de la revision Carlo y al an´alisis de un conjunto de datos reales para ejemplificar el uso de las t´ecnicas. Para efectuar las simulaciones de Monte Carlo, el proceso consistio´ en simular observaciones de un conjunto de poblaciones completamente conocidas; a e´ stas observaciones 24
se les aplicaron los estad´ısticos incluidos en el an´alisis para cada una de las poblaciones ´ se observo´ el comporsimuladas, y posteriormente, con el resultado de la simulacion, tamiento de los estad´ısticos (AIC y Vuong). ´ de los criterios descritos, se considero´ un conjunto de di´amePara mostrar la aplicacion tros normales conformado por n = 67 observaciones el cual se muestra en el Cuadro 4.1. Las observaciones, de di´ametros normales (medidos a 1.3 m de altura y expresadas en ´ metros) corresponden al mismo numero de a´ rboles seleccionados de manera aleatoria. El a´ rea en la cual fueron medidos e´ stos a´ rboles se ubica dentro del Parque Nacional Zoquia´ Experimental de Ensenanza ˜ ´ de Zoquiapan, y es conocida como Estacion e Investigacion ´ experimental se localiza en los l´ımites del Estado de M´exico y Puebla, en pan. La estacion ´ montanosa ˜ la region conocida como Sierra Nevada, aproximadamente entre los paralelos 19◦ 12’30” y 19◦ 20’00” de latitud Norte, y entre los meridianos 98◦ 42’30” y 98◦ 30’00” de longitud Oeste, en la cual predomina bosque boreal y en el subtipo bosque de pino, los arboles predominantes del bosque boreal son el oyamel (Abies sp.), ocote (Pinus spp.), aile (Alnus spp.), pinabete (Pseudotsuga spp.) y enebros alpinos (Juniperus spp.), el subpiso esta constituido por zacatones (Festuca, Muhlenbergia y Agrostis) y se encuentra confinada en la zona alta de la cordillera volc´anica que cruza al pa´ıs, ver Jos´e et al. (1981).
25
´ forestal experiCuadro 4.1: Di´ametros normales de n = 67 a´ rboles de pino en la estacion mental Zoquiapan. Obs.
Di´ametro
Obs.
Di´ametro
(m)
Obs.
(m)
Di´ametro (m)
1
0.22
24
0.19
47
0.235
2
0.25
25
0.17
48
0.385
3
0.13
26
0.23
49
0.36
4
0.13
27
0.13
50
0.175
5
0.3
28
0.13
51
0.22
6
0.21
29
0.65
52
0.49
7
0.185
30
0.17
53
0.36
8
0.155
31
0.205
54
0.2
9
0.14
32
0.13
55
0.16
10
0.385
33
0.165
56
0.58
11
0.105
34
0.14
57
0.21
12
0.93
35
0.13
58
0.2
13
0.245
36
0.275
59
0.38
14
0.17
37
0.6
60
0.4
15
0.225
38
0.6
61
0.17
16
0.23
39
0.18
62
0.2
17
0.41
40
0.18
63
0.13
18
0.375
41
0.17
64
0.26
19
0.245
42
0.21
65
0.59
20
0.2
43
0.16
66
0.3
21
0.15
44
0.225
67
0.12
22
0.19
45
0.135
23
0.135
46
0.25
26
Cap´ıtulo 5 Aplicacion, ´ Simulacion ´ y Discusion ´ 5.1.
Aplicacion ´
´ de los criterios de discriminacion ´ en una aplicacion ´ forestal real Para una aplicacion considere los datos contenidos en el Cuadro 4.1. En la Figura 5.1 se observa el compor-
4 0
2
Frecuencia
6
8
tamiento de tales di´ametros.
0.0
0.2
0.4
0.6
0.8
1.0
X
Figura 5.1: Histograma de n = 67 observaciones de di´ametros normales 27
´ del histograma presenta cierta asimetr´ıa y una Como es posible notar, la distribucion forma semejante a la que las distribuciones Weibull y Gausiana Inversa pueden presentar, ´ y por el hecho de que con ciertos valores particulares de los par´ametros. Es por esta razon la variable aleatoria que representa Di´ametro es estrictamente positiva que los modelos Weibull y Gausiano inverso se consideran buenos modelos candidatos.
5.1.1.
Estimacion ´ Distribucion ´ Weibull
´ de la distribucion ´ Weibull se realiza mediante el procedimiento de La estimacion m´axima verosimilitud. Considere la muestra aleatoria X1 , X2 , ..., Xn donde Xi ∼ Weibull (c, α), ´ de densidad como (3.28). La funcion ´ de verosimilitud para de la muestra ancon funcion terior es
`(θ ) = L( f X ( x )) =
³ c ´n α
Ã
n
∏ i =1
xi α
! c −1
Ã
!
n
exp − ∑ ( xi /α)c i =1
donde θ = (c, α) ´ de verosimilitud se obtiene Tomando el logaritmo natural de la funcion
`(θ ) = log( L( f X ( x ))) = n log
³c´ α
n
+ (c − 1) ∑ log i =1
³x ´ i
α
n
−∑
i =1
³ x ´c i
α
´ anterior se deriva con respecto a c y α, y las expresiones resultantes se Si la expresion ´ de e´ stas generan los estimadores de m´axima verosimilitud de igual a cero, la solucion ´ de c y α, como se muestra en (3.29) y (3.30). Para maximizar el logaritmo de la funcion verosimilitud se utilizo´ el procedimiento nlm() del software R, el cual minimiza funciones continuas, diferenciables (hasta de orden dos) y con un solo punto critico, por lo ´ de verosimilitud se proporciona − log(`(θ )) y as´ı de esta que para optimizar la funcion ´ manera el procedimiento nlm() maximiza tal funcion. ´ a optimizar El procedimiento nlm() necesita que se le especifique e ingrese la funcion y los valores de inicio de los par´ametros de la siguiente manera, nlm( f , p) 28
´ a minimizar (en este caso se maximizo´ la forma donde el t´ermino f representa la funcion ´ Para ya mencionada) y el t´ermino p, los valores de inicio para realizar la optimizacion. ´ de la funcion ´ de verosimilitud de la siguiente manera esto se realizo´ la programacion > lfmv<-function(p,X) + { +
cc<-p[1]
+
a<-p[2]
+
f<-n*(log(cc/a))+(cc-1)*sum(log(X/a))-sum((X/a)**cc)
+
return(-f)
+ } > ´ en la donde el t´ermino p es un vector que contiene ambos par´ametros de la distribucion, ´ se encuentra el par´ametro c y en la segunda el par´ametro α; el objeto X primer posicion ´ Finalmente representan los valores de la muestra con los que se realizar´a la optimizacion. ´ el termino return(− f ) proporciona el valor de la funcion. ´ a maximizar y los valores de inicio iguales a (0.2, 0.1) para c y α Dada la funcion ´ nlm() de la siguiente forma respectivamente, se evaluo´ la funcion > estimacion<-nlm(lfmv,c(0.2,0.1)) donde el objeto estimacion contiene los par´ametros ajustados por m´axima verosimilitud, adem´as de otros componentes que el software proporciona. Las estimaciones obtenidas se extraen del objeto estimacion agregando al final del objeto el signo $ y la palabra estimate, de la siguiente manera, >estimacion$estimate, con este procedimiento se obtuvo las siguientes estimaciones basada en la muestra observada cˆ = 1.8291474
29
y αˆ = 0.2925577. ´ de densidad estimada Weibull se De acuerdo con estos valores estimados la funcion
2.5
3.0
muestra en la Figura 5.2
1.5 0.0
0.5
1.0
fX(x)
2.0
Weibull(c = 1.829147 α = 0.292557)
0.0
0.2
0.4
0.6
0.8
1.0
x
´ estimada Weibull (cˆ = 1.8291474, αˆ = 0.292557) Figura 5.2: Distribucion
´ con el comportamiento de la muestra se tiene la Figura 5.3 que conEn comparacion ´ del modelo Weibull, con el comportamiento de la muestra trasta la estimacion
30
8 6 4
fX(x)
0
2
Weibull(c = 1.829147 α = 0.292557)
0.0
0.2
0.4
0.6
0.8
1.0
x
´ estimada Weibull (cˆ = 1.8291474, αˆ = 0.292557) vs histograma muestral Figura 5.3: Distribucion
´ realizada por el modelo Weibull aparentemente Se puede observar que la estimacion resulta malo. Lo cual sera determinado por el criterio AICc y Vuong m´as adelante.
5.1.2.
Estimacion ´ Distribucion ´ Gausiana Inversa
´ de la distribucion ´ Gausiana Inversa basada en la muestra observada se La estimacion ´ Weibull, es decir, por m´axima verosimilobtuvo de igual forma que para la distribucion ´ de itud. Considere la muestra aleatoria X1 , X2 , ..., Xn , donde Xi ∼ IG (µ, λ), con funcion ´ de verosimilitud est´a dada densidad como en (3.2), entonces la correspondiente funcion por
µ L( f X ( x )) =
λ 2π
¶n/2 Ã
n
∏ i =1
1 xi3
!1/2
"
λ exp − 2 2µ
n
( x i − µ )2 ∑ xi i =1
#
´ con respecto a µ y λ respectivamente se obtiene dos expreDerivando la expresion siones que al ser igualadas a cero y resolver para µ y λ generan los estimadores de m´axiˆ ma verosimilitud µˆ y λ. 31
∑n X µˆ = X¯ = i=1 i n y λˆ =
1 n−1 ∑in=1 ( Xi−1
− X¯ −1 )
.
Sustituyendo valores para la muestra (di´ametros normales) se obtiene µˆ = 0.2576866 y λˆ = 0.9933674 en consecuencia, la gr´afica de la densidad estimada Gausiana Inversa se muestra en la
3
4
figura (5.4)
2 0
1
fX(x)
IG(µ = 0.2576866 λ = 0.9933674)
0.0
0.2
0.4
0.6
0.8
1.0
x
´ estimada IG (µˆ = 0.2576866, λˆ = 0.9933674) Figura 5.4: Distribucion
32
5.1.3.
Criterio de Informacion ´ AICc
´ en el presente an´alisis har´a uso de criterio AIC modificado por el La discriminacion ˜ de muestra (llamada AICc ) y se define como tamano 2K (K + 1) AICc = −2 log(L(θˆ) | y) + 2K + n−K−1 ´ de verosimilen donde el t´ermino log(L(θˆ) | y) representa el punto m´aximo de la funcion itud evaluada en los estimadores de m´axima verosimilitud. El valor m´aximo del logar´ de verosimilitud, as´ı como los par´ametros estimados se obtienen al itmo de la funcion ejecutar el procedimiento nlm() del software R. El valor correspondiente del AICc para ´ Gausiana Inversa y el correspondiente codigo ´ la distribucion en R se muestra en el Anexo ´ Weibull es AICc = −78.7921. 1. El valor de AICc para la distribucion Cuadro 5.1: Valores estimados de AICc para cada una de las densidades Densidad
AICc
Estimaciones
Estimaciones
Weibull
-78.7921
cˆ = 1.8291474
αˆ = 0.292557
Gausiana I.
-104.0795
µˆ = 0.2576866
λˆ = 0.9933674
5.1.4.
Estad´ıstico de Vuong
´ ´ de modelos es contrastar El proposito b´asico del estad´ıstico de Vuong en la seleccion ´ la hipotesis
· H0 : Eh
¸ f ( y | xi , θ ) =0 ln g ( y | xi , γ )
´ ´ Gausiana Inversa Para llevar a cabo la prueba de hipotesis denotaremos a la distribucion ´ como f (y | xi , θ ), y al modelo Weibull como g(y | xi , γ). Para contrastar la hipotesis se calculo´ el valor del estad´ıstico TLR,NN = 22.49689, y el valor del percentil z0.05 = 1.644854; de ´ se rechazo´ la hipotesis ´ acuerdo con la regla de decision nula, de igualdad de modelos, en favor del modelo Fθ , que para este caso representa el modelo determinado por la distribu´ Gausiana Inversa con par´ametros estimados µˆ = 0.2576866 y λˆ = 0.9933675 como cion 33
´ Weibull con par´ametros cˆ = 1.8291474 y αˆ = 0.292557 mejor modelo, sobre la distribucion ´ ´ realizada por el critedenotada por g(y | xi , γ). Notese que al igual que la discriminacion ´ estimada rio de AICc , el criterio de Vuong selecciono´ como mejor modelo la distribucion ´ Gausiana Inversa. El codigo propio al software R con el cual se realizo´ la prueba se presenta en el Anexo 1.
5.1.5.
Resultados
De acuerdo con el criterio de AICc , el cual establece que el mejor modelo para aproximarse al que en realidad genero´ los datos, es aquel que tenga el menor valor de AICc , de acuerdo con esto y con nuestros resultados concentrados en el Cuadro 5.1, el mejor ´ Gausiana Inversa con par´ametros µˆ = 0.2576866 y modelo resulta ser la distribucion λˆ = 0.9933674. Al igual que el criterio del AICc el estad´ıstico de Vuong selecciono´ como mejor modelo ´ del modelo Gausiano al Gausiana inverso. Como se muestra en la Figura 5.5, la eleccion ´ con la estimacion ´ de la distribucion ´ Weibull, es inverso como el mejor, en comparacion ´ realizada por la distribucion ´ Gausiana Injustificado gr´aficamente, es decir, la estimacion versa resulta ser semejante al comportamiento de las frecuencias observadas en la muestra.
34
8 6 4
fX(x)
0
2
IG(µ = 0.2576866 λ = 0.9933674)
0.0
0.2
0.4
0.6
0.8
1.0
x
´ estimada IG (µˆ = 0.2576866, λˆ = 0.9933674) vs histograma muestral Figura 5.5: Distribucion
´ considere la prueba de bondad de ajuste Con la finalidad de validar la discriminacion, ´ ´ de la de Kolmogorov-Smirnov, que establece, bajo la hipotesis nula, que la distribucion muestra es F0 . ´ La hipotesis que se prueba es: ´ de distribucion ´ de donde la muestra aleatoria proviene es F0 H0 : La funcion ´ de distribucion ´ de donde la muestra aleatoria proviene no es F0 H1 : La funcion El procedimiento es el siguiente; 1) A partir de los valores observados x1 , x2 , ..., xn se calcula
dn = m´ax | Fˆn ( xi ) − F0 ( xi ) | donde 1 n Fˆn ( xi ) = ∑ I(−∞,xi ] ( xi ) n i =1 35
´ de distribucion ´ especificada en la hipotesis ´ y F0 es la funcion nula. 2) Se rechaza H0 con un nivel de significancia α, si dn > k1−α , donde k1−α es tal que se cumple n √ 2 2 α ≈ 1 − H ( nk1−α ) = 1 − 2 ∑ (−1)i−1 e−2ni k1−α i =1
Considerando F0 = IG (µˆ = 0.2576866, λˆ = 0.9933674) y la muestra conformada por los n = 67 di´ametros normales, se tiene que dn = m´ax | Fˆn ( xi ) − F0 ( xi ) |= 0.1445098. Los valores cr´ıticos para k1−α , cuando n > 50 y α = 0.05, est´an determinados por
1.360 √ , n
por lo que k1−α = 0.1856975 ´ de donde, como k1−α > dn , entonces no se rechaza la hipotesis nula, es decir, el modelo estimado Gausiano inverso ajusta satisfactoriamente al conjunto de observaciones diametricas.
5.2.
Simulacion ´ y Evaluacion ´ del Comportamiento de los Estad´ısticos
5.2.1.
Simulacion ´ del criterio AICc
´ ´ es evaluar el comportamiento de cualquier estad´ıstico; El proposito de la simulacion en ese sentido, se realizaron una serie de simulaciones, en las cuales se supone una dis´ completamente conocida. Se llevo´ a cabo un proceso de muestreo simulado a tribucion partir del modelo conocido y posteriormente se realizo´ el calculo de los valores del AICc para ambas distribuciones (Gausiana Inversa y Weibull). ˜ 10,20,30,50 y 60 y para En el presente estudio se simularon 10,000 muestras de tamano ˜ de muestra se calcularon los 10,000 valores del estad´ıstico AICc . Adem´as cada tamano 36
´ la simulacion ´ se consideraron tres distribuciones Gausiana Inversa para la realizacion completamente conocidas, IG (µ = 0.5, λ = 1),IG (µ = 0.5, λ = 0.8) y IG (µ = 0.5, λ = ´ ya descrito. 2.5), para las cuales se realizo´ el proceso de simulacion ´ Gausiana Inversa y se Es natural pensar que si se simulan muestras de una poblacion ´ el criterio siempre elecalculan los valores de AICc para las distribuciones en cuestion, ´ gir´a el modelo Gausiano inverso como el mejor modelo, comparando con la distribucion Weibull. En este sentido despu´es de generar las 10,000 muestras para cada uno de los 5 ˜ de muestra, se calculo´ la fraccion ´ de rechazos observado, es decir, la fraccion ´ tamanos de veces en las que el criterio AICc no eligio´ a el modelo Gausiano inverso como mejor modelo, siendo este el que en realidad genero´ las observaciones. Para el caso en el que se generaron las muestras con el modelo IG (µ = 0.5, λ = 1) los resultados se concentran en el Cuadro 5.2. Cuadro 5.2: Fracciones de rechazo observadas. Modelo verdadero IG (µ = 0.5, λ = 1), y ´ basada en el criterio AICc con nivel de significancia nominal de α = 0.05 discriminacion ˜ de Muestra Tamano
´ de rechazo observada Fraccion
10
0.2779
20
0.1679
30
0.1310
40
0.0642
60
0.0463
´ en la cual el modelo verdadero es el modelo determinado por IG (µ = Para la situacion 0.5, λ = 0.8), las fracciones de rechazo observadas se concentran en el Cuadro 5.3.
37
Cuadro 5.3: Fracciones de rechazo observadas. Modelo verdadero IG (µ = 0.5, λ = 0.8), ´ basada en el criterio AICc con nivel de significancia α = 0.05 y discriminacion ˜ de Muestra Tamano
´ de rechazo observada Fraccion
10
0.0702
20
0.0115
30
0.0051
40
0.0013
60
0.0037
Finalmente, para el caso en el que el modelo verdadero est´a determinado por la dis´ IG (µ = 0.5, λ = 2.5), las fracciones de rechazo se presentan en el Cuadro 5.4 tribucion Cuadro 5.4: Fracciones de rechazo observadas. Modelo verdadero IG (µ = 0.5, λ = 2.5), ´ basada en el criterio AICc con nivel de significancia α = 0.05 y discriminacion ˜ de Muestra Tamano
´ de rechazo observada Fraccion
10
0.0621
20
0.0100
30
0.0044
40
0.0093
60
0.0045
´ est´a determinada por Tambi´en se considero´ el caso en el que la pseudo poblacion ´ se consideraron los siguientes tres modelos, una densidad Weibull, para esta situacion Weibull (c = 1.5, α = 0.5),Weibull (c = 1.5, α = 0.95) y Weibull (c = 1.5, α = 2.5), para los ´ de rechazo observada, es decir, la fraccion ´ cuales se calculo´ su correspondiente fraccion de veces (de las 10,000 muestras) en las que el criterio AICc no eligio´ al modelo Weibull ´ Gausiana Inversa, siendo que e´ ste es el que en como mejor modelo sobre la distribucion ´ est´a deterrealidad genero´ las observaciones. Para el caso en el que la pseudo poblacion 38
minada por el modelo Weibull (c = 1.5, α = 0.5), las fracciones de rechazo se muestran en el Cuadro 5.5. Cuadro 5.5: Fracciones de rechazo observadas. Modelo verdadero Weibull (c = 1.5, α = ´ basada en el criterio AICc con nivel de significancia de α = 0.05 0.5), y discriminacion ˜ de Muestra Tamano
´ de rechazo observada Fraccion
10
0.3544
20
0.2269
30
0.1576
40
0.0744
60
0.0519
Para el caso en el que las muestras fueron generadas bajo el modelo Weibull (c = 1.5, α = 0.95), las fracciones de rechazo se muestran en el Cuadro 5.6. Cuadro 5.6: Fracciones de rechazo observadas. Modelo verdadero Weibull (c = 1.5, α = ´ basada en el criterio AICc con nivel de significancia α = 0.05 0.95), y discriminacion ˜ de Muestra Tamano
´ de rechazo observada Fraccion
10
0.3532
20
0.2218
30
0.1487
40
0.0771
60
0.0520
Finalmente, para el caso en el que se considero el modelo determinado por Weibull (c = 1.5, α = 2.5), las fracciones de rechazo se muestran en el Cuadro 5.7.
39
Cuadro 5.7: Fracciones de rechazo observadas. Modelo verdadero Weibull (c = 1.5, α = ´ basada en el criterio AICc con nivel de significancia α = 0.05 2.5) y discriminacion ˜ de Muestra Tamano
´ de rechazo observada Fraccion
10
0.3515
20
0.2260
30
0.1506
40
0.0728
60
0.0532
Con la finalidad de observar el comportamiento del criterio AICc en el caso en el que ´ est´a determinada por un modelo diferente a las distribuciones inla pseudo poblacion volucradas en el an´alisis (Weibull y Gausiana Inversa), se realizaron tres simulaciones ´ Gamma de dos par´ametros. Para estas simulaciones se consideraron bajo la distribucion ´ Gamma con par´ametro tres casos, el primer caso est´a determinado por una distribucion de forma igual a 3 y de escala a 0.9, el segundo caso est´a determinado por una distribu´ Gamma con par´ametro de forma igual a 1.5 y de escala a 0.9; finalmente, el tercer cion ´ Gamma con par´ametro de forma 2 y escala 0.9. caso consiste en una distribucion ˜ de muestra se Para cada uno de los tres casos ya mencionados y para cada tamano ´ de seleccion, ´ es decir, en la que la distribucion ´ Gausiana Inversa fue calculo´ la fraccion elegida como mejor modelo sobre el modelo Weibull, basada en el criterio AICc , dichas fracciones se concentran en el Cuadro 5.8. ´ Cuadro 5.8: Fracciones de Seleccion ˜ Tamano de
10
20
30
50
60
Caso 1
0.4881
0.4402
0.4008
0.3594
0.3402
Caso 2
0.3996
0.2839
0.2084
0.1273
0.0956
Caso 3
0.4430
0.3703
0.3033
0.2103
0.1822
Muestra
40
´ para la fraccion ´ de seleccion ´ es, por ejemplo, para el caso numero ´ Una interpretacion ˜ de muestra 30, aproximadamente el 40.08 % de las 10000 simulaciones uno y el tamano fueron discriminadas en favor del modelo Gausiano inverso; de manera semejante se puede interpretar cada una de las fracciones. Las instrucciones propias al lenguaje R con las que se realizo´ el c´alculo de los valores ´ as´ı como para de las fracciones de rechazo observadas y las fracciones de discriminacion, ´ se presenta en el Anexo 2. general la simulacion
5.2.2.
Simulacion ´ del criterio de Vuong
Al igual que con el estad´ıstico AICc , con el de Vuong se realizo´ un proceso de simu´ por el m´etodo de Monte Carlo, para lo cual se consideraron dos casos; el primero lacion ´ consiste en probar el criterio de Vuong considerando la hipotesis. · H0 : Eh
¸ f ( y | xi , θ ) =0 ln g ( y | xi , γ )
que contrasta la equivalencia de modelos. El correspondiente estad´ıstico de contraste es 1 ˆ γˆ )/wˆ 2 TLR,NN = √ LR(θ, n en donde ˆ γˆ ) = L f (θˆ) − L g (γˆ ) LR(θ, es la diferencia de las funciones de verosimilitud en su punto m´aximo, y à !2 à !2 n n ˆ ˆ f ( y | x , θ ) 1 f ( y | x , θ ) 1 i i − log wˆ 2 = ∑ log ˆ) n i =1 g(y | xi , γˆ ) n i∑ g ( y | x i, γ =1 para el cual el modelo representado por f (y | xi , θ ) es el modelo descrito por la dis´ Gausiana Inversa, mientras que el modelo representado por g(y | xi , γ) es la tribucion ´ Weibull. distribucion
41
Para el an´alisis del comportamiento del estad´ıstico de Vuong, bajo este primer caso ´ consistio´ en generar datos a partir se realizaron tres simulaciones, la primer simulacion ´ IG (µ = 0.8, λ = 1.5) con la que se genero´ la muestra; la segunda de una distribucion ´ se realizo´ bajo el modelo Weibull (c = 4, α = 1.5) y finalmente una tercera simulacion ´ se genero´ bajo una distribucion ´ Ji-cuadrada central con 15 grados de libertad, simulacion esto con la finalidad de analizar el comportamiento del criterio de Vuong bajo un modelo ´ plantea que no es alguno de los dos involucrados en el estudio. Como la regla de decision ´ rechazar al nivel de significancia de 0.05 la hipotesis nula de equivalencia de modelos, en favor de Fθ , es decir, es mejor modelo que Gγ si TLR,NN > z0.05 , es decir, mayor al cauntil ´ N (0, 1), entonces al realizar el muestreo de una poblacion ´ superior z0.05 de la distribucion determinada por el modelo IG (µ = 0.8, λ = 1.5), las fracciones para las cuales el criterio ´ de Vuong rechazo´ la hipotesis nula se esperan sean muy cercanas a uno, e´ stas fracciones de resumen en el Cuadro 5.9. Cuadro 5.9: Fracciones de rechazo, modelo verdadero IG (µ = 0.8, λ = 1.5), criterio de ´ de Vuong con nivel de significancia α = 0.05 discriminacion ˜ de Muestra Tamano
´ de rechazo Fraccion
10
0.7057
20
0.8417
30
0.9035
40
0.9616
60
0.9762
´ en la que la muestra proviene de una poblacion ´ determinada por Para la situacion Weibull (c = 4, α = 1.5), las fracciones de rechazo se presentan en el Cuadro 5.10.
42
Cuadro 5.10: Fracciones de rechazo, modelo verdadero Weibull (c = 4, α = 1.5), criterio ´ de Vuong con nivel de significancia α = 0.05 de discriminacion ˜ de Muestra Tamano
´ de rechazo Fraccion
10
0.7283
20
0.7620
30
0.7819
40
0.8062
60
0.8150
´ en la que la muestra proviene Finalmente, las fracciones de rechazo, para la situacion ´ simulada determinada por una distribucion ´ Ji-cuadrada central con 15 de una poblacion grados de libertad, se presentan en el Cuadro 5.11. ´ Cuadro 5.11: Fracciones de rechazo, modelo verdadero es χ215 , criterio de discriminacion de Vuong con nivel de significancia α = 0.05 ˜ de Muestra Tamano
´ de rechazo Fraccion
10
0.6991
20
0.6944
30
0.6618
40
0.6156
60
0.5963
El segundo caso planteado del criterio de Vuong fue realizado bajo el cambio de papeles de los modelos, es decir, ahora el modelo determinado por f (y | xi , θ ) paso´ a ser descrito por el modelo Weibull, y el modelo representado por g(y | xi , γ) se convirtio´ en la densidad Gausiana Inversa. Para este caso las poblaciones simuladas fueron las mismas que para el primer caso, este cambio de papeles de las distribuciones lo permite el criterio de Vuong y este otro enfoque puede ser otra forma de trabajar. Las fracciones de 43
´ en la que la poblacion ´ est´a determinada por la distribucion ´ rechazo para las situacion IG (µ = 0.8, λ = 1.5) se presentan en el cuadro 5.12. Cuadro 5.12: Fracciones de rechazo, modelo verdadero IG (µ = 0.8, λ = 1.5), criterio de ´ de Vuong con nivel de significancia α = 0.05 discriminacion ˜ de Muestra Tamano
Fracciones de rechazo
10
0.7225
20
0.7682
30
0.7929
40
0.8005
60
0.8081
´ en la que la poblacion ´ esta determinada por el modelo Weibull (c = Ahora para la situacion 4, α = 1.5), las fracciones de rechazo se presentan en el Cuadro 5.13. Cuadro 5.13: Fracciones de rechazo, modelo verdadero Weibull (c = 4, α = 1.5), criterio ´ Vuong de discriminacion ˜ de Muestra Tamano
Fracciones de rechazo
10
0.7118
20
0.8430
30
0.9084
40
0.9011
60
0.9116
´ est´a determinada Finalmente, las fracciones de rechazo cuando la pseudo poblacion ´ Ji-cuadrada central con 15 g.l. se presentan en el Cuadro 5.14. por una distribucion
44
´ Cuadro 5.14: Fracciones de rechazo, modelo verdadero χ215 , criterio de discriminacion Vuong ˜ de Muestra Tamano
Fracciones de rechazo
10
0.3030
20
0.2993
30
0.3369
40
0.3749
60
0.3945
Como se observa, el comportamiento del estad´ıstico de Vuong para ambos casos es semejante, es decir, cuando se espera que las proporciones sean grandes (cercanas a uno) ˜ (cercanas a cero) lo son, e´ stas lo son, y cuando las proporciones se esperan sean pequenas ´ similarmente al caso salvo ciertas situaciones que se comentaran en la siguiente seccion, ´ es tomada de un modelo distinto a los involucrados en que la muestra de la simulacion en el an´alisis. Cabe mencionar que para el segundo caso del criterio de Vuong, e´ ste discriminar´a ahora en favor del modelo Weibull. Las l´ıneas que ejecutan y calculan cada una de los simulaciones anteriores son presentados en el Anexo 2.
5.3.
Discusion ´
´ de par´ametEs importante tener presente la complejidad que el proceso de estimacion ros implica, en este sentido, computacionalmente es m´as sencillo obtener las estimaciones ´ con el modelo Weibull, debido a propias al modelo Gausiano nverso en comparacion ´ Gausiana Inversa son cerradas, que las expresiones correspondientes a la distribucion es decir, pueden ser calculadas anal´ıticamente, en contraste con las estimaciones correspondientes al modelo Weibull que deben de ser obtenidas mediante un proceso de 45
´ num´erica. Es precisamente en este punto en el cual algunas estimaciones optimizacion ´ bibliogr´afica, y el proceso de estipueden no existir, como ya se planteo´ en la revision ´ resulte aun ´ m´as complicado. Por estas razones es que desde el punto de vista macion ´ de par´ametros propia al modelo Gausiano invercomputacional y pr´actico, la estimacion so resulta m´as sencilla de realizar. Como se puede notar en los cuadros de resumen de las fracciones de rechazo obser´ que realiza el criterio AICc es cada vez mejor, cuando el tamano ˜ vadas, la discriminacion de muestra es mayor, es decir, el modelo estimado ajusta mejor a los datos observados ´ se basa en un numero ´ cuando la estimacion mayor de observaciones. Por ejemplo, con´ simulada est´a determinada por el modelo sidere el Cuadro 5.3 para el cual la poblacion ˜ de muestra n = 20, y el nivel de significancia observado GI (µ = 0.5, λ = 0.8) y el tamano es menor (0.0115) que el nivel de significancia nominal (α = 0.05) con el que el criterio ´ AICc realiza la discriminacion. ´ fue realizada bajo el modelo Weibull (c = 1.5, α = Para el caso en el que la simulacion 0.5), el nivel de significancia observado m´as cercano al nivel de significancia nominal ˜ de muestra n = 60 con un valor del criterio AICc (0.05) se observo cuando el tamano de 0.0519, es decir, el modelo Weibull requiere de mayor cantidad de datos para poder discriminar efectivamente en favor de tal modelo cuando se simula del mismo. ´ se pudo notar que el comportamiento del criterio de Durante el proceso de simulacion ´ (AICc ) es como se espera, es decir, cuando se simulan datos provenientes discriminacion ´ Gausiana Inversa esperamos que el criterio discrimine siempre en de una distribucion ´ es este el comportamiento. De igual manera ocurrio´ con este sentido, y como se observo, ´ Weibull. el proceso cuando se simularon datos que proven´ıan de una distribucion ´ del criterio de discriminacion ´ de Vuong, para el primer caRespecto a la simulacion ´ so, en el cual se rechaza en favor del modelo Gausiano inverso, y en que la poblacion ˜ de muestra es n = 40, est´a determinada por el modelo IG (µ = 0.8, λ = 1.5) y el tamano ´ de rechazo, es decir, las veces que el estad´ıstico de Vuong rechazo´ la igualdad la fraccion 46
de modelos en favor del modelo Gausiano inverso, fue de 0.9616, es decir, el nivel de significancia observado es 0.0384, que es menor que el nivel de significancia nominal de ´ ´ de 0.05 con el cual se realiza la prueba de hipotesis propia al criterio de discriminacion Vuong. ´ para el segundo caso considerado en la simulacion ´ de Vuong, para el En comparacion, ´ en la cual la poblacion ´ cual se rechaza en favor del modelo Weibull, considere la situacion ˜ de muestra n = 60, la simulada proviene del modelo Weibull (c = 4, α = 1.5) y tamano ´ de rechazo fue de 0.9116, es decir, un nivel de significancia estimado de 0.0884, el fraccion ´ cual es mayor al nivel de significancia nominal de la prueba de Vuong. En comparacion ´ del criterio de Vuong, se puede notar que con el primer caso considerado en la simulacion el modelo Gausiano inverso en contraste con el modelo Weibull requiere menor cantidad de observaciones para presentar un nivel de significancia observado menor que el nivel de significancia nominal (0.05) propio a la prueba de Vuong.
47
Cap´ıtulo 6 Conclusiones y Recomendaciones ´ de Vuong y AICc , Del presente trabajo se concluye que los criterios de discriminacion puede considerarse como criterios simples de usar para discriminar entre modelos candidatos a describir cierto conjunto de datos. Los resultados de las simulaciones de Monte Carlo muestran que el comportamiento de los criterios de Vuong y AICc es congruente con lo que se espera cuando hay una dis´ efectiva: cuando se opera bajo cualquiera de los modelos, Weibull o Gausiano criminacion ´ la discriminacion ´ en un alto porcentaje es realizada en inverso, como pseudo poblacion, favor del que genera la muestra. ´ entre los modelos Weibull y GauComo se puede notar al realizar la discriminacion siano inverso mediante los criterios de Vuong y AICc , ambos discriminaron en favor de ´ Gausiana inversa como mejor modelo para describir el comportamiento la distribucion del conjunto de observaciones consideradas, es decir las bservaciones diam´etricas. Esto es respaldado por la prueba de Kolmogorov-Smirnov, mediante la cual se concluye que el modelo seleccionado por los criterios de Vuong y AICc ajusta satisfactoriamente al conjunto de observaciones. ´ bibliogr´afica, el criterio AICc corregido por el tamano ˜ de De acuerdo con la revision ´ n/K < 40, y el criterio de informacion ´ no muestra, es decir, para el caso en que fraccion 48
corregido AIC, al igual que la prueba de Vuong, representan buenos criterios de discrimi´ para seleccionar entre modelos candidatos a ajustar un conjunto de observaciones nacion diam´etricas. Del presente trabajo se concluye tambi´en que debe tenerse en cuenta que el uso del ˜ de la muestra, es decir, si la fraccion ´ n/K < 40; criterio AICc o AIC, est´a sujeto al tamano ´ que es recomendable el uso alternativo de ambos criterios cuando la es por esta razon ´ del tamano ˜ de la muestra as´ı lo requiera. situacion La prueba de Vuong puede ser utilizada en ambos sentidos para poder concluir si cierto modelo es el mejor dentro de los candidatos a modelar nuestros datos; es decir, se puede probar la igualdad de modelos en el sentido de que el mejor modelo sea el ´ establecido en la hipotesis nula y viceversa. ´ a los niveles de significancia observados, la tendencia de la fracciones de En relacion rechazo es mejor con el modelo Gausiano inverso que con el Weibull; es decir, para re´ entre ambos modelos basados en el criterio de AICc , la estimacion ´ alizar la discriminacion ˜ de muestra n = 20, esdel modelo Gausiano inverso se realiza mejor a partir del tamano ˜ de muestra menor el modelo estimado ajusta mejor a los to significa que con un tamano ´ se recomienda que el tamano ˜ de muestra sea mayor o datos observados. Por esta razon, ´ que realiza el modelo Gausiano inverso es cada vez igual a 20, puesto que la estimacion ˜ de muestra, y es cuando el nivel de significancia observado mejor a partir de este tamano comienza a ser menor que el nivel de significancia nominal. Considerando el criterio de Vuong, para el caso en el que se rechaza en favor del modelo Gausiano inverso, el nivel de significancia observado es menor que el nominal ˜ de muestra n = 40, en contraste con el modelo Weibull, para el cual a partir del tamano ´ con un tamano ˜ de muestra n = 60, el nivel de significancia observado es mayor aun al nominal; es decir, el modelo Weibull requiere una mayor cantidad de observaciones que el modelo Gausiano inverso para que el modelo estimado ajuste mejor a los datos observados. 49
´ el modelo Gausiano inverso requiere menor cantidad de observaciones En conclusion, ´ en comparacion ´ con el modelo Weibull. En este senpara realizar una buena estimacion ˜ de muestra para estimar un modelo Gausiano inverso tido, se recomienda que el tamano ˜ de muestra el nivel de signifisea mayor o igual a 40, puesto que a partir de este tamano cancia observado para los criterios AICc y Vuong es menor que el nominal. ´ del modelo Weibull debe Respecto a la estimabilidad del los modelos, la estimacion ´ num´erica, en contraste con el modeefectuarse mediante un proceso de optimizacion lo Gausiano inverso, en el cual las estimaciones provienen de ecuaciones anal´ıticas, es decir, expresiones en forma cerrada, por lo que las estimaciones siempre existen. Por lo ´ el modelo Gausiano inverso es tanto, desde el punto de vista de facilidad de estimacion, m´as f´acil de estimar que el Weibull; por estas razones, adem´as de que el modelo Gausiano inverso puede adoptar una gran diversidad de formas, bajo ciertos valores de sus par´ametros, y requiere menor cantidad de observaciones para realizar una mejor esti´ se recomienda el uso de e´ ste para modelar estructuras diam´etricas. macion,
50
´ APENDICES
51
Ap´endice A Codigo ´ R usado para calcular el criterio AICc Instrucciones propias del software R con las cuales se calculo´ el criterio de Informa´ de Akaike corregido por el tamano ˜ de muestra AICc cion datos<-read.table("C:/----/) X<-datos[,4] n<-length(X) lfmvig<-function(~ n) { mu<-~ n[1] lambda<-~ n[2] ff<-(n/2)*log(lambda/(2*pi))+(1/2)*sum(log(1/(X**3))) -(lambda/(2*mu**2))*sum(((X-mu)**2)/X) return(-ff) } uno<-c(rep(1/n,n)) MU<- uno%*%X 52
lambdaa<-n/(sum(X^{-1}-MU^{-1})) parametros_i<-c(MU,lambdaa) ln_i<-lfmvig(c(MU,lambdaa)) AIC_i<--2*ln_i+2*k*(n/(n-k-1)) ´ Gausiana Inversa no es neceEs importante notar que para el caso de la distribucion ´ nlm(), puesto que las ecuaciones para estimar los par´ametros sario utilizar la funcion ´ por m´axima verosimilitud son anal´ıticas (cerradas). El valor de AICc para la distribucion Gausiana Inversa es AICc = −104.0795. De manera similar se obtuvo el valor de AICc ´ Weibull, a continuacion ´ se presenta el codigo ´ para la distribucion correspondiente lfmvw<-function(p) { cc<-p[1] a<-p[2] f<-n*(log(cc/a))+(cc-1)*sum(log(X/a))-sum((X/a)**cc) return(-f) } parametros_w<-nlm(lfmvw,c(0.2,0.1))$estimate k<-length(parametros_w) ln_w<--nlm(lfmvw,c(0.2,0.1))$minimum AIC_w<--2*ln_w+2*k*(n/(n-k-1)) ´ se presentan las instrucciones para el c´alculo de la prueba de Vuong, A continuacion en el lenguaje R.
#
M´ aximo del logaritmo Gausiana Inversa
X<-datos 53
n<-length(X) lfmvig<-function(~ n) { mu<-~ n[1] lambda<-~ n[2] ff<-(n/2)*log(lambda/(2*pi))+(1/2)*sum(log(1/(X**3)))(lambda/(2*mu**2))*sum(((X-mu)**2)/X) return(ff) } uno<-c(rep(1/n,n)) MU<- uno%*%X lambdaa<-n/(sum(X^{-1}-MU^{-1})) parametros_i<-c(MU,lambdaa) ln_i<-lfmvig(c(MU,lambdaa)) #-----------------------------------------------------------# #
Maximo del logaritmo Weibull
#
lfmvw<-function(p) { cc<-p[1] a<-p[2] f<-n*(log(cc/a))+(cc-1)*sum(log(X/a))-sum((X/a)**cc) return(-f) } parametros_w<-nlm(lfmvw,c(0.2,0.1))$estimate k<-length(parametros_w) ln_w<--nlm(lfmvw,c(0.2,0.1))$minimum #----------------------------------------------------------# #
W estimada^2
# 54
w_2<-(1/n)*sum((ln_i-ln_w)^2)-((1/n)*(ln_i-ln_w))^2 Tlrnn<-((1/(sqrt(n)))*(ln_i-ln_w))/(w_2)
55
Ap´endice B Codigo ´ R usado para la simulacion ´ de Monte Carlo ´ Codigos propios al lenguaje R con los cuales fueron realizadas cada una de las simulaciones tanto para el criterio AICc como para la prueba de Vuong. ´ ´ del criterio AICc , bajo el muestreo de una El siguiente codigo realiza la simulacion ´ determinada por la densidad Gausiana Inversa con los siguientes par´ametros poblacion µ = 0.5 y λ = 1, µ = 0.5 y λ = 0.8, µ = 0.5 y λ = 2.5 respectivamente, se puede realizar ´ de par´ametros solo ´ haciendo un cambio en la cada uno de las pruebas para combinacion l´ınea X<-rinverse.gaussian(tm,0.5,1) ˜ que es la que genera las muestras aleatorias, en la cual el t´ermino tm representa el tamano de muestra, el siguiente objeto representa el valor del par´ametro µ y el tercer t´ermino ´ Gausiana Inversa. representa el valor del par´ametro λ de la distribucion #----------------------------------------------------------# #
Simulaci´ on Gausiana Inversa 56
#----------------------------------------------------------# #
Genera Muestras Aleatorias de Inversa Gausiana
rinverse.gaussian<-function(n,mu,lambda) { k<-0 for (i in 1:n) { d<-0 y<-rchisq(1,1) u<-runif(1,0,1) r1<-(mu/2*lambda)*(2*lambda+mu*y-sqrt(4*lambda*mu*y+mu**2 *y**2 )) r2<-((mu**2)/r1) if (u<(mu/(mu+r1))) d<-r1 else d<-r2 k[i]<-d } return(k) } #--------------------------------------------------------------#
Simulaci´ on criterio de Vuong
M<-10000 res<-matrix(0,M,10) l<-0 for (i in c(1:5)) { pos<-c(10,20,30,50,60) tm<-pos[i] l<-l+2 57
for (j in c(1:M)) { X<-rinverse.gaussian(tm,0.8,1.5) n<-length(X) #--------------------------------------------------------lfmvig<-function(~ n) { mu<-~ n[1] lambda<-~ n[2] ff<-(n/2)*log(lambda/(2*pi))+(1/2)*sum(log(1/(X**3))) -(lambda/(2*mu**2))*sum(((X-mu)**2)/X) ff2<-((1/2)*log(lambda/(2*pi))+(1/2)*(log(1/(X**3))) -(lambda/(2*mu**2))*(((X-mu)**2)/X)) return(cbind(ff,ff2)) } #-------------------------------------------------------uno<-c(rep(1/n,tm)) MU<- uno%*%X lambdaa<-n/(sum(X^{-1}-MU^{-1})) parametros_i<-c(MU,lambdaa) ln_i<-lfmvig(parametros_i)[1]
#-------------------------------------------------------------#
Maximo del logaritmo Weibull lfmvw<-function(p,X) { cc<-p[1] a<-p[2] 58
f<-n*(log(cc/a))+(cc-1)*sum(log(X/a))-sum((X/a)**cc) return(-f) } parametros_w<-nlm(lfmvw,c(0.2,0.1),X=X)$estimate k<-length(parametros_w) ln_w<--nlm(lfmvw,c(0.2,0.1),X=X)$minimum #--------------------------------------------------------------#
W estimada^2 fwe<--(apply(matrix(X,nc=1),1,lfmvw,p=parametros_w)+(n-1)*
(log(parametros_w[1]/parametros_w[2]))) w_2<-(1/n)*sum((fwe-lfmvig(parametros_i)[,2])^2)-((1/n)*(ln_w-ln_i))^2 Tlrnn<-1/sqrt(n)*(ln_w-ln_i)/w_2 res[j,c(l-1,l)]<-cbind(abs(Tlrnn),qnorm(0.975)) } } ´ efectuan ´ las siguientes l´ıneas solo el c´alculo de los niveles de significancia observados ´ de valores de AIC para cada una de las muestras para mediante una simple comparacion ˜ de muestra, adem´as del calculo de fraccion. ´ los correspondientes tamanos z<-0 nob<-0 for (i in c(1:5)) { z<-z+2 f<-0 for (j in c(1:M)) { if (res[j,z-1]>res[j,z]) f<-f+1 59
else f<-f+0 } nob[i]<-(f/M) } ´ de AICc cuando la muestra Para el resto de las simulaciones, es decir, para la simulacion ´ determinada por una distribucion ´ Weibull el cambio en el proviene de una poblacion ´ ´ es en la l´ınea en la que se genera la muestra aleatoria, en lugar de presentar codigo solo la l´ınea X<-rinverse.gaussian(tm,0.5,1) el cambio se realiza por X<-rweibull(tm,0.5,1) ˜ tm con par´ametros c = 0.5 y α = 1, al igual la cual genera muestras aleatorias de tamano ´ de la distribucion ´ Gamma, el cambio solo ´ es por la l´ınea que para la simulacion X<-rgamma(tm,shape=2,scale=0.9) ´ del criterio de Vuong bajo el proceso de Monte Carlo, en Para el proceso de simulacion ´ en la que el modelo determinado por f (yi | xi , θ ) es la distribucion ´ Gausiana la situacion ´ Weibull y simulando una Inversa y el modelo nombrado por g(yi | xi , γ) la distribucion ´ Gausiana Inversa el codigo ´ poblacion es el siguiente: #------------------------------------------------------------#
Genera Muestras Aleatorias de Gausiana Inversa
rinverse.gaussian<-function(n,mu,lambda) { k<-0 60
for (i in 1:n) { d<-0 y<-rchisq(1,1) u<-runif(1,0,1) r1<-(mu/2*lambda)*(2*lambda+mu*y-sqrt(4*lambda*mu*y+mu**2 *y**2 )) r2<-((mu**2)/r1) if (u<(mu/(mu+r1))) d<-r1 else d<-r2 k[i]<-d } return(k) } #----------------------------------------------------------#
Simulacion criterio de Vuong
M<-1000 res<-matrix(0,M,10) l<-0 for (i in c(1:5)) { pos<-c(10,20,30,50,60) tm<-pos[i] l<-l+2 for (j in c(1:M)) { X<-rweibull(tm,3.5,5.5) n<-length(X) #-----------------------------------------------------61
#
Maximo de la Verosimilitud IG
lfmvig<-function(~ n) { mu<-~ n[1] lambda<-~ n[2] ff<-(n/2)*log(lambda/(2*pi))+(1/2)*sum(log(1/(X**3))) -(lambda/(2*mu**2))*sum(((X-mu)**2)/X) return(ff) } uno<-c(rep(1/n,tm)) MU<- uno%*%X lambdaa<-n/(sum(X^{-1}-MU^{-1})) parametros_i<-c(MU,lambdaa) ln_i<-lfmvig(c(MU,lambdaa)) #-----------------------------------------------------#
Maximo del logaritmo Weibull
lfmvw<-function(p) { cc<-p[1] a<-p[2] f<-n*(log(cc/a))+(cc-1)*sum(log(X/a))-sum((X/a)**cc) return(-f) } parametros_w<-nlm(lfmvw,c(0.2,0.1))$estimate k<-length(parametros_w) ln_w<--nlm(lfmvw,c(0.2,0.1))$minimum #-----------------------------------------------------#
W estimada^2 62
w_2<-(1/n)*sum((ln_i-ln_w)^2)-((1/n)*sum(ln_i+ln_w))^2 Tlrnn<-((1/(sqrt(n)))*(ln_i-ln_w))/(w_2) res[j,c(l-1,l)]<-cbind(Tlrnn,qnorm(0.05,0,1)) } } ´ ´ ˜ esta ultima parte del codigo calcula la promociona para cada uno de los cinco tamanos de muestra en los cuales el criterio de Vuong discrimino correctamente cuando el proceso ´ se realizo bajo un modelo Gausiano Inversa. de simulacion z<-0 nob<-0 for (i in c(1:5)) { z<-z+2 f<-0 for (j in c(1:M)) { if (res[j,z-1]>res[j,z]) f<-f+1 else f<-f+0 } nob[i]<-(f/M) } para las dem´as simulaciones, es decir, para las cuales las muestran provienen de pobla´ se realiza el cambio en la l´ınea ciones Weibull y Ji-Cuadrada solo X<-rinverse.gaussian(tm,3.5,5.5) ´ Weibull la l´ınea para las correspondientes, en concreto cuando se simulo´ una poblacion ´ anterior del codigo cambia por 63
X<-rweibull(tm,3.5,5.5) ´ Ji-cuadrada no central y en el caso cuando se generaron las muestras de una poblacion ´ con 15 grados de libertad el codigo para las muestras cambia por X<-rchisq(tm,15) con los respectivos cambios en los valores de par´ametros para los cuales se realizo´ el ´ proceso de simulacion. ´ del criterio de Vuong, es decir, Para el segundo caso considerado en la simulacion cuando los papeles de las distribuciones determinadas por f (yi | xi , θ ) y g(yi | xi , γ) son ´ Weibull es el siguiente: intercambiados y cuando la muestra proviene de una poblacion M<-10000 res<-matrix(0,M,10) l<-0 for (i in c(1:5)) { pos<-c(10,20,30,50,60) tm<-pos[i] l<-l+2 for (j in c(1:M)) { X<-rweibull(tm,4,1.5) n<-length(X) #------------------------------------------------------------#
Maximo de la verosimilitud IG
lfmvig<-function(~ n) { 64
mu<-~ n[1] lambda<-~ n[2] ff<-(n/2)*log(lambda/(2*pi))+(1/2)*sum(log(1/(X**3)))(lambda/(2*mu**2))*sum(((X-mu)**2)/X) return(ff) } #------------------------------------------------------------uno<-c(rep(1/n,tm)) MU<- uno%*%X lambdaa<-n/(sum(X^{-1}-MU^{-1})) parametros_i<-c(MU,lambdaa) ln_i<-lfmvig(c(MU,lambdaa)) #-----------------------------------------------------------#
Maximo de la verosimilitud
Weibull
lfmvw<-function(p) { cc<-p[1] a<-p[2] f<-n*(log(cc/a))+(cc-1)*sum(log(X/a))-sum((X/a)**cc) return(-f) } parametros_w<-nlm(lfmvw,c(0.2,0.1))$estimate k<-length(parametros_w) ln_w<--nlm(lfmvw,c(0.2,0.1))$minimum #-----------------------------------------------------------#
W estimada^2
w_2<-(1/n)*sum((ln_w-ln_i)^2)-((1/n)*sum(ln_w-ln_i))^2 Tlrnn<-((1/(sqrt(n)))*(ln_w-ln_i))/(w_2) 65
res[j,c(l-1,l)]<-cbind(Tlrnn,qnorm(0.05,0,1)) } } z<-0 nob<-0 for (i in c(1:5)) { z<-z+2 f<-0 for (j in c(1:M)) { if (res[j,z-1]>res[j,z]) f<-f+1 else f<-f+0 } nob[i]<-(f/M) } ´ del criterio de Vuong, la parte final calcula Igual que para el primer caso de simulacion ´ de casos en los que el modelo Weibull es seleccionado como mejor con la proporcion respecto al modelo Gausiano inverso; de igual forma que para el primer caso cuando se ´ generadora de la muestra el cambio en desea probar el estad´ıstico con distinta poblacion ´ ´ se efectua ´ en la l´ınea el codigo solo X<-rweibull(tm,4,1.5) de manera similar que en el primer caso.
66
Referencias Akaike, H. (1973). Information theory as an extension of the maximum likelihood principle. Akademiai Kiado, Budapest. Bain, T. D. R. and Antle, C. E. (1969). Inferences on the parameters of the weibull distribution. Technometrics, 11:445–460. Bedrick, E. J. and Tsai, C.-L. (1994). Model selection for multivariate regression in small samples. Biometrics, 50. Bertram, H., Thomas, W. B., and John, A. K. (2003). Forest Mensuration. John Wiley and Sons Inc., New York. Bozdogan, H. (1987). Model selection and akaike information criterion (aic): The general theory and its analytical extensions. Psychometrika, 52. Burnham, K. and Anderson, D. R. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach, volume Second Edition. Springer. Cameron, A. C. and Trivedi, P. K. (1998). Regression Analysis of Count Data. Cambridge University Press. Cox, D. (1961). Test of separate families of hypotheses. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, volume I. New York. Dobson, A. (2001). An Introduction to Generalized Linear Models. Ed. Chapman and Hall.
67
Fahrmeir, L. and Tutz, W. (2001). Multivariate statistical Modelling Based on Generalized Linear Models. Ed. Springer Verlag, New York. Hurvich, C. M., Shumway, R., and Tsai, C.-L. (1990). Improved estimators of kullbackleibler information for autoregressive model selection in small samples. Biometrika, 77. Hurvich, C. M. and Tsai, C.-L. (1991). Bias of the corrected aic criterion for underfitted regression and time series models. Biometrika, 78. Hurvich, C. M. and Tsai, C.-L. (1995a). Relative rates of convergence for efficient model selection criteria in linear regression. Biometrika, 82. Hurvich, C. M. and Tsai, C.-L. (1995b). Model selection for extended quasi-likelihood models in small samples. Biometrics, 51. Johnson, N. L., Kotz, S., and Balakrishna, N. (1994). Continuous Univariate Distribution, volume 1, Second Edition. John Wiley and Sons Inc. Jorgensen, B. (1982). Statistical Properties of the Generalized Inverse Gaussian Distribution, Lecture Notes in Statistics 9. Ed. Springer Verlag, New York. Jos´e, M., Rafael, P., Ana, S., and Segundo, B. (1981). Ecolog´ıa de la estaci´on experimental ´ zoquiapan. Universidad Autonoma Chapingo, M´exico. Kalbfleisch, J. D. and Prentice, R. L. (2002). The Statistical Analysis of Failure Time Data. Ed. John Wiley and Sons, New York. Klein, J. P. and Moeschberger, M. L. (2003). Survival Analysis: Techniques for censured and Truncates Data. Ed. Springer Verlag. Lawless, J. F. (1982). Statistical Models and Methods for Lifetime Data. John Wiley and Sons, New York. Linhart, H. and Zuchini, W. (1986). Model Selection. John Wiley and Sons, New York. 68
Roy, L. K. and Wasan, M. T. (1968). The first passage time distribution of brownian motion with positive drift. Mathematical Biosciences, 3. Sakamoto, Y. (1991).
Categorical data analysis by AIC.
KTK Scientific Publishers,
Tokyo,Japan. Seshadri, V. (1993). The Inverse Gaussian Distribution: A Case Study in Exponential Families. Oxford Science Publications. Shibata, R. (1983). A theoretical view of the use of AIC, Time series analysis: theory and practice. Elsevier Science Publication, North-Holland. Smith, R. L. and Weissman, I. (1985). Maximum likehood estimation of the lower tail of a probability distribution. Biometrika, 47:285–298. Sugiura, N. (1978). Further analysis of the data by akaike´s information criterion and the finite corrections. Communications in Statistics, Theory and Methods, A7. Takeuchi, K. (1976). Distribution of informational statistics and a criterion of model fitting. Suri-Kagaku. Tweedie, M. C. K. (1957). Statistics properties of inverse gaussian distributions. Annals of Mathematical Statistical, 28. Vuong, Q. H. (1989). Likelihood ratio tests for model selection and non-nested hypothesis. Econometrica, 57. Wendy, L. M. and Angel, R. M. (2002). Computational Statistics Handbook with MATLAB. CHAPMAN AND HALL/CRC.
69