Modelos de distribuciones obtenidas de la distribución normal

MODELOS DE DISTRIBUCIONES OBTENIDAS DE LA DISTRIBUCIÓN NORMAL María F. Martínez de Paz MODELOS DE DISTRIBUCIONES OBTENIDAS DE LA DISTRIBUCIÓN NORMAL María F. Martínez de Paz Memoria presentada como parte de los requisitos para la obtención del título de Grado en Matemáticas por la Universidad de Sevilla. Tutorizada por Inmaculada Barranco Chamorro Índice general Resumen 1 English Abstract 3 1. Herramientas básicas 5 1.1. Conceptos básicos en distribuciones de tiempos de vida . . . . . . . . 5 1.2. Inferencia en muestras completas . . . . . . . . . . . . . . . . . . . . 9 1.3. Inferencia en muestras censuradas . . . . . . . . . . . . . . . . . . . . 12 1.4. DistribuciónNormal............................ 14 2. Distribución Log Normal 21 2.1. Propiedadesbásicas............................ 22 2.2. Inferencia de la distribución Log Normal . . . . . . . . . . . . . . . . 25 2.3. Aplicaciones................................ 29 3. Distribución Inversa Gaussiana 37 3.1. Propiedadesbásicas............................ 39 3.2. Inferencia de la distribución Inversa Gaussiana . . . . . . . . . . . . . 44 3.3. Aplicación ................................. 47 ii modelos de distribuciones obtenidas de la distribución normal 4. Distribución Birnbaum-Saunders 51 4.1. Introducción................................ 51 4.2. Interpretación física . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.3. Propiedadesbásicas............................ 53 4.3.1. Características de las distribuciones de tiempos de vida . . . . 56 4.4. Momentos y momentos inversos . . . . . . . . . . . . . . . . . . . . . 58 4.5. Relación con la distribución Inversa Gaussiana . . . . . . . . . . . . . 63 4.6. Inferencia.................................. 64 4.6.1. Método de los momentos . . . . . . . . . . . . . . . . . . . . . 64 4.6.2. Método de los momentos modificado . . . . . . . . . . . . . . 65 4.6.3. Método de máxima verosimilitud . . . . . . . . . . . . . . . . 66 4.6.4. Máxima verosimilitud para muestras censuradas . . . . . . . 68 4.7. Aplicaciones................................ 69 Bibliografía 77 Resumen El objetivo de este trabajo es el estudio de modelos continuos de tiempo de vida relacionados con la distribución Normal. Específicamente, nos centramos en las distribuciones Log-Normal, Inversa Gaussiana y Birnbaum-Saunders. Todos ellos son modelos biparamétricos en que se darán propiedades que permiten estudiar la forma de la distribución, obtener características de interés en teoría de fiabilidad y análisis de supervivencia, y realizar inferencia. Se ha estudiado la estimación de los parámetros del modelo, tanto para datos completos como para muestras censuradas. Cabe destacar que en la mayoría de las ocasiones, para el proceso de estimación habrá que recurrir a técnicas numéricas. El estudio teórico de cada modelo se ha completado con aplicaciones a datos reales (completos y censurados) utilizando el software RStudio. La memoria se ha estructurado en cuatro capítulos. En el Capítulo 1 se introducen herramientas básicas en modelos de tiempo de vida, como son la función de supervivencia, función de riesgo o hazard, y cuantiles. Para la estimación de parámetros se propone el método de máxima verosimilitud, introduciendo las modificaciones adecuadas para tratar con muestras censuradas. El capítulo se completa con una recopilación de resultados en la distribución normal que se usarán posteriormente. El Capítulo 2 se dedica a la distribución Log-Normal. Se obtienen su función de supervivencia, función de hazard y momentos, los cuáles nos serán útiles para estudiar la forma de la distribución y propiedades de los estimadores. Destacar que en la estimación para datos completos se tienen fórmulas cerradas, y para datos censurados se han de aplicar métodos numéricos (paquete survival de R). En el Capítulo 3 se estudia la distribución Inversa Gaussiana siguiendo la parametrización propuesta por Tweedie. Se recogen propiedades, función generatriz de momentos y de cumulantes y resultados sobre la interpretación de los parámetros. Cabe destacar que al realizar inferencia para muestras completas se tiene un resultado análogo al Teorema de Fisher para poblaciones normales. 2modelos de distribuciones obtenidas de la distribución normal Finalmente, en el Capítulo 4 se realizará un estudio más extenso de la distribución Birnbaum-Saunders, debido a su importante base física entre otras propiedades interesantes. Veremos cómo surge la distribución a partir de un proceso de fatiga cíclica, y las propiedades fundamentales de la misma. Asímismo, se considerará la notable relación que guarda con la distribución Inversa Gaussiana, además de resultados de inferencia mediante diferentes métodos que se verán reflejados en las aplicaciones a datos reales, las cuales darán por concluido este trabajo. English Abstract The aim of this dissertation is to study continuous lifetime models related to the Normal distribution. Specifically, we will focus on the Log-Normal, Gaussian Inverse and Birnbaum Saunders distribution. All of them are two-parameter models in which properties will be given that allow us to study the shape of the distribution, to obtain summaries of interest in Reliability Theory and Survival Analysis, and to carry out Inference. The issue of estimation of parameters in the model is studied for complete data and for censored data. It is important to note that, quite often, numerical techniques are necessary in the process of estimation. Applications to real dataset (complete and censored data) ha been carried out by using the software RStudio. The document is divided into four chapters. In Chapter 1, the basic tools in lifetime models are introduced. These are: survival function, hazard function and quantiles. Estimation of parameters is carried out by using maximum likelihood. Specific results for censored samples are given. To conclude the chapter, a section is included dealing with relevant results in the normal distribution, which will be used later. Chapter 2 is devoted to Log-Normal distribution. Its survival function, hazard function and moments are given. Moments will be useful to study the shape of the distribution and to obtain properties of estimators. As for the estimation of parameters, we highlight that, in the case of complete data we have closed expressions. On the other hand, for censored data numerical techniques must be applied (R survival package). In Chapter 3, the Gaussian Inverse distribution is studied with the parameterization proposed by Tweedie. Properties, moment and cumulant generating functions, and results, which allow us to interpret the parameters in this model, are given. As for the inference in this model for complete data, we highlight a result analogous to Fisher Theorem in normal populations. 10 modelos de distribuciones obtenidas de la distribución normal Método. El método para hallar el estimador de 𝜃consistirá en hallar aquel valor  𝜃= 𝜃(𝑇1,…, 𝑇𝑛)tal que maximiza la función de verosimilitud, es decir, que verifica 𝐿(𝑡𝑖,…, 𝑡𝑛; 𝜃)=m ax 𝜃∈Θ 𝐿(𝑡;𝜃). Este valor recibirá el nombre de estimador máximo verosímil oestimador de máxima verosimilitud del parámetro 𝜃. Sin embargo, es común que la función de verosimilitud en (1.4) tenga una expresión complicada, dificultando los cálculos. Esta función es positiva y sus máximos se mantienen al aplicar un logaritmo neperiano, al ser esta función monótona creciente. Por ello, es frecuente considerar su logaritmo neperiano. Esta nueva función recibe el nombre de función de log-verosimilitud, cuya expresión es 𝑙(𝜃) = log 𝐿(𝑡;𝜃) = 𝑛 ∑ 𝑖=1 log 𝑓(𝑡𝑖;𝜃),∀𝜃∈ℝ.(1.5) Generalmente, la ecuación o sistema de ecuaciones de verosimilitud se puede resolver sin grandes inconvenientes. Sin embargo, en ocasiones hay que recurrir a métodos numéricos iterativos para poder hallar el estimador. Caso uniparamétrico. En el caso de que tengamos un solo parámetro 𝜃 ⊆ Θ⊆ℝ y se verifiquen las siguientes condiciones de regularidad: El campo de variación de 𝜃es un intervalo abierto del eje real. El campo de variación de la variable aleatoria poblacional no depende de 𝜃. La fdd 𝑓(𝑡𝑖, 𝜃)es positiva y derivable respecto de 𝜃. Se verifica la condición del máximo 𝜕2𝑙 𝜕𝜃2|||𝜃= 𝜃<0. Entonces, se cumplen las siguientes propiedades. 1. Consistencia. Los EMV son consistentes, lo que significa que ∀𝜖 > 0, se verifica l ım 𝑛→∞𝑃(| 𝜃−𝜃|< 𝜖)=1,∀𝜃∈ Θ. 2. Insesgadez. En general, los EMV no serán insesgados. En caso de serlo, serán asintóticamente insesgados. 1. herramientas básicas 11 3. Normalidad asintótica. Al aumentar el tamaño muestral, la distribución de  𝜃 puede aproximarse por una distribución normal. Para un parámetro 𝜃∈ Θ, se verifica  𝜃 →𝑁(𝜃, √𝑉 𝑎𝑟( 𝜃), con 𝑉 𝑎𝑟( 𝜃)la varianza asintótica de  𝜃. Caso multiparamétrico 𝑘≥2.En el caso de que la función de densidad dependa de 𝑘parámetros, los EMV se obtienen resolviendo el sistema de ecuaciones de verosimilitud en 𝜃1,…, 𝜃𝑘. Las ecuaciones de este sistema vendrán dadas por 𝜕𝑙(𝑡1,…, 𝑡𝑛;𝜃1,…, 𝜃𝑘) 𝜕𝜃𝑗 = 𝑛 ∑ 𝑖=1 𝜕log 𝑓(𝑡𝑖;𝜃1,…, 𝜃𝑘) 𝜕𝜃𝑗 = 0,∀𝑗= 1,…, 𝑘. En nuestro caso, será 𝑘= 2. Bajo condiciones de regularidad los EMVs obtenidos tendrán buenas propiedades asintóticas, que son una generalización al caso multiparamétrico de las expuestas anteriormente. En particular, será de gran interés la matriz de información de Fisher, pues nos permite obtener la matriz de varianza asintótica de los EMVs. Matriz de información de Fisher En este trabajo consideramos modelos biparamétricos para 𝑇. Dada una muestra aleatoria de T, sus parámetros 𝜃= (𝜃1, 𝜃2)se pueden estimar por máxima verosimilitud. En aquellas situaciones en que se pueda suponer la normalidad asintótica de los estimadores MV, nos interesa obtener la matriz de información de Fisher, pues a partir de ella se pueden proponer aproximaciones asintóticas para las varianzas de  𝜃1y 𝜃2, así como para la covarianza 𝐶𝑜𝑣( 𝜃1, 𝜃2). Denotando por (𝜃)a la matriz de información de Fisher, se tiene que el elemento (𝑖, 𝑗)de esta matriz viene dado por −𝜕𝑙(𝜃) 𝜕𝜃𝑖𝜕𝜃𝑗 , 𝑖, 𝑗 = 1,2, donde 𝑙(𝜃)es la función de logverosimilitud de la muestra en cuestión, dada por (1.5). (𝜃)depende de los valores de los parámetros que son desconocidos. Estos se suelen 12 modelos de distribuciones obtenidas de la distribución normal estimar por MV, obteniendo así la matriz de información de Fisher observada; ( 𝜃). La matriz inversa de ( 𝜃), 𝑉=[( 𝜃)]−1, nos da estimaciones de la varianza y covarianza de los estimadores de los parámetros. En concreto, los elementos de la diagonal nos dan las varianzas, y los de fuera de la diagonal, la covarianza. Este resultado se utilizará para aproximar el error estándar de los EMV de los 𝜃𝑖,𝑖= 1,2, y para obtener intervalos de confianza aproximados. 1.3 Inferencia en muestras censuradas En el estudio de datos de supervivencia, es muy común encontrar que las muestras de datos puedan ser censuradas, esto es, nos dan información incompleta por algún motivo, el cual dependerá del tipo de censura. De esta manera, se necesita un mecanismo diferente para tratar con estos datos a los usados con las muestras completas. Hay diferentes tipos de censura, como podemos ver a continuación. Censura de tipo I: Se da cuando el investigador fija un tiempo máximo para la observación del tiempo de vida de los datos. Entonces, cuando un individuo observado no presenta una falla dentro de ese tiempo, no aporta la información necesaria al investigador y constituye una observación censurada. Censura de tipo II: Supongamos que se observa una muestra de 𝑛elementos y se fija un número 𝑟tal que 1≤𝑟≤𝑛. En este caso, el experimento termina cuando 𝑟de los individuos fallan, resultando por tanto 𝑟datos observados tales que 𝑡1∶𝑛≤⋯≤𝑡𝑟∶𝑛,y𝑛−𝑟datos censurados, puesto que no han fallado al acabar la investigación. Por 𝑡𝑖∶𝑛se ha denotado al estadístico de orden 𝑖en una muestra de tamaño 𝑛. Censura aleatoria: En este tipo, la censura no depende de ninguna decisión del investigador. Se da cuando cada individuo bajo estudio tiene un tiempo de vida 𝑇𝑖y un tiempo de censura 𝐶𝑖, siendo 𝑇𝑖y𝐶𝑖variables aleatorias continuas e independientes entre ellas. Para modelizar este tipo de datos, se utilizará la notación que sigue. Supongamos una muestra aleatoria de 𝑛individuos, cuyos tiempos de vida vienen dados por las variables aleatorias 𝑇1,…, 𝑇𝑛, respectivamente. Llamamos 𝑡𝑖a lo que puede ser o bien el tiempo de vida observado o bien el tiempo de censura según lo que ocurra primero, esto es, 𝑡𝑖= m ın{𝑇𝑖, 𝐶𝑖}, siendo 𝐶𝑖el tiempo de censura en la observación 𝑖. Definimos 1. herramientas básicas 13 la variable 𝛿𝑖,𝑖= 1,…, 𝑛, tal que 𝛿𝑖={1 si 𝑡𝑖dato observado 0 si 𝑡𝑖dato censurado De esta manera, las muestras censuradas se escriben como (𝑡𝑖, 𝛿𝑖),𝑖= 1,…, 𝑛. Método de máxima verosimilitud Para poder hacer inferencia sobre este tipo de observaciones, se suele utilizar el método de máxima verosimilitud, que se verá alterado por el carácter de estos datos. De forma general, denotaremos al vector con los parámetros del modelo bajo estudio como 𝜃= (𝜃1, 𝜃2). La función de verosimilitud cuando partimos de una muestra de datos censurados (𝑡𝑖, 𝛿1),𝑖= 1,…, 𝑛, viene dada por la siguiente expresión 𝐿(𝜃) = 𝑛 ∏ 𝑖=1 𝑓(𝑡𝑖)𝛿𝑖𝑆(𝑡𝑖)1−𝛿𝑖,(1.6) donde 𝑓es la fdd de la muestra y 𝑆denota la función de supervivencia correspondiente. Esta forma de la función de verosimilitud es común a los tres tipos de datos censurados nombrados anteriormente (ver Lawless [11]). También puede ser de utilidad la función de logverosimilitud, cuya expresión para este tipo de muestras es 𝑙(𝜃) = log 𝐿(𝜃) = 𝑛 ∑ 𝑖=1 𝛿𝑖log 𝑓(𝑡𝑖) + 𝑛 ∑ 𝑖=1 (1 − 𝛿𝑖) log 𝑆(𝑡𝑖).(1.7) Matriz de información de Fisher En el caso de muestras censuradas, sigue teniendo sentido considerar la matriz de información de Fisher. Al igual que en el caso de muestras completas, viene definida como la matriz donde cada elemento (𝑖, 𝑗)viene dado por −𝜕𝑙(𝜃) 𝜕𝜃𝑖𝜕𝜃𝑗 , 𝑖, 𝑗 = 1,2. En este caso, 𝑙(𝜃)representa la función de log-verosimilitud de la muestra censurada, dada por (1.7). 14 modelos de distribuciones obtenidas de la distribución normal Esta matriz se utilizará para aproximar el error estándar y covarianza asintótica de los estimadores de máxima verosimilitud. 1.4 Distribución Normal En esta sección se van a recoger algunos resultados importantes de la distribución normal, que nos serán útiles para abordar el estudio de modelos más complejos en capítulos posteriores. Definición 1.6. Sean 𝜇∈ℝy𝜎 > 0, y sea 𝑋una variable aleatoria (v.a.). Decimos que 𝑋sigue una distribución normal univariante de parámetros 𝜇y𝜎, y se denota por 𝑋∼𝑁(𝜇, 𝜎2), si su función de densidad (fdd) se expresa como sigue 𝑓(𝑥;𝜇, 𝜎) = 1 √2𝜋𝜎 exp {−(𝑥−𝜇)2 2𝜎2},−∞ <𝑥<∞.(1.8) Es conveniente remarcar el caso particular de la distribución normal donde 𝜇= 0 y𝜎= 1, debido a su gran utilidad en múltiples campos de la estadística y de este mismo trabajo. Definición 1.7. Sea una v.a. 𝑍tal que 𝑍∼𝑁(0,1). Entonces, se dice que 𝑍sigue una distribución normal estándar univariante La fdd de 𝑍se denota por 𝜙, y viene dada por 𝜙(𝑥) = 1 √2𝜋 𝑒−𝑥2 2,−∞ <𝑥<∞,(1.9) mientras que su FdD, denotada como Φ, es de la siguiente forma Φ(𝑥) = ∫𝑥 −∞ 𝜙(𝑢)𝑑𝑢. (1.10) La siguiente proposición nos da algunas de las propiedades de la distribución normal estándar, que serán importante para el posterior estudio de las propiedades tanto de la normal como de otras distribuciones. 1. herramientas básicas 15 Proposición 1.3. Si 𝑍∼𝑁(0,1), entonces se verifica 1. 𝜙es una función par, i.e., 𝜙(𝑥) = 𝜙(−𝑥),∀𝑥∈ℝ. 2. Φ(−𝑥) = 1 − Φ(𝑥) 3. Función generatriz de momentos de Z: 𝑀𝑍(𝑡) = 𝑒1 2𝑡2, 𝑡 ∈ℝ.(1.11) 4. 𝐸(𝑍)=0;𝑉 𝑎𝑟(𝑍)=1. Demostración. 1. Trivial a partir de (1.9). 2. Inmediato por la simetría de 𝜙. 3. La FGM se define como 𝑀𝑍(𝑡) = 𝐸(𝑒𝑡𝑍).Entonces 𝑀𝑍(𝑡) = 1 √2𝜋∫∞ −∞ 𝑒𝑡𝑧𝑒−𝑧2 2𝑑𝑧 =𝑒𝑡2 2 √2𝜋∫∞ −∞ 𝑒𝑡𝑧−𝑧2 2−𝑡2 2𝑑𝑧 =𝑒𝑡2 2 √2𝜋∫∞ −∞ 𝑒−1 2(𝑧−𝑡)2𝑑𝑧. Tomamos el cambio de variable 𝑟=𝑧−𝑡, resultando así 𝑀𝑍(𝑡) = 𝑒𝑡2 2∫∞ −∞ 1 √2𝜋 𝑒−1 2(𝑟)2𝑑𝑟 =𝑒𝑡2 2∫∞ −∞ 𝜙(𝑟)𝑑𝑟 =𝑒𝑡2 2, 𝑡 ∈ℝ. 4. Recordemos que 𝐸(𝑍𝑘) = 𝑀𝑘)(𝑡)|𝑡=0. Entonces 𝑑𝑀𝑍 𝑑𝑡 (𝑡) = 𝑡𝑒𝑡2 2⇒𝐸(𝑍)=0. Sabemos que 𝑉 𝑎𝑟(𝑍) = 𝐸(𝑍2) − 𝐸(𝑍)2. Acabamos de probar que 𝐸(𝑍)=0, luego basta con hallar el valor de 𝐸(𝑍2). 𝐸(𝑍2) = 𝑑2𝑀𝑍 𝑑𝑡2(𝑡) = 𝑒𝑡2 2+𝑡2𝑒𝑡2 2⇒𝑉 𝑎𝑟(𝑍)=1. Definición 1.8. Sea 𝑍∼𝑁(0,1). Para 0<𝑝<1, se denotará por 𝑧𝑝al cuantil de orden p, es decir, al valor 𝑧𝑝tal que Φ(𝑧𝑝) = 𝑝. Sea 𝑋∼𝑁(𝜇, 𝜎2). Denotamos por 𝑥𝑝al cuantil de orden p. Puesto que podemos expresar 𝑋=𝜇+𝜎𝑍 con 𝑍∼𝑁(0,1), se tiene la siguiente relación entre los cuantiles de ambas distribuciones 𝑥𝑝=𝜇+𝜎𝑧𝑝. 16 modelos de distribuciones obtenidas de la distribución normal La siguiente proposición es esencial en el estudio de la distribución normal, ya que nos da una relación entre esta y la distribución normal estándar. A esta relación se le llama estandarización. Proposición 1.4 (Estandarización). Sea 𝑋∼𝑁(𝜇, 𝜎2)y sea 𝑍=𝑋−𝜇 𝜎. Entonces, 𝑍∼𝑁(0,1). Demostración. Sea 𝑍=𝑋−𝜇 𝜎. Podemos calcular su FdD como sigue 𝐹𝑍(𝑧) = 𝑃(𝑍≤𝑧) = 𝑃(𝑋−𝜇 𝜎≤𝑧)=𝑃(𝑋≤𝜎𝑧 +𝜇)=𝐹𝑋(𝜎𝑧 +𝜇). A partir de esto, podemos obtener su fdd derivando 𝐹𝑍y usando (1.8). 𝑓𝑍(𝑧) = 𝑑𝐹𝑋 𝑑𝑧 (𝜎𝑧 +𝜇) = 𝜎𝑓𝑋(𝜎𝑧 +𝜇) = 𝜎1 √2𝜋𝜎 exp {−(𝜎𝑧 +𝜇−𝜇)2 2𝜎2}. Se tiene entonces 𝑓𝑍(𝑧) = 𝜙(𝑧) ∀𝑥∈ℝ, luego 𝑍∼𝑁(0,1) Estas propiedades nos dan las herramientas necesarias para ver distintas características de la distribución normal, empezando por el siguiente corolario, que nos da una caracterización de la FdD de la distribución. Corolario 1.1.Dada 𝑋∼𝑁(𝜇, 𝜎2), esta puede expresarse como 𝑋=𝜇+𝜎𝑍, y su FdD vendrá dada por 𝐹(𝑥)=Φ(𝑥−𝜇 𝜎).(1.12) Proposición 1.5 (Propiedades). Si 𝑋∼𝑁(𝜇, 𝜎2), entonces se verifica 1. 𝑓(𝑥)es simétrica respecto de 𝜇:𝑓(𝑥) = 𝑓(2𝜇−𝑥),∀𝑥∈ℝ 2. 𝐸(𝑋) = 𝜇,𝑉 𝑎𝑟(𝑋) = 𝜎2 3. Función generatriz de momentos de X: 𝑀𝑋(𝑡) = exp {𝑡𝜇 +1 2𝜎2𝑡2}, 𝑡 ∈ℝ.(1.13) Demostración. 1. Tomando 𝑥= 2𝜇−𝑥en (1.8), se observa 𝑓(2𝜇−𝑥) = 1 √2𝜋𝜎 exp {−(2𝜇−𝑥−𝜇)2 2𝜎2}=𝜙(𝑥), 𝑥 ∈ℝ. 1. herramientas básicas 17 2. Sabemos que se cumple la relación 𝑋=𝜇+𝜎𝑍, donde 𝑍∼𝑁(0,1). Entonces, mediante propiedades conocidas de la esperanza y la varianza se tiene que 𝐸(𝑋) = 𝐸(𝜇+𝜎𝑍) = 𝐸(𝑋) = 𝜇+𝜎𝐸(𝑍) = 𝜇. 𝑉 𝑎𝑟(𝑋) = 𝑉 𝑎𝑟(𝜇+𝜎𝑍) = 𝜎2𝑉 𝑎𝑟(𝑍) = 𝜎2. 3. Aplicando de nuevo la transformación lineal 𝑋=𝜇+𝜎𝑍, 𝑀𝑋(𝑡) = 𝐸(𝑒𝑡𝑋) = 𝐸(𝑒𝑡(𝜇+𝜎𝑍)=𝑒𝑡𝜇𝐸(𝑒𝑡𝜎𝑍). Tomamos el cambio de variable 𝑡𝜎 =𝑢, resultando 𝑀𝑋(𝑡) = 𝑒𝑡𝜇𝐸(𝑒𝑢𝑍)=𝑒𝑡𝜇𝑒−1 2𝑢2=𝑒𝑡𝜇𝑒−1 2𝑢2. Finalmente, deshaciendo el cambio, 𝑀𝑋(𝑡) = exp {𝑡𝜇 +1 2𝜎2𝑡2}, 𝑡 ∈ℝ. La función de densidad de 𝑋∼𝑁(𝜇, 𝜎2)tiene forma de campana (campana de Gauss). El parámetro 𝜇es la media, mediana y moda de la distribución. Por otro lado, 𝜎es la desviación típica, y también un parámetro de escala. En la figura (1.1) pueden observarse estos hechos junto con el efecto de variar los valores de estos parámetros. Proposición 1.6. Sea 𝑋∼𝑁(𝜇, 𝜎2). Entonces, 𝑎𝑋 +𝑏∼𝑁(𝑎𝜇 +𝑏, 𝑎2𝜎2)para 𝑎, 𝑏 ∈ℝ. Demostración. Sea 𝑋∼𝑁(𝜇, 𝜎2). Veamos que 𝑌=𝑎𝑋 +𝑏sigue una distribución normal, usando la FGM como caracterización de la misma. 𝑀𝑌(𝑡) = 𝐸(𝑒𝑡𝑌 ) = 𝑒𝑡𝑏𝑀𝑋(𝑎𝑡) = exp {𝑡(𝑎𝜇 +𝑏) + 1 2(𝑎𝜎)2𝑡2}, luego 𝑋∼𝑁(𝑎𝜇 +𝑏, 𝑎2𝜎2). Proposición 1.7 (Reproductividad). Sean 𝑋1∼𝑁(𝜇1, 𝜎2 1)y𝑋2∼𝑁(𝜇2, 𝜎2 2) independientes entre sí. Entonces 𝑋1+𝑋2∼𝑁(𝜇1+𝜇2, 𝜎2 1+𝜎2 2). 18 modelos de distribuciones obtenidas de la distribución normal −6 −2 0 2 4 6 0.0 0.1 0.2 0.3 0.4 Comparativa N(0,1) y N(0,9) x σ = 1 σ = 3 −6 −2 0 2 4 6 0.0 0.1 0.2 0.3 0.4 Comparativa N(1,1) y N(−1,1) x µ = −1µ = 1 Figura 1.1: fdd de la distribución Normal con distintos parámetros Demostración. Tomamos las expresiones de las FGM de 𝑋1y𝑋2vistas en (1.13). Debido a la independencia entre ambas variables, se cumple lo siguiente 𝑀𝑋1+𝑋2(𝑡) = 𝑀𝑋1(𝑡)𝑀𝑋2(𝑡) = exp {𝑡(𝜇1+𝜇2) + 1 2𝑡2(𝜎2 1+𝜎2 2)}, que se corresponde con la expresión de la FGM de una distribución normal de media 𝜇1+𝜇2y varianza 𝜎2 1+𝜎2 2. Por tanto, 𝑋1+𝑋2∼𝑁(𝜇1+𝜇2, 𝜎2 1+𝜎2 2). Esta propiedad se generaliza de forma inmediata a la suma de 𝑛v.a.s independientes. Destacamos el caso en que además de independientes, las v.a.s son idénticamente distribuidas. Corolario 1.2.Sean 𝑋1, 𝑋2, ..., 𝑋𝑛v.a.s independientes e idénticamente distribuidas (i.i.d.) según 𝑁(𝜇, 𝜎2).Entonces 1. ∑𝑛 𝑖=1 𝑋𝑖∼𝑁(𝑛𝜇, 𝑛𝜎2). 2.  𝑋=1 𝑛∑𝑛 𝑖=1 𝑋𝑖∼𝑁(𝜇, 𝜎2 𝑛),donde  𝑋es la media muestral. La propiedad anterior es de gran utilidad al hacer inferencia en poblaciones normales, junto con los dos resultados que siguen. En el primero, destacamos que el mo- 1. herramientas básicas 19 delo normal pertenece a la familia exponencial de distribuciones, y en el Teorema de Fisher, se establece la independencia de la media muestral,  𝑋, y la cuasivarianza muestral , 𝑆2 𝑐=∑𝑛 𝑖=1(𝑋𝑖− 𝑋)2∕(𝑛− 1), junto con la distribución de estos dos estadísticos en el caso normal. Proposición 1.8. Sea 𝑋∼𝑁(𝜇, 𝜎2). Nótese que podemos escribir su fdd como 𝑓(𝑥) = 1 𝜎√2𝜋 exp {−𝑥2 2𝜎2−𝜇2 2𝜎2+2𝜇𝑥 2𝜎2}=1 𝜎√2𝜋 exp {𝜇𝑥 𝜎2−𝑥2 2𝜎2−𝜇2 2𝜎2}, donde 𝜇∈ℝ, 𝜎 > 0. Pertenece a la familia exponencial biparamétrica de distribuciones. Dada una muestra aleatoria simple (variables aleatorias i.i.d.) 𝑋1, 𝑋2, ..., 𝑋𝑛∼𝑁(𝜇, 𝜎2), entonces (∑𝑛 𝑖=1 𝑋𝑖,∑𝑛 𝑖=1 𝑋2 𝑖)es un estadístico suficiente y completo para (𝜇, 𝜎2). Equivalentemente, ( 𝑋, 𝑆2 𝑐)también es un estadístico suficiente y completo para (𝜇, 𝜎2). Teorema 1.1 (Fisher). Sean 𝑋1, ..., 𝑋𝑛variables aleatorias i.i.d. según 𝑁(𝜇, 𝜎2). Entonces 1.  𝑋y𝑆2 𝑐son independientes, siendo  𝑋es la media muestral y 𝑆2 𝑐la cuasivarianza muestral. 2.  𝑋=1 𝑛∑𝑛 𝑖=1 𝑋𝑖∼𝑁(𝜇, 𝜎2 𝑛). 3. (𝑛−1) 𝜎2𝑆2 𝑐=1 𝜎2∑𝑛 𝑖=1(𝑋𝑖− 𝑋)2∼𝜒2 𝑛−1. Este resultado es la base para realizar inferencia en la distribución normal cuando tenemos muestras completas. 26 modelos de distribuciones obtenidas de la distribución normal 0 2 4 6 8 10 12 0.0 0.5 1.0 1.5 Funciones de densidad Log−Normal x µ = 0 µ = 1 µ = 2 Figura 2.3: Función de densidad Log Normal con 𝜎= 0.25 y𝜇=1,2 y 3. máxima verosimilitud de los parámetros 𝜇y𝜎de la distribución. Muestras completas Para muestras aleatorias completas, es fácil obtener 𝜇 y𝜎, ya que tomando la transformación 𝑦𝑖= log 𝑡𝑖, reducimos la muestra a una que sigue una distribución normal. Usando (2.1), calculamos la función de verosimilitud y la función de log-verosimilitud. 𝐿(𝜇, 𝜎) = 𝑛 ∏ 𝑖=1 𝑓(𝑡𝑖) = 𝑛 ∏ 𝑖=1 1 √2𝜋𝜎𝑡𝑖 exp {−1 2(log 𝑡𝑖−𝜇 𝜎)2},(2.8) 𝑙(𝜇, 𝜎) = log 𝐿(𝜇, 𝜎)=−𝑛 2log 2𝜋−𝑛log 𝜎− 𝑛 ∑ 𝑖=1 (log 𝑡𝑖+1 2(log 𝑡𝑖−𝜇 𝜎)2).(2.9) Calculamos las derivadas de (2.9) e igualamos a cero para obtener los estimadores 𝜇 y𝜎. 𝜕𝑙 𝜕𝜇(𝜇, 𝜎) = 𝑛 ∑ 𝑖=1 log 𝑡𝑖−𝜇 𝜎2= 0 ⇒ 𝑛 ∑ 𝑖=1 log 𝑡𝑖=𝑛𝜇 ⇒𝜇 =1 𝑛 𝑛 ∑ 𝑖=1 log 𝑡𝑖. También puede expresarse como 𝜇 =𝑦 =1 𝑛 𝑛 ∑ 𝑖=1 𝑦𝑖,(2.10) 2. distribución log normal 27 donde 𝑦 es la media muestral de (𝑦1, ..., 𝑦𝑛) = (log 𝑡1, ..., log 𝑡𝑛). Por otra parte, 𝜕𝑙 𝜕𝜎(𝜇, 𝜎)=−𝑛 𝜎+ 𝑛 ∑ 𝑖=1 (log 𝑡𝑖−𝜇)2 𝜎3= 0 ⇒𝜎 =(𝑛 ∑ 𝑖=1 (log 𝑡𝑖−𝜇)2 𝑛)1∕2. Calculando la matriz hessiana y evaluándola en (𝜇, 𝜎), se tiene que (𝜇, 𝜎)corresponde a un máximo. Teniendo en cuenta la expresión alternativa de 𝜇, la desviación típica de las 𝑦𝑖se puede escribir también como 𝜎 =(𝑛 ∑ 𝑖=1 (𝑦𝑖−𝑦)2 𝑛)1∕2 (2.11) o bien, como 𝜎 =[(𝑛− 1)𝑠2 𝑡∕𝑛]1∕2 donde 𝑠2 𝑡=∑(𝑦𝑖−𝑦)2∕(𝑛− 1) es la cuasivarianza muestral de las 𝑦𝑖. Aplicando el Teorema de Fisher, se tiene que 𝜇 y 𝜎2son independientes. Además 𝜇 ∼𝑁(𝜇, 𝜎2 𝑛),(𝑛− 1)  𝜎2 𝜎2∼𝜒2 𝑛−1. Por tanto, pueden aplicarse resultados para poblaciones normales al realizar inferencia sobre estos parámetros. Muestras censuradas Sea una muestra aleatoria censurada de tiempos de vida (𝑡𝑖, 𝛿𝑖), 𝑖 = 1, ..., 𝑛, tal que sigue una distribución Log Normal. Para trabajar con estos datos, utilizaremos la fdd de la distribución, 𝑓(𝑡), ya vista en (2.1), y la función de supervivencia 𝑆(𝑡), cuya expresión se encuentra en (2.3). Si denotamos 𝑦𝑖= log 𝑡𝑖, 𝑧𝑖=𝑦𝑖−𝜇 𝜎, 𝑆0(𝑧) = 1 − Φ(𝑧), 𝑓0(𝑧) = 𝜙(𝑧) = 1 √2𝜋 𝑒−𝑧2∕2, entonces se tiene que 𝑓(𝑡𝑖) = 1 𝑡𝑖 1 𝜎𝜙(log 𝑡𝑖−𝜇 𝜎)=1 𝑡𝑖 1 𝜎𝑓0(𝑧), 𝑆(𝑡𝑖) = 1 − Φ(log 𝑡𝑖−𝜇 𝜎)=𝑆0(𝑧𝑖). 28 modelos de distribuciones obtenidas de la distribución normal Usando (1.6), la función de verosimilitud vendrá entonces dada por 𝐿(𝜇, 𝜎) = 𝑛 ∏ 𝑖=1 [1 𝑡𝑖]𝑑𝑖𝑛 ∏ 𝑖=1 [1 𝜎𝑓0(𝑦𝑖−𝜇 𝜎)]𝛿𝑖 𝑆0(𝑦𝑖−𝜇 𝜎)1−𝛿𝑖 .(2.12) Como ∏𝑛 𝑖=1 [1 𝑡𝑖]𝑑𝑖no depende de ninguno de los parámetros, para calcular la función de log-verosimilitud nos basta con tener en cuenta la siguiente relación 𝐿(𝜇, 𝜎) ∝ 𝑛 ∏ 𝑖=1 [1 𝜎𝑓0(𝑦𝑖−𝜇 𝜎)]𝛿𝑖 𝑆0(𝑦𝑖−𝜇 𝜎)1−𝛿𝑖 . Si llamamos 𝑟al número de datos realmente observados, esto es, 𝑟=∑𝑛 𝑖=1 𝛿𝑖, entonces la función de log-verosimilitud para 𝜇y𝜎viene dada por 𝑙(𝜇, 𝜎)=−𝑟log 𝜎+ 𝑛 ∑ 𝑖=1 [𝛿𝑖log 𝑓0(𝑧𝑖) + (1 − 𝛿𝑖) log 𝑆0(𝑧𝑖)],(2.13) donde al sustituir 𝑓0, es lo mismo que 𝑙(𝜇, 𝜎)=−𝑟log 𝜎−1 2 𝑛 ∑ 𝑖=1 𝛿𝑖𝑧2 𝑖+ 𝑛 ∑ 𝑖=1 (1 − 𝛿𝑖) log 𝑆0(𝑧𝑖).(2.14) Es claro que 𝜕𝑧𝑖 𝜕𝜇 = − 1 𝜎,𝜕𝑧𝑖 𝜕𝜎 =−(𝑦𝑖−𝜇) 𝜎2=−𝑧𝑖 𝜎. Usamos esto para calcular sus derivadas a partir de (2.14) y obtener así las ecuaciones de verosimilitud, formadas por las siguientes expresiones igualadas a cero. 𝜕𝑙 𝜕𝜇 =1 𝜎 𝑛 ∑ 𝑖=1 (𝛿𝑖𝑧𝑖+ (1 − 𝛿𝑖)𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖)),(2.15) 𝜕𝑙 𝜕𝜎 = − 𝑟 𝜎+1 𝜎 𝑛 ∑ 𝑖=1 (𝛿𝑖𝑧2 𝑖+ (1 − 𝛿𝑖)𝑧𝑖 𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖)).(2.16) Los estimadores 𝜇 y𝜎 no tienen una forma explícita para los datos censurados, siendo necesario un método numérico para resolver las ecuaciones de verosimilitud y obtener así los estimadores que maximizan (2.14). 2. distribución log normal 29 A continuación, calculamos las derivadas segundas, ya que son necesarias para calcular la matriz de información observada. Estas vienen dadas por 𝜕2𝑙 𝜕𝜇2=1 𝜎2 𝑛 ∑ 𝑖=1 (−𝛿𝑖+ (1 − 𝛿𝑖)𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖)(𝑧𝑖−𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖))),(2.17) 𝜕2𝑙 𝜕𝜎2=𝑟 𝜎2−2 𝜎2 𝑛 ∑ 𝑖=1 (𝛿𝑖𝑧2 𝑖+ (1 − 𝛿𝑖)𝑧𝑖 𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖)) +1 𝜎2 𝑛 ∑ 𝑖=1 (−𝛿𝑖𝑧𝑖+ (1 − 𝛿𝑖)𝑧2 𝑖 𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖)(𝑧𝑖−𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖))), (2.18) 𝜕2𝑙 𝜕𝜇𝜕𝜎 = − 1 𝜎2 𝑛 ∑ 𝑖=1 (𝛿𝑖𝑧2 𝑖+ (1 − 𝛿𝑖)𝑧𝑖 𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖)) +1 𝜎2 𝑛 ∑ 𝑖=1 (−𝛿𝑖𝑧2 𝑖+ (1 − 𝛿𝑖)𝑧2 𝑖 𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖)(𝑧𝑖−𝑓0(𝑧𝑖) 𝑆0(𝑧𝑖))). (2.19) Entonces, para obtener la matriz de información observada de Fisher basta con sustituir en la siguiente expresión, donde todos los valores son conocidos. 𝐼(𝜇, 𝜎) = (−𝜕2𝑙∕𝜕𝜇2−𝜕2𝑙∕𝜕𝜇𝜕𝜎 −𝜕2𝑙∕𝜕𝜎𝜕𝜇 −𝜕2𝑙∕𝜕𝜎2).(2.20) La matriz 𝐼(𝜇, 𝜎)nos da la matriz de covarianzas asintótica  𝑉=[𝐼(𝜇, 𝜎)]−1 de los estimadores (𝜇, 𝜎). 2.3 Aplicaciones Aplicación 1. Muestra completa Los datos utilizados a continuación han sido obtenidos del trabajo de Overduin [15]. Esta serie de datos es completa, y representa los tiempos de vida en meses de 184 pacientes fallecidos por cáncer de pulmón. Estos datos son los descritos en el Cuadro 2.1, y su representación gráfica puede observarse en el histograma de la Figura 2.4. Usaremos las expresiones calculadas a lo largo de este capítulo para estudiar e interpretar los datos dados. 30 modelos de distribuciones obtenidas de la distribución normal Datos de cáncer de pulmón 4,04 4,70 5,82 6,15 7,07 7,36 7,56 7,76 7,82 7,86 7,86 7,89 8,15 8,19 8,84 9,04 9,17 9,24 9,47 9,67 10,03 10,06 10,13 10,26 10,32 10,36 10,42 10,42 10,45 10,52 10,52 10,72 10,75 10,75 11,15 11,18 11,28 11,34 11,47 11,77 11,80 11,93 12,03 12,30 12,39 12,53 12,53 12,53 12,56 12,56 12,82 12,95 13,05 13,12 13,15 13,18 13,28 13,32 13,74 13,91 14,04 14,17 14,33 14,33 14,93 14,93 14,99 15,02 15,01 15,12 15,35 15,52 15,58 15,88 15,95 15,95 16,01 16,11 16,14 16,27 16,41 16,41 16,60 16,67 16,77 17,13 17,16 17,23 17,52 17,79 17,82 17,98 18,02 18,02 18,48 18,61 18,81 18,81 19,13 19,17 19,20 19,20 19,30 19,46 19,53 19,63 19,73 19,82 19,86 19,89 20,05 20,12 20,19 10,22 20,28 20,32 20,65 20,65 20,68 20,68 20,78 20,81 20,84 21,11 21,14 21,47 21,50 21,70 21,80 21,90 22,45 22,62 23,31 23,54 23,57 23,64 23,70 23,70 23,70 23,84 24,03 24,16 24,20 24,46 24,46 24,69 24,72 24,79 25,18 25,35 25,45 25,97 25,97 27,12 27,16 27,48 27,65 28,04 28,27 28,64 29,10 29,98 30,02 30,05 30,97 31,27 32,55 32,61 33,83 34,88 35,38 36,62 38,37 42,38 43,00 44,42 44,65 47,28 47,64 53,82 55,69 58,50 58,82 64,64 Cuadro 2.1: Tiempos de vida en meses de pacientes de cáncer de pulmón. Tenemos entonces una muestra aleatoria (𝑡1,…, 𝑡𝑛)de 𝑛= 184 observaciones. El histograma de la muestra nos da una idea de que los datos siguen una distribución log-normal, ya que se puede ver una asimetría muy fuerte a la derecha, lo cual es muy característico de esta distribución. Podemos comprobar esto mediante un test de hipótesis, ya que si llamamos 𝑌= log 𝑇y tomamos 𝑦𝑖= log 𝑡𝑖,𝑖= 1,…, 𝑛, entonces la muestra (𝑦1,…, 𝑦𝑛)debe seguir una distribución normal. Se puede aplicar a la variable Y el test de normalidad de Shapiro-Wilk, que contrasta la hipótesis 𝐻0∶𝑌∼𝑁𝑜𝑟𝑚𝑎𝑙 frente a la hipótesis alternativa 𝐻1∶𝑌≁𝑁𝑜𝑟𝑚𝑎𝑙. Utilizando el programa R (ver [16]) para implementar este contraste, se obtiene un p-valor igual a 0,5595, por lo que podemos suponer la 2. distribución log normal 31 Cáncer de pulmón (n=184) Tiempo de vida en meses Frecuencia 0 10 20 30 40 50 60 0 10 20 30 40 Figura 2.4: Histograma obtenido de los datos del Cuadro 2.1. hipótesis 𝐻0. Esto implica de forma inmediata que, en efecto, la variable 𝑇sigue una distribución log-normal. Para calcular los estimadores de los parámetros, partimos de los resultados de máxima verosimilitud dados por (2.10) y (2.11). 𝜇 = 184 ∑ 𝑖=1 1 184 log 𝑡𝑖= 2,855664. 𝜎 =(184 ∑ 𝑖=1 (log 𝑡𝑖−𝜇)2 184 )1∕2 = 0,4965. Podemos suponer de esta manera que la muestra aleatoria dada sigue una distribución 𝑇∼𝐿𝑜𝑔𝑁(𝜇, 𝜎2).En la Figura 2.5 se observa que el modelo propuesto con estos estimadores se ajusta correctamente a los datos. Podemos estimar la media utilizando su expresión en (2.5), resultando  𝐸(𝑇) = 19,66692. La información que nos da este resultado es que al diagnosticar a un paciente de cáncer de pulmón, se puede esperar que su tiempo de vida desde entonces sea de 19 meses aproximadamente. También resulta interesante el estudio de algunos de los cuantiles de los datos. Para estimarlos, usamos (2.7) con los estimadores 𝜇 y𝜎 ya calculados. 32 modelos de distribuciones obtenidas de la distribución normal Cáncer de pulmón (n=184) Tiempo de vida en meses Frecuencias relativas 0 10 20 30 40 50 60 0.00 0.02 0.04 0.06 Histograma relativo Densidad Log N Figura 2.5: Histograma relativo de los datos y curva de la fdd 𝐿𝑜𝑔𝑁(2′8556,0′4965). Cuantil de orden 0,05  𝑡0,05 = 7,682425 Cuantil de orden 0,1  𝑡0,1= 9,201144 Mediana  𝑡0,5= 17,38598 Cuantil de orden 0,9  𝑡0,9= 32,85161 Cuantil de orden 0,95  𝑡0,95 = 39,34596 Las estimaciones obtenidas se corresponden con las conclusiones que podríamos sacar de la Figura 2.5. Podemos ver que la mitad de los individuos bajo estudio llegan a un tiempo de vida aproximado de 17 meses. Por otro lado, la información que nos aportan los demás cuantiles estimados es que se puede esperar que los pacientes vivan más de 7 y 9 meses tras el diagnóstico, ya que los casos en que el tiempo de vida es menor corresponde a apenas el 5 y el 10 % de los pacientes. Por otra parte, observamos que de un estudio de unos 60 meses, el 90 % de los individuos estudiados fallecen antes de los 32 meses, y el 95 %; antes de los 39 meses. Siendo así, podemos concluir que son muy excepcionales los casos en que un paciente sobrepase este límite. 2. distribución log normal 33 Aplicación 2. Muestra censurada El siguiente ejemplo se ha obtenido del trabajo de Lawless [11]. Cada dato de este conjunto recoge el número de miles de millas recorridas por una locomotora antes de que se produzca un fallo en sus controles. El estudio se realiza durante 135.000 millas y sobre una muestra de 96 locomotoras diferentes, de las cuales 37 presentan un fallo antes del fin del estudio. Se tienen entonces 37 muestras observadas y 59 muestras censuradas iguales a 135. En el Cuadro 2.2 podemos observar los datos, donde el asterisco (*) señala un dato censurado. Su representación gráfica en forma de histograma se encuentra en la Figura 2.6. Datos locomotoras Miles de millas recorridas Frecuencia 20 40 60 80 100 120 140 0 20 40 60 Figura 2.6: Histograma de datos Cuadro 2.2 Utilizamos el paquete "survival" de R (ver [18]) para calcular los estimadores, resultando 𝜇 = 5,116925 𝜎 = 0,705494 En la Figura 2.7 se puede ver el ajuste resultante mediante estos estimadores. La densidad calculada se ajusta bien hasta los últimos valores, debido a que estos últimos datos se corresponden con los de censura. Podemos considerarlo, por tanto, un buen ajuste. 34 modelos de distribuciones obtenidas de la distribución normal Millas recorridas por las locomotoras 22,5 37,5 46,0 48,5 51,5 53,0 54,5 57,5 66,5 68,0 69,5 76,5 77,0 78,5 80,0 81,5 82,0 83,0 84,0 91,5 93,5 102,5 107,0 108,5 112,5 113,5 116,0 117,0 118,5 119,0 120,0 122,5 123,0 127,5 131,0 132,5 134,0 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* 135,0* Cuadro 2.2: Miles de millas recorridas antes de un fallo en la locomotora. Histograma y densidad Miles de millas recorridas Frecuencias relativas 20 40 60 80 100 120 140 0.000 0.015 0.030 Histograma relativo Densidad Log N Figura 2.7: Histograma relativo de datos y curva de densidad 𝐿𝑁(𝜇, 𝜎). Se puede obtener también la matriz de covarianzas asintóticas estimada de esta muestra utilizando R u otro software, resultando  𝑉=(0,01085 0,005729 0,005729 0,008686) Los valores de la diagonal de  𝑉se corresponden con la varianza estimada de 𝜇 y𝜎, respectivamente, mientras que el valor de la diagonal inversa nos da la estimación 2. distribución log normal 35 de la covarianza de dichos parámetros. Vemos que son valores muy próximos a cero, mostrando la poca variabilidad de estos valores y es coherente con la independencia de los estimadores, ya que su covarianza estimada es prácticamente cero. 42 modelos de distribuciones obtenidas de la distribución normal La función generatriz de cumulantes se corresponde con el logaritmo neperiano de la función característica (3.7). De esta manera, su expresión es 𝜓(𝑟) = log 𝜑𝑇(𝑟) = 𝜆 𝜇(1 − √1 − 2𝑖𝜇2𝑟 𝜆), 𝑟 ∈ℝ.(3.12) La expresión de los cumulantes de orden 𝑡será por tanto la que sigue 𝜅𝑡=1 𝑖𝑡 𝑑𝑡 𝑑𝑟𝑡(log 𝜑(𝑟))|||𝑟=0 (3.13) Los cumulantes nos dan una forma más sencilla de calcular propiedades de interés de la distribución, como la media, la varianza y los coeficientes de asimetría y de kurtosis. Proposición 3.3. Sea 𝑇∼𝐼𝐺(𝜇, 𝜆). Se verifica 1. 𝜅1=𝜇. 2. 𝜅2=𝜇3∕𝜆. 3. En general, para 𝑡≥2, se tiene 𝜅𝑡= 1 ⋅3⋅5⋅⋯⋅(2𝑡− 3)𝜇2𝑡−1𝜆1−𝑡. Demostración. 1. Usamos la función generatriz de cumulantes (3.12). 𝑑𝜓 𝑑𝑟 =𝜆 𝜇 2𝜇2∕𝜆 2√1 − 2𝜇2𝑟 𝜆 . Evaluando en 𝑟= 0, se tiene 𝜅1=𝜇. 2. De igual manera, 𝑑2𝜓 𝑑𝑟2= −𝜇 2(1 − 2𝜇2𝑟 𝜆)−3∕2(−2𝜇2 𝜆). Tomando 𝑟= 0, efectivamente 𝜅2=𝜇3∕𝜆. 3. Ver [19]. Observación 3.3.Se tiene que 𝜅2 𝜅3 1 =1 𝜆. Por ello, se dice que 𝜆es una medida inversa de precisión relativa, lo cual nos puede ser útil para posibles interpretaciones de 𝜆. 3. distribución inversa gaussiana 43 Proposición 3.4. Sea 𝜙=𝜆 𝜇. Entonces, los coeficientes de asimetría y de kurtosis de la distribución vienen dados por 1. Coeficiente de asimetría 𝛾1= 3𝜙−1∕2.(3.14) 2. Coeficiente de kurtosis 𝛾2= 15𝜙−1.(3.15) Demostración. Se tiene 𝜅3= 3𝜇5∕𝜆2y𝜅4= 15𝜇2∕𝜆3. Usando esto: 𝛾1=𝜅3∕𝜅3∕2 2= 3(𝜇∕𝜆)1∕2 𝛾2=𝜅4∕𝜅−2 2= 15(𝜇∕𝜆) Cuando 𝜙aumenta, la distribución se asemeja a la normal. Es por tanto un parámetro de forma o parámetro de normalidad. El siguiente resultado especifica esta relación. Proposición 3.5. Si 𝑇∼𝐼𝐺(𝜇, 𝜆), se tiene que cuando 𝜙→∞, entonces 𝑇 ←←←←←←→ 𝑁(𝜇, 𝜇3 𝜆). Un cambio de escala en la distribución conserva el modelo, salvo porque los parámetros se ven también multiplicados por el mismo factor que la variable. Lo vemos en la siguiente proposición. Proposición 3.6. Sea 𝑇∼𝐼𝐺(𝜇, 𝜆). Entonces, 𝑎𝑇 ∼𝐼𝐺(𝑎𝜇, 𝑎𝜆). Demostración. Llamamos 𝑌=𝑎𝑇 . Entonces, su FGM vendrá dada por 𝑀𝑌(𝑟) =𝐸[𝑒𝑟𝑌 ]=𝑀𝑇(𝑎𝑟) = exp {𝜆 𝜇(1 − √1 − 2𝜇2𝑎𝑟 𝜆)} = exp {𝑎𝜆 𝑎𝜇(1 − √1 − 2𝜇2𝑎2𝑟 𝑎𝜆 )}, de donde 𝑌∼𝐼𝐺(𝑎𝜇, 𝑎𝜎) 44 modelos de distribuciones obtenidas de la distribución normal Teorema 3.1. Si 𝑇∼𝐼𝐺(𝜇, 𝜆), su fdd 𝑓(𝑡)pertenece a la familia exponencial, y se puede expresar como 𝑓(𝑡, 𝜇, 𝜆) = exp {−𝜆(𝑡2− 2𝑡𝜇 +𝜇2) 2𝑡𝜇2+1 2log(𝜆) − 1 2log(2𝜋𝑡3)}= = exp {−𝜆 2𝜇2[𝑡− 2𝜇+𝜇2 𝑡]+1 2log(𝜆) − 1 2log(2𝜋𝑡3)}. 3.2 Inferencia de la distribución Inversa Gaussiana Sea una v.a. 𝑇que sigue una distribución 𝐼𝐺(𝜇, 𝜆), co 𝜇y𝜆desconocidos. En este apartado vamos a calcular 𝜇 y 𝜆, los estimadores de los parámetros, mediante el método de máxima verosimilitud, diferenciando entre los casos de muestras completas y censuradas. Muestras completas Sea la muestra aleatoria completa (𝑡1,…, 𝑡𝑛). Tomando la fdd de la distribución dada por (3.2), calculamos la función de verosimilitud 𝐿(𝜇, 𝜆) = 𝑛 ∏ 𝑖=1 𝑓𝑖(𝑡) = 𝑛 ∏ 𝑖=1 [(𝜆 2𝜋𝑡3 𝑖)1∕2]exp {𝑛 ∑ 𝑖=1 −𝜆(𝑡𝑖−𝜇)2 2𝜇2𝑡𝑖},(3.16) y a partir de ella, la función de log-verosimilitud 𝑙(𝜇, 𝜆) = 𝑛 2log (𝜆 2𝜋)−3 2 𝑛 ∑ 𝑖=1 log(𝑡𝑖) + 𝑛 ∑ 𝑖=1 −𝜆(𝑡𝑖−𝜇)2 2𝜇2𝑡𝑖 .(3.17) Calculamos las derivadas de (3.17) respecto de 𝜇y𝜆y las igualamos a cero para obtener 𝜇 y 𝜆respectivamente: 𝜕𝑙 𝜕𝜇(𝜇, 𝜆) = 𝑛 ∑ 𝑖=1 𝜆𝜇2(𝑡𝑖−𝜇) + 𝜆𝜇(𝑡𝑖−𝜇)2 𝜇4𝑡𝑖 = 0, de donde 𝜇 =1 𝑛 𝑛 ∑ 𝑖=1 𝑡𝑖=𝑡. (3.18) 3. distribución inversa gaussiana 45 Por otro lado, 𝜕𝑙 𝜕𝜆 =𝑛 2𝜆− 𝑛 ∑ 𝑖=1 (𝑡𝑖−𝜇)2 2𝜇2𝑡𝑖 = 0 ⇒ 𝜆=𝑛(𝑛 ∑ 𝑖=1 (𝑡𝑖−𝜇)2 𝜇2𝑡𝑖)−1, donde usando que 𝜇 =𝑡, 𝑛 ∑ 𝑖=1 (𝑡𝑖−𝜇)2 𝜇2𝑡𝑖 = 𝑛 ∑ 𝑖=1 (𝑡𝑖−𝑡)2 𝑡2𝑡𝑖 =1 𝑡2 𝑛 ∑ 𝑖=1 (𝑡𝑖− 2𝑡) + 𝑡 𝑡2 𝑛 ∑ 𝑖=1 𝑡2 𝑡𝑖 = − ∑1 𝑡+ 𝑛 ∑ 𝑖=1 1 𝑡𝑖 . Por tanto, se tiene  𝜆=𝑛 ∑𝑛 𝑖=1(𝑡−1 𝑖−𝑡−1) .(3.19) Los cálculos para obtener la matriz de información de Fisher son muy extensos. Se ha obtenido su expresión del trabajo de Henze [9], y es la siguiente 𝐼(𝜇, 𝜆) = (𝜆𝜇−3 0 0 (2𝜆2)−1)(3.20) Del trabajo de Tweedie [19], obtenemos algunos resultados interesantes relacionados con los estimadores de los parámetros, que enunciaremos a continuacion. Estos resultados pueden ser útiles para hacer inferencia sobre esta distribución. Dada una muestra aleatoria (𝑇1,…, 𝑇𝑛)de una distribución 𝐼𝐺(𝜇, 𝜆)con 𝜇 > 0, 𝜆 > 0, se tiene que el siguiente es un estadístico suficiente y completo para (𝜇, 𝜆). 𝑊(𝑇) = (𝑛 ∑ 𝑖=1 𝑇𝑖, 𝑛 ∑ 𝑖=1 1 𝑇𝑖). Teorema 3.2. Para los estimadores MV de los parámetros 𝜇y𝜆,𝜇 y 𝜆, se verifica 1. Los estimadores 𝜇 y 𝜆son independientes. 2.  𝜆es sesgado. 3.  𝜆verifica que 𝑛𝜆  𝜆∼𝜒2 𝑛−1, o lo que es lo mismo, 𝑛  𝜆∼1 𝜆𝜒2 𝑛−1. 4. El estimador 𝜇 cumple 𝜇 = 𝑡∼𝐼𝐺(𝜇, 𝑛𝜆). 46 modelos de distribuciones obtenidas de la distribución normal Muestras censuradas Consideramos la muestra aleatoria censurada (𝑡𝑖, 𝛿𝑖),𝑖= 1,…, 𝑛, que sigue una distribución 𝐼𝐺(𝜇, 𝜆).Para calcular los estimadores de máxima verosimilitud 𝜇 y 𝜆 utilizamos la función de verosimilitud dada por (1.6) y la fdd y función de supervivencia de la distribución, dadas por (3.2) y (3.4) respectivamente. Sea entonces la función de logverosimilitud 𝑙(𝜇, 𝜆) = 𝑛 ∑ 𝑖=1 𝛿𝑖{1 2log 𝜆 2𝜋−3 2log 𝑡𝑖−𝜆(𝑡𝑖−𝜇)2 2𝜇2𝑡𝑖} + 𝑛 ∑ 𝑖=1 (1 − 𝛿𝑖) log {Φ(𝐴𝑖) − exp (2𝜆 𝜇Φ(𝐵𝑖)}, (3.21) donde 𝐴𝑖= −√𝜆 𝑡𝑖(𝑡𝑖 𝜇− 1)y𝐵𝑖= −√𝜆 𝑡𝑖(𝑡𝑖 𝜇+ 1). A partir de esto, podemos calcular sus ecuaciones de verosimilitud. La primera de ellas vendría dada por 𝜕𝑙 𝜕𝜇(𝜇, 𝜆) = 𝑛 ∑ 𝑖=1 𝛿𝑖 𝜆(𝑡𝑖−𝜇) 𝜇3+ 𝑛 ∑ 𝑖=1 (1 − 𝛿𝑖)𝜕log 𝑆(𝑡𝑖) 𝜕𝜇 ,(3.22) donde 𝜕log 𝑆(𝑡) 𝜕𝜇 =1 𝜇2𝑆(𝑡){√𝜆𝑡𝜙[√𝜆 𝑡(𝑡 𝜇− 1)]+𝑒2𝜆∕𝜇[2𝜆Φ(−√𝜆 𝑡(𝑡 𝜇+ 1)) −√𝜆𝑡𝜙(√𝜆 𝑡(𝑡 𝜇+ 1))]}. Por otra parte, la segunda ecuación de verosimilitud es 𝜕𝑙 𝜕𝜆(𝜇, 𝜆) = 𝑛 ∑ 𝑖=1 𝛿𝑖(1 2𝜆−(𝑡𝑖−𝜇)2 2𝜇2𝑡𝑖)+ 𝑛 ∑ 𝑖=1 (1 − 𝛿𝑖)𝜕log 𝑆(𝑡𝑖) 𝜕𝜆 ,(3.23) donde 𝜕log 𝑆(𝑡) 𝜕𝜆 =1 𝑆(𝑡){−1 2√𝜆𝑡(𝑡 𝜇−1)𝜙(√𝜆 𝑡(𝑡 𝜇−1))−𝑒2𝜆∕𝜇[2 𝜇Φ(−√𝜆 𝑡(𝑡 𝜇+1)) −1 2√𝜆𝑡(𝑡 𝜇+ 1)𝜙(√𝜆 𝑡(𝑡 𝜇+ 1))]}. 3. distribución inversa gaussiana 47 Las ecuaciones de verosimilitud (3.22) y (3.23) nos darán la expresión de los estimadores de máxima verosimilitud 𝜇 y 𝜆utilizando algún método numérico, ya que no tienen una forma explícita. 3.3 Aplicación Para aplicar los resultados de este capítulo, consideramos un conjunto de datos obtenido del paquete "invGauss" de R (ver [8]), que originalmente fueron cedidos por el Radiation Therapy Oncology Group, en Estados Unidos. En la muestra, se estudia a 192 pacientes de cáncer de orofaringe, y se recoge el tiempo de vida en años de los pacientes desde su diagnóstico. Estos datos se encuentran en el Cuadro 3.1, donde los asteriscos (*) representan los datos censurados. Para esta aplicación vamos a tomar solamente los datos observados. Lamentablemente, no hemos encontrado un paquete de R que nos permita calcular los estimadores para datos censurados con la parametrización de la distribución utilizada en este trabajo, por lo que no ha sido posible realizar una aplicación sobre este tipo de datos. Aplicación. Muestra completa Tomamos los datos no censurados del Cuadro 3.1, formando así una muestra completa de 139 observaciones que suponemos que es una m.a.s. de una distribución Inversa Gaussiana. El histograma de esta muestra, en la Figura 3.4, nos permite observar una asimetría muy grande a la derecha. Procedemos a estimar los parámetros 𝜇y𝜆por el método de máxima verosimilitud, cuya expresión viene dada en (3.18) y (3.19). Resulta entonces 𝜇 = 1,067311  𝜆= 0,888078 Por tanto, podemos suponer que la muestra sigue aproximadamente una distribución 𝐼𝐺(1′07,0′89). Viendo su representación en la Figura 3.5, podemos tomar el 48 modelos de distribuciones obtenidas de la distribución normal Años de vida tras el diagnóstico 1,73 0,74 0,90 0,67 2,51 4,99* 1,75 0,64 0,70 0,50 2,92 1,13 0,59 0,89 1,32 0,67 4,29* 1,53 1,03 2,50 0,76 0,39 2,99 0,26 0,48 4,03* 1,44 0,47 1,58 0,61 0,46 4,29 0,70 0,37 1,11 4,10* 0,44 0,72 0,84 2,14 1,81 1,50 4,84* 1,02 4,08* 3,96* 0,20 4,41* 0,82 0,90 1,26 1,22 4,50* 1,35 0,76 2,51 0,62 0,35 4,31 1,54 1,01 2,21 0,53 0,75 3,77* 1,12 2,55 1,50 3,61* 3,61* 1,42 3,58* 0,63 2,09 0,47 3,99* 3,38* 1,49 2,19 4,00* 2,15 1,96 0,93 1,18 3,59* 0,96 0,56 3,34* 0,03 1,82 0,40 2,90* 1,31 2,90* 3,59* 1,91 0,31 0,84 0,04 0,36 0,81 0,80 1,49 2,98* 3,42* 0,40 1,99 0,85 1,64 2,73* 2,98* 1,05 2,55* 0,72 0,03 2,50* 0,24 1,44 1,46* 1,75 0,31 3,00* 0,47 2,58* 0,52 2,54* 2,52* 2,26* 0,27 0,27 2,56* 1,26 0,95 1,02 2,00* 0,65 1,62* 0,60 1,27 1,52 1,46 0,42 1,01 1,48 0,29 2,34* 2,25* 2,12 0,92 1,41 2,50* 2,07 2,18* 0,29 2,01* 1,64* 0,87 1,12 0,95 1,42 1,08 0,22 1,67 2,08* 0,94 0,89 0,70 2,06* 0,92 0,75 1,50* 0,31 0,50* 0,57 0,57 0,48 1,78* 1,84* 0,80 1,98* 1,36 0,76* 0,25* 0,58 0,10 0,35 1,22* 0,44 0,60 0,47 1,13* 0,75 Cuadro 3.1: Datos de pacientes de cáncer de orofaringe. ajuste como bueno. Podemos estimar la moda de la distribución utilizando (3.11), resultando  𝑡𝑚𝑜𝑑𝑎 = 0,2762016. Lo que podemos concluir a partir de la moda es que usualmente, ese es el tiempo que tarda en fallecer un paciente de este tipo de cáncer, ya que es el valor más repetido. Por otra parte, la varianza de la distribución se puede estimar usando (3.9), de donde obtenemos  𝜎2= 𝑉 𝑎𝑟(𝑇)=1,369056. 3. distribución inversa gaussiana 49 Cáncer de orofaringe Tiempo en años Frecuencia 01234 0 10 30 Figura 3.4: Histograma de datos no censurados del Cuadro 3.1. Cáncer de orofaringe Tiempo en años Frecuencia relativa 01234 0.0 0.4 0.8 Histograma relativo Densidad IG Figura 3.5: Histograma relativo de datos y densidad 𝐼𝐺(𝜇,  𝜆). La varianza tiene un valor bajo, ya que tal y como se ve en el histograma de la Figura 3.4, no hay mucha variabilidad de los datos. De (3.20), podemos calcular la matriz de información de Fisher estimada, que sería  𝐼=(0,7304302 0 0 0,6339686). 50 modelos de distribuciones obtenidas de la distribución normal Calculando su inversa, se obtiene la matriz de covarianzas estimada, resultando  𝑉=(1,369056 0 0 1,577365), donde los elementos de la diagonal son las varianzas estimadas de 𝜇 y 𝜆, respectivamente, lo cual nos da una idea de la variabilidad de estos estimadores. Es fácil estimar también el error estándar de cada estimador, tomando la raíz cuadrada de la diagonal de  𝑉. Se tiene entonces 𝜎𝜇= 1,170067, 𝜎𝜆= 1,255932. Son también valores bajos, lo que nos da la idea de que las estimaciones no se alejan mucho del valor real de los parámetros. 4 Distribución Birnbaum-Saunders 4.1 Introducción La distribución Normal es la más conocida y utilizada de todas las distribuciones estadísticas. Ya hemos visto en capítulos anteriores dos ejemplos de distribuciones de tiempo de vida relacionadas con la Normal. Otra de ellas es la distribución de Birnbaum-Saunders, que abreviaremos como BS. Es una distribución de dos parámetros que puede expresarse como una transformación monótona de la distribución Normal Estándar. Aunque ya se conocía esta distribución, en 1969 fueron Birnbaum y Saunders quienes le dieron una justificación natural física: modelar el tiempo de vida de un metal sujeto a estrés periódico. Por este motivo, la distribución recibe su nombre. Desde entonces se ha trabajado mucho en esta distribución, hallando así más campos de aplicación y nuevas propiedades de esta, tal y como veremos en este capítulo. Entre estas nuevas conclusiones, se encuentra una importante relación entre la distribución Inversa Gaussiana (IG) y la BS, ya que esta última se puede obtener como una aproximación de la primera. También se puede ver la distribución BS como una mezcla entre la IG y su recíproca, lo cual provoca que la distribución BS herede ciertas propiedades ya conocidas de la IG y es útil para estudiar su inferencia. La función de densidad y la función de hazard de esta distribución también son objetos interesantes de estudio. Se ha observado que la fdd de la BS es unimodal, al igual que su función de hazard. La función de hazard de las distribuciones IG y BS tienen comportamientos parecidos, ya que cuando 𝑡→∞, tras llegar a su máximo ambas decrecen hacia una constante, mientras que la de la distribución log-normal tiende a cero. En este capítulo estudiaremos también los estimadores de máxima verosimilitud de los parámetros de forma y escala (que notaremos 𝛼y𝛽respectivamente) basados en muestras completas y censuradas, y veremos que estos no pueden obtenerse de forma explícita, sino resolviendo una ecuación no lineal. 58 modelos de distribuciones obtenidas de la distribución normal Como 𝛽es el parámetro de escala, el punto de cambio de la función de hazard en (4.13), que notaremos como 𝑐𝛼,𝛽 se puede calcular de forma inmediata como 𝑐𝛼,𝛽 =𝛽𝑐𝛼. 4.4 Momentos y momentos inversos Proposición 4.4. Sea 𝑇una v.a. tal que 𝑇∼𝐵𝑆(𝛼, 𝛽)y sea un entero r no negativo. Entonces, el momento no central de orden r vendrá dado por 𝐸(𝑇𝑟) = 𝛽𝑟 𝑟 ∑ 𝑗=0 (2𝑟 2𝑗)𝑗 ∑ 𝑘=0 (𝑗 𝑘)(2𝑟− 2𝑗+ 2𝑘)! 2𝑟−𝑗+𝑘(𝑟−𝑗+𝑘)!(𝛼 2)2𝑟−2𝑗+2𝑘.(4.15) Demostración. Partiendo de (4.8), obtenemos la expresión 𝐸[(𝑇 𝛽)𝑟]=𝐸((𝛼𝑍 2+√(𝛼𝑍 2)2 + 1)2𝑟). Si aplicamos a esta expresión el binomio de Newton, esto es, (𝑎+𝑏)𝑚= 𝑚 ∑ 𝑘=0 (𝑚 𝑘)𝑎𝑚−𝑘𝑏𝑘, se tiene que 𝐸[(𝑇 𝛽)𝑟]= 2𝑟 ∑ 𝑘=0 (2𝑟 𝑘)𝐸[((𝛼𝑍 2)2 + 1)𝑘∕2(𝛼𝑍 2)2𝑟−𝑘]. Observamos que 𝐸{[(𝛼𝑍 2)2 + 1]𝑡(𝛼𝑍 2)𝑠}= 0 si 𝑠es impar. Por tanto, 𝐸[(𝑇 𝛽)𝑟]= 𝑟 ∑ 𝑘=0 (2𝑟 2𝑘)𝐸{[(𝛼𝑍 2)2 + 1]𝑘(𝛼𝑍 2)2(𝑟−𝑘)}. Usando de nuevo el binomio de Newton: [(𝛼𝑍 2)2 + 1]𝑗 = 𝑗 ∑ 𝑘=0 (𝑗 𝑘)(𝛼𝑍 2)2𝑘 , 4. distribución birnbaum-saunders 59 de donde 𝐸[(𝑇 𝛽)𝑟]= 𝑟 ∑ 𝑗=0 (2𝑟 2𝑗)𝐸[𝑗 ∑ 𝑘=0 (𝑗 𝑘)(𝛼𝑍 2)2𝑘(𝛼𝑍 2)2(𝑟−𝑗)] = 𝑟 ∑ 𝑗=0 (2𝑟 2𝑗)𝑗 ∑ 𝑘=0 (𝛼 2)2(𝑟−𝑗+𝑘)𝐸[𝑍2(𝑟−𝑗+𝑘)]. Como 𝑍∼𝑁(0,1), entonces 𝐸(𝑍2(𝑟−𝑗+𝑘))=(2(𝑟− +𝑘))! 2𝑟−𝑗+𝑘(𝑟−𝑗+𝑘)!, de donde obtenemos finalmente que 𝐸[(𝑇 𝛽)𝑟]= 𝑟 ∑ 𝑗=0 (2𝑟 2𝑗)𝑗 ∑ 𝑘=0 (𝑗 𝑘)(2𝑟− 2𝑗+ 2𝑘)! 2𝑟−𝑗+𝑘(𝑟−𝑗+𝑘)!(𝛼 2)2𝑟−2𝑗+2𝑘. La expresión de los momentos no centrales nos permite calcular fácilmente el valor de la esperanza y varianza de la distribución, así como los coeficientes de asimetría y kurtosis. Proposición 4.5. Sea 𝑇∼𝐵𝑆(𝛼, 𝛽). Entonces, 1. Esperanza 𝐸(𝑇) = 𝛽(1 + 𝛼2 2)(4.16) 2. Varianza 𝑉 𝑎𝑟(𝑇) = 𝛽2𝛼2(1 + 5 4𝛼2)(4.17) 3. Coeficiente de asimetría 𝛾1(𝑇) = 16𝛼2(11𝛼2+ 6)2 (5𝛼2+ 4)3(4.18) 4. Coeficiente de kurtosis 𝛽2(𝑇) = 3 + 6𝛼2(93𝛼2+ 40) (5𝛼2+ 4)2(4.19) Demostración. 60 modelos de distribuciones obtenidas de la distribución normal 1. Tomamos 𝑟= 1 en (4.15): 𝐸(𝑇 𝛽)= 1 ∑ 𝑗=0 (2 2𝑗)𝑗 ∑ 𝑘=0 (𝑗 𝑘)(2 − 2𝑗+ 2𝑘)! 21−𝑗+𝑘(1 − 𝑗+𝑘)!(𝛼 2)2−2𝑗+2𝑘 =(2 0)2! 1!2(𝛼 2)2 +(2 2)1 ∑ 𝑘=0 (1 𝑘)(2𝑘)! 2𝑘𝑘!(𝛼 2)2𝑘 = 1 + 𝛼2 2, de donde 𝐸(𝑇) = 𝛽(1 + 𝛼2 2). 2. Sabemos que 𝑉 𝑎𝑟(𝑇) = 𝐸(𝑇2) − 𝐸(𝑇)2, luego calculamos 𝐸(𝑇2)usando (4.15): 𝐸[(𝑇 𝛽)2]= 2 ∑ 𝑗=0 (4 2𝑗)𝑗 ∑ 𝑘=0 (𝑗 𝑘)(4 − 2𝑗+ 2𝑘)! 22−𝑗+𝑘(2 − 𝑗+𝑘)!(𝛼 2)4−2𝑗+2𝑘 = 3(𝛼 2)4 + 6(𝛼 2)2 + 18(𝛼 2)4 + 1 + 2(𝛼 2)2 + 3(𝛼 2)4 =3 2𝛼4+ 2𝛼2+ 1. Entonces, 𝐸(𝑇2) = 𝛽2(3 2𝛼4+ 2𝛼2+ 1), luego 𝑉 𝑎𝑟(𝑇) = 𝛽2(3 2𝛼4+ 2𝛼2+ 1)−(𝛽(1 + 𝛼2 2))2 =𝛽2𝛼2(1 + 5 4𝛼2). 3. El coeficiente de asimetría se define de forma general como 𝛾1(𝑇) = 𝜇3 (𝜇2)3∕2 ,(4.20) donde 𝜇𝑘=𝐸[(𝑇−𝐸(𝑇))𝑘]es el momento centrado de orden k. Entonces, solo necesitamos hallar 𝜇3, ya que se tiene 𝜇2=𝑉 𝑎𝑟(𝑇). 𝐸[(𝑇 𝛽)3]= 3 ∑ 𝑗=0 (6 2𝑗)𝑗 ∑ 𝑘=0 (𝑗 𝑘)(4 − 2𝑗+ 2𝑘)! 22−𝑗+𝑘(2 − 𝑗+𝑘)!(𝛼 2)4−2𝑗+2𝑘 = 1 + 9 2𝛼2+ 9𝛼4+15 2𝛼6, luego 𝐸(𝑇3) = 𝛽3 2(2+9𝛼2+ 18𝛼4+ 15𝛼6). 4. distribución birnbaum-saunders 61 𝜇3se puede expresar en función de momentos no centrales como sigue 𝜇3=𝐸[𝑇3+ 3𝑇 𝐸(𝑇)2− 3𝑇2𝐸(𝑇) − 𝐸(𝑇)3], y aplicando propiedades conocidas de la esperanza, resulta 𝜇3=𝐸(𝑇3)+2𝐸(𝑇)3− 3𝐸(𝑇)𝐸(𝑇2). Aplicando estos resultados a la definición de 𝛾1en (4.20), se tiene el resultado. 4. El coeficiente de kurtosis se expresa como 𝛽2(𝑇) = 𝜇4 𝜇2 2 .(4.21) Siguiendo el mismo proceso de antes, calculamos 𝐸(𝑇4). 𝐸[(𝑇 𝛽)4]= 4 ∑ 𝑗=0 (6 2𝑗)𝑗 ∑ 𝑘=0 (𝑗 𝑘)(4 − 2𝑗+ 2𝑘)! 22−𝑗+𝑘(2 − 𝑗+𝑘)!(𝛼 2)4−2𝑗+2𝑘 = 1 + 8𝛼2+ 30𝛼4+ 60𝛼6+105 2𝛼8. Entonces, 𝐸(𝑇4) = 𝛽4 2(2 + 16𝛼2+ 60𝛼4+ 120𝛼6+ 105𝛼8). Por otro lado, se cumple 𝜇4(𝑇) = 𝐸[𝑇4+ 6𝑇2𝐸(𝑇)2− 4𝑇3𝐸(𝑇)−4𝑇 𝐸(𝑇)3+𝐸(𝑇)4], de donde 𝜇4(𝑇) = 𝐸(𝑇4)+6𝐸(𝑇2)𝐸(𝑇)2− 4𝐸(𝑇3)𝐸(𝑇)−3𝐸(𝑇)4. Entonces, sustituyendo los valores necesarios en (4.21), finaliza la prueba. Para obtener los momentos inversos de la BS, vamos a utilizar la propiedad descrita a continuación, que a su vez nos será útil posteriormente para la inferencia de la distribución. Proposición 4.6. Sea 𝑇una v.a. tal que 𝑇∼𝐵𝑆(𝛼, 𝛽).Entonces, se cumple que tanto 𝑇∕𝛽como 𝛽∕𝑇siguen la misma distribución, esto es, 𝑇 𝛽∼𝐵𝑆(𝛼, 1),𝛽 𝑇∼𝐵𝑆(𝛼, 1). 62 modelos de distribuciones obtenidas de la distribución normal Demostración. Es claro que al ser 𝛽el parámetro de escala y ser 𝑇∼𝐵(𝛼, 𝛽),entonces se tiene que 𝑇∕𝛽∼𝐵(𝛼, 1). Para la segunda parte de la prueba, hacemos el cambio de variable 𝑅=𝛽∕𝑇. Entonces, 𝐹𝑅(𝑟) = 𝑃(𝛽 𝑇≤𝑟)=𝑃(𝛽 𝑟≤𝑇)= 1 − Φ[1 𝛼{(𝛽∕𝑟 𝛽)1∕2 −(𝛽 𝛽∕𝑟)1∕2}]= = Φ[1 𝛼{(𝑟 1)1∕2 −(1 𝑟)1∕2}], luego 𝑅=𝛽 𝑇∼𝐵(𝛼, 1). A partir de esta última proposición, es fácil ver que 𝐸(𝑇−𝑟) = 𝐸(𝑇𝑟) 𝛽2𝑟,(4.22) relación que nos permitirá calcular de forma simple los momentos no centrados de la v.a. 𝑇−1.En la siguiente proposición se darán la media y la varianza de la misma. Proposición 4.7. Dada una v.a. 𝑇∼𝐵𝑆(𝛼, 𝛽), se tiene que la media y la varianza de 𝑇−1 vienen dadas por 1. Media 𝐸(1 𝑇)=1 𝛽(1 + 𝛼2 2)(4.23) 2. Varianza 𝑉 𝑎𝑟(1 𝑇)=𝛼2 𝛽2(5 4𝛼2+ 1).(4.24) Demostración. 1. Por (4.22), se tiene que 𝐸(𝑇−1) = 𝐸(𝑇) 𝛽2, donde el valor de 𝐸(𝑇)viene dado por (4.16). Entonces, 𝐸(𝑇−1) = 1 𝛽2𝛽(𝑎+𝛼2 2)=1 𝛽(1 + 𝛼2 2). 4. distribución birnbaum-saunders 63 2. De igual manera, por (4.22), 𝐸(𝑇−2) = 𝐸(𝑇2) 𝛽4=1 𝛽2(3 2𝛼4+ 2𝛼2+ 1). Entonces, como 𝑉 𝑎𝑟(𝑇−1) = 𝐸(𝑇−2) − 𝐸(𝑇−1)2, se tiene finalmente 𝑉 𝑎𝑟(𝑇−1) = 1 𝛽2(3 2𝛼4+ 2𝛼2+ 1)−[1 𝛽(1 + 𝛼2 2)]2, de donde se tiene el resultado. 4.5 Relación con la distribución Inversa Gaussiana Consideramos dos variables aleatorias 𝑋1y𝑋2tales que 𝑋1∼𝐼𝐺(𝜇, 𝜆)y𝑋−1 2∼ 𝐼𝐺(𝜇−1, 𝜆𝜇2),donde 𝜇 > 0y𝜆 > 0. Definimos la v.a. 𝑇como 𝑇={𝑋1con probabilidad 1∕2 𝑋2con probabilidad 1∕2 Entondes, la fdd de 𝑇vendrá dada por 𝑓𝑇(𝑡;𝜇, 𝜆) = 1 2𝑓𝑋1(𝑡;𝜇, 𝜆) + 1 2𝑓𝑋2(𝑡;𝜇, 𝜆), 𝑡 > 0,(4.25) donde 𝑓𝑋1y𝑓𝑋2son las funciones de densidad de 𝑋1y𝑋2, respectivamente. La primera viene dada por la expresión (3.2), mientras que la segunda será de la forma 𝑓𝑋2(𝑡;𝜇, 𝜆) = 𝑡𝑓𝑋1(𝑡;𝜇, 𝜆)∕𝜇. Tomando 𝛼=√𝜇∕𝜆y𝛽=𝜇, obtenemos que la v.a. 𝑇definida anteriormente sigue una distribución 𝐵𝑆(𝛼, 𝛽).Además, podemos obtener transformaciones de muchas de las propiedades de la 𝐼𝐺 en las mismas de la 𝐵𝑆 de forma inmediata, usando la relación vista en (4.25). Un ejemplo de esto es que podemos usarla para obtener la FGM de la distribución 𝐵𝑆 a partir de las FGMs de la 𝐼𝐺. Veámoslo: Si 𝑇∼𝐵𝑆(𝛼, 𝛽), su FGM será de la forma 𝑀𝑇(𝑟) = 1 2[𝑀𝑋1(𝑟) + 𝑀𝑋2(𝑟)], 64 modelos de distribuciones obtenidas de la distribución normal siendo 𝑀𝑋1y𝑀𝑋2son las FGM de 𝑋1y𝑋2, respectivamente, las cuales tienen la siguiente expresión 𝑀𝑋1(𝑟) = exp {𝜆 𝜇(1 − (1 − 2𝜇2𝑟 𝜆)1∕2)}, 𝑟 ∈ℝ, 𝑀𝑋2(𝑟) = exp {𝜆 𝜇(1 − (1 − 2𝜇2𝑟 𝜆)1∕2)}(1 − 2𝜇2𝑡 𝜆)−1∕2, 𝑟 ∈ℝ. De manera inmediata, se obtiene 𝑀𝑇(𝑟) = 1 2exp {𝜆 𝜇(1 − (1 − 2𝜇2𝑟 𝜆)1∕2)}[1 + (1 − 2𝜇2𝑡 𝜆)−1∕2]. 4.6 Inferencia En esta sección, estudiaremos distintos métodos de estimación de los parámetros 𝛼y𝛽basados primero en una muestra aleatoria completa 𝑇1, ..., 𝑇𝑛de la v.a. 𝑇∼𝐵𝑆(𝛼, 𝛽)con fdd 𝑓𝑇y cuyas observaciones son {𝑡1, ..., 𝑡𝑛}.Posteriormente, estudiaremos el cálculo de estos estimadores para el caso de muestras censuradas. 4.6.1 Método de los momentos Los estimadores de 𝛼y𝛽por el método de los momentos se obtienen igualando la media muestral y la varianza muestral a 𝐸(𝑇)y𝑉 𝑎𝑟(𝑇), respectivamente, cuya expresión se vio en (4.16) y (4.17). Denotamos la media y la varianza muestrales como 𝑠y𝑣, respectivamente. Entonces, los estimadores por el método de los momentos de 𝛼y𝛽se obtienen resolviendo el siguiente sistema: 𝑠=𝛽 2(𝛼2+ 2) 𝑣=𝛽2 4(5𝛼4+ 4𝛼2) De aquí, se tiene de forma inmediata que el estimador por el método de los momentos de 𝛼se puede obtener como raíz de la siguiente ecuación no lineal 𝛼4(5𝑠2−𝑣)+4𝛼2(𝑠2− 4) − 4𝑣= 0. 4. distribución birnbaum-saunders 65 Entonces, si el coeficiente de variación muestral, 𝜎∕𝑡, es menor que √5, el estimador de 𝛼tendrá la siguiente expresión 𝛼𝑀=[−2(𝑠2− 4) + 2√(𝑠2− 4)2+𝑣(5𝑠2−𝑣) 5𝑠2−𝑣]1∕2 ,(4.26) mientras que el de 𝛽será entonces  𝛽𝑀=2𝑠 𝛼2 𝑀+ 2.(4.27) Sin embargo, si el coeficiente de variación muestral 𝑣∕𝑠fuera mayor que √5, el estimador de 𝛼por el método de los momentos podría no existir, y por tanto, ocurriría lo mismo con el de 𝛽. 4.6.2 Método de los momentos modificado Como veremos más tarde, los estimadores MV de 𝛼y𝛽no pueden obtenerse de forma explícita, y a su vez, los estimadores por el método de los momentos no tienen por qué existir. Este es el motivo por el que se propone el método de los momentos modificado (MM) para calcular los estimadores de 𝛼y𝛽. Para calcularlos, usamos el hecho de que si 𝑇∼𝐵𝑆(𝛼, 𝛽), entonces 𝑇−1 ∼𝐵𝑆(𝛼, 𝛽−1)visto en la Proposición 4.6. Igualamos la media muestral aritmética 𝑠y la media muestral armónica 𝑟a sus correspondientes expresiones poblacionales en (4.16) y (4.23), respectivamente, donde 𝑠=1 𝑛 𝑛 ∑ 𝑖=1 𝑡𝑖y𝑟=[1 𝑛 𝑛 ∑ 𝑖=1 𝑡−1 𝑖]−1. Resultan entonces las ecuaciones 𝑠=𝛽(1 + 1 2𝛼2), 𝑟−1 =𝛽−1(1 + 1 2𝛼2). Así, los estimadores MM de 𝛼y𝛽poseen la siguiente expresión 𝛼𝑀𝑀 =(2[(𝑠 𝑟)1∕2 − 1])1∕2 y 𝛽𝑀𝑀 =√𝑠𝑟. (4.28) 66 modelos de distribuciones obtenidas de la distribución normal Los estimadores MM siempre existen y tienen forma explícita. Además, utilizando el TCL en 1 𝑛∑𝑛 𝑖=1 𝑇𝑖y1 𝑛∑𝑛 𝑖=1 𝑇−1 𝑖, se puede probar que 𝛼𝑀𝑀 y 𝛽𝑀𝑀 siguen una distribución normal bivariante (ver [6]), de manera que (𝛼𝑀𝑀  𝛽𝑀𝑀 )∼𝑁[(𝛼 𝛽),Σ], donde Σ = ⎛⎜⎜⎜⎝ 𝛼2 2𝑛0 0𝛼𝛽2 𝑛(1+ 3 4𝛼2 (1+ 1 2𝛼2)2)⎞⎟⎟⎟⎠ .(4.29) es la matriz de covarianzas asintótica. Observación 4.1.Los estimadores MM son asintóticamente independientes. 4.6.3 Método de máxima verosimilitud La función de log-verosimilitud de (𝛼, 𝛽)basada en las observaciones 𝑡1, ..., 𝑡𝑛es la siguiente: 𝑙(𝛼, 𝛽) = 𝑐1+𝑛 𝛼2−1 2𝛼2 𝑛 ∑ 𝑖=1 (𝑡𝑖 𝛽+𝛽 𝑡𝑖)−𝑛log(𝛼) − 𝑛 2log(𝛽) + 𝑛 ∑ 𝑖=1 log(𝑡𝑖+𝛽),(4.30) donde 𝑐1es una constante que no depende de 𝛼ni de 𝛽. Para maximizar (4.30), necesitamos calcular sus primeras derivadas respecto de 𝛼y de 𝛽, obteniendo así las siguientes expresiones: 𝜕𝑙(𝛼, 𝛽) 𝜕𝛼 = −2𝑛 𝛼3+1 𝛼3 𝑛 ∑ 𝑖=1 (𝑡𝑖 𝛽+𝛽 𝑡𝑖)−𝑛 𝛼, 𝜕𝑙(𝛼, 𝛽) 𝜕𝛽 =1 2𝛼2 𝑛 ∑ 𝑖=1 (𝑡𝑖 𝛽2−1 𝑡𝑖)−𝑛 2𝛽+ 𝑛 ∑ 𝑖=1 1 𝑡𝑖+𝛽. A continuación, igualamos ambas expresiones a cero, y evaluamos en 𝛼=𝛼𝑀𝑉 y 𝛽= 𝛽𝑀𝑉 .Para el caso de 𝛼, obtenemos 𝛼𝑀𝑉 =(1 𝑛 𝑛 ∑ 𝑖=1 (𝑡𝑖  𝛽𝑀𝑉 +  𝛽𝑀𝑉 𝑡𝑖))1∕2 ,(4.31) 4. distribución birnbaum-saunders 67 que también podemos expresar como 𝛼𝑀𝑉 =(𝑠  𝛽𝑀𝑉 +  𝛽𝑀𝑉 𝑟− 2)1∕2 ,(4.32) donde 𝑠y𝑟son las medias aritmética y armónica, respectivamente, de las observaciones 𝑡1, ..., 𝑡𝑛. Por otro lado, al repetir el procedimiento con 𝛽, y sustituyendo 𝛼𝑀𝑉 según la expresión (4.32), obtenemos la siguiente ecuación no lineal: 𝛽2−𝛽(2𝑟+𝐾(𝛽)) + 𝑟(𝑠+𝐾(𝛽)) = 0,(4.33) donde hemos llamado 𝐾(𝑥)a 𝐾(𝑥) = [1 𝑛 𝑛 ∑ 𝑖=1 (𝑥+𝑡𝑖)−1]−1 para 𝑥≥0. El estimador de 𝛽, 𝛽𝑀𝑉 , será la raíz positiva de (4.33). Birnbaum y Saunders [3] ya probaron que  𝛽es la única raíz positiva de (4.33), y además, 𝑟 <  𝛽 < 𝑠. Para poder resolver la ecuación (4.33), se necesita un método numérico iterativo. Una vez calculado  𝛽, se puede sustituir en (4.32), obteniendo así el estimador 𝛼. El método iterativo propuesto por Leiva [13] se compone de las siguientes ecuaciones iterativas para 𝑚= 0,1,2,…. 𝛼(𝑚+1) =1 √𝑛(𝑛 ∑ 𝑖=1 (𝑡𝑖  𝛽(𝑚)+ 𝛽(𝑚) 𝑡𝑖 − 2))1∕2 ,(4.34)  𝛽(𝑚+1) =√ √ √ √ √ 1 2(𝛼(𝑚))2∑𝑛 𝑖=1 𝑡𝑖 1 2(𝛼(𝑚))2∑𝑛 𝑖=1 1 𝑡𝑖 −∑𝑛 𝑖=1 1 𝑡𝑖+ 𝛽(𝑚)+𝑛 2 𝛽(𝑚) .(4.35) Son necesarios unos valores iniciales 𝛼(0) y 𝛽(0). Birnbaum y Saunders propusieron el valor inicial  𝛽(0) = (𝑠𝑟)1∕2, a partir del cual se puede obtener el valor inicial 𝛼(0) correspondiente mediante la expresión (4.32). Los estimadores 𝛼 y 𝛽calculados por máxima verosimilitud (MV) son estimadores consistentes, para más detalles ver [6]. Además, estudios más recientes han demostrado que para un tamaño muestral 𝑛= 1, los estimadores MV de 𝛼y𝛽no existen, 74 modelos de distribuciones obtenidas de la distribución normal Comparando este resultado con las densidades ajustadas de la muestra completa en la Figura 4.3, vemos que el ajuste en este caso es peor, aunque podemos tomarlo igualmente como válido. Histograma de datos censurados Fatiga Frecuencia relativa 100 150 200 0.000 0.010 Histograma relativo Densidad BS Figura 4.5: Histograma de datos de fatiga censurados y densidad 𝐵𝑆(𝛼,  𝛽). Procediendo como en la aplicación anterior, calculamos los cuantiles de interés de la distribución obtenida. Cuantil de orden 0,05  𝑡0,05 = 96,17202 Cuantil de orden 0,10  𝑡0,1= 103,60073 Mediana  𝑡0,5= 𝛽= 134,99664 Cuantil de orden 0,90  𝑡0,9= 175,90700 Cuantil de orden 0,95  𝑡0,95 = 189,49476 En vista de estos resultados, podemos esperar que las láminas de aluminio puedan resistir unos 135.000 ciclos bajo estrés antes de fallar. Podemos garantizar con un 5 % de significación que este tipo de producto puede soportar 96.000 ciclos sin romperse, 4. distribución birnbaum-saunders 75 mientras que no podrán durar más de 190.000. Trabajar con datos censurados es más complejo que en el caso completo, por lo que no podemos extender mucho más el estudio de esta aplicación. Bibliografía [1] Balakrishnan, N., and Kundu, D. Birnbaum-saunders distribution: A review of models, analysis, and applications. Appl Stochastic Models Bus Ind. 35 (2019), 4–49. [2] Barros, M., Leiva, V., and Paula, G. gbs: Package for Generalized Birnbaum- Saunders Distributions, 2012. R package version 1.0. [3] Birnbaum, Z., and Saunders, S. Estimation for a family of life distributions with applications to fatigue. Journal of Applied Probability 6, 2 (1969), 328–347. doi:10.2307/3212004. [4] Chang, D., and Tang, L. Reliability bounds and critical time for the birnbaumsaunders distribution. IEEE Trans. Device Mater. Reliab. 42 (1993), 464–469. [5] Chhikara, R., and Folks, J. The inverse gaussian distribution as a lifetime model. Technometrics 19, 4 (1977), 461–468. [6] Cox, D., and Hinkley, D. Theoretical Statistics. Chapman and Hall, 1974. [7] Cristóbal, J. Inferencia estadística. Ed. Prensas de la Universidad de Zaragoza, 2000. [8] Gjessing, H. K. invGauss: Threshold regression that fits the (randomized drift) inverse Gaussian distribution to survival data., 2014. R package version 1.1. [9] Henze, N., and Klar, B. Goodness-of-fit tests for the inverse gaussian distribution based on the empirical laplace transform. Annals of the Institute of Statistical Mathematics 54 (2002), 425—-444. [10] Kundu, D., Kannan, N., and Balakrishnan, N. On the hazard function of birnbaum–saunders distribution and associated inference. Comput Stat Data Anal 52, 5 (2008), 2692–2702. 78 modelos de distribuciones obtenidas de la distribución normal [11] Lawless, J. Statistical Models and Methods for Lifetime Data. Ed. Wiley- Inters- cience, 2003. [12] Lee, E., and Wang, J. Statistical Methods for Survival Data Analysis. Ed. John Wiley Sons, 2003. [13] Leiva, V. The Birnbaum-Saunders Distribution. Ed. Academic Press, 2016. [14] Lemeshko, B., et al. Inverse Gaussian Model and Its Applications in Reliability and Survival Analysis. 10 2010, pp. 433–453. [15] Overduin, S. Use of the lognormal distribution for survival data: Inference and robustness. Master of science, Simon Fraser University, 2004. [16] R Core Team.R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2019. [17] Rohatgi, V., and Ehsanes Saleh, A. K. M. An Introduction to Probability and Statistics. Ed. John Wiley Sons, 2011. [18] Therneau, T. M. A Package for Survival Analysis in S, 2015. version 2.38. [19] Tweedie, M. C. K. Statistical properties of inverse gaussian distributions. i. The Annals of Mathematical Statistics 28, 2 (1957), 362–377.