Un recorrido por las distribuciones notables en la Inferencia Estadística con aplicaciones
Abstract
Las distribuciones de probabilidad son una de las herramientas fundamentales dentro de la Inferencia Estadística. El objetivo que se persigue en este trabajo es llevar a cabo un paseo por algunas de las distribuciones de probabilidad unidimensionales más notables, tanto continuas como discretas, estudiándolas desde el punto de vista de sus aplicaciones y mostrando algunas de las relaciones que hay entre ellas a través de resultados como el Teorema Central del Límite. Además, se lleva a cabo una pequeña introducción a la simulación de distribuciones de variables aleatorias y se muestran algunos de los principales algoritmos incluyendo algunos ejemplos prácticos. Por último, para enfatizar la idea de que las distribuciones de probabilidad no son solo una construcción analítica sino que muchos fenómenos de la realidad se ajustan a esas leyes se estudia la bondad de ajuste de dos distribuciones, una continua y una discreta, a dos conjuntos de datos no simulados.
Full text
Traballo Fin de Grao UN RECORRIDO POR LAS DISTRIBUCIONES NOTABLES EN LA INFERENCIA ESTADÍSTICA CON APLICACIONES Pablo Carballo Fraguas 2021/2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA
GRAO DE MATEMÁTICAS Traballo Fin de Grao UN RECORRIDO POR LAS DISTRIBUCIONES NOTABLES EN LA INFERENCIA ESTADÍSTICA CON APLICACIONES Pablo Carballo Fraguas julio, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA
Trabajo propuesto Área de Coñecemento: Estadística e Investigación Operativa Título: Un recorrido por las distribuciones notables en la Inferencia Estadística con aplicaciones. Breve descrición do contido Este trabajo tiene por objetivo realizar un mapa completo de las distribuciones notables discretas y continuas y sus relaciones en la Inferencia Estadística. La motivación del estudio surge por las importantes conexiones que existen entre las distribuciones notables, su motivación para la modelización en varios contextos de la vida real y la posibilidad de usar estas relaciones para la simulación articial de las distribuciones diversas, discretas o continuas. El trabajo se estructurará como sigue: 1) Un recorrido por las distribuciones notables discretas. 2) Un recorrido por las distribuciones notables continuas. 3) Simulación articial de las diversas distribuciones. 4) Aplicaciones con datos simulados y/o reales. iii
Índice Resumen vii Introducción ix 1 Un recorrido por las distribuciones notables discretas 1 1.1 Uniformediscreta.................................... 1 1.2 Bernoulli, Binomial e Hipergeométrica . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Geométrica ....................................... 9 1.4 Poisson ......................................... 10 1.5 Benford ......................................... 14 1.6 Zipf ........................................... 17 2 Un recorrido por las distribuciones notables continuas 21 2.1 Weibull ......................................... 22 2.1.1 Exponencial................................... 24 2.1.2 Weibull-Geométrica .............................. 25 2.1.3 Weibull-Geométrica complementaria . . . . . . . . . . . . . . . . . . . . . 27 2.2 Distribución de Valores Extremos . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3 Normal ......................................... 35 2.3.1 Ji-cuadrado .................................. 38 2.3.2 t deStudent................................... 39 2.3.3 F deSnédecor ................................. 41 3 Simulación articial de distribuciones 43 3.1 Introducción a la simulación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2 Simulación de distribuciones continuas . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2.1 Métododeinversión .............................. 44 3.2.2 Método de aceptación/rechazo . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2.3 Algoritmos para distribuciones especícas . . . . . . . . . . . . . . . . . . 47 3.3 Simulación de distribuciones discretas . . . . . . . . . . . . . . . . . . . . . . . . 48 3.3.1 Método de la transformación cuantil . . . . . . . . . . . . . . . . . . . . . 48 v
vi ÍNDICE 3.3.2 Métodos de truncamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3.3 Algoritmos para distribuciones especícas . . . . . . . . . . . . . . . . . . 51 4 Aplicaciones con datos reales 53 Anexos I Resultados auxiliares 57 II Figuras 59 IIITablas 69 IV Código de 91 IV.1 Código relativo a los test del Capítulo 4 . . . . . . . . . . . . . . . . . . . . . . . 91 IV.2 Código relativo a guras del Anexo II . . . . . . . . . . . . . . . . . . . . . . . . 95 IV.3 Código relativo a las tablas del Apéndice III . . . . . . . . . . . . . . . . . . . . . 104 Bibliografía 125
Resumen Las distribuciones de probabilidad son una de las herramientas fundamentales dentro de la Inferencia Estadística. El objetivo que se persigue en este trabajo es llevar a cabo un paseo por algunas de las distribuciones de probabilidad unidimensionales más notables, tanto continuas como discretas, estudiándolas desde el punto de vista de sus aplicaciones y mostrando algunas de las relaciones que hay entre ellas a través de resultados como el Teorema Central del Límite. Además, se lleva a cabo una pequeña introducción a la simulación de distribuciones de variables aleatorias y se muestran algunos de los principales algoritmos incluyendo algunos ejemplos prácticos. Por último, para enfatizar la idea de que las distribuciones de probabilidad no son solo una construcción analítica sino que muchos fenómenos de la realidad se ajustan a esas leyes se estudia la bondad de ajuste de dos distribuciones, una continua y una discreta, a dos conjuntos de datos no simulados. Abstract Probability distributions are one of the most important tools in Statistical Inference. The aim of this work is to take a journey through some of the most relevant one-dimensional probability distributions, both continuous and discrete, by studying them from the point of view of their applications and showing some of the relationships between them through results such as the Central Limit Theorem. In addition, a short introduction to the simulation of distributions of random variables is given and some of the main algorithms are shown, including some practical examples. Finally, to emphasize the idea that probability distributions are not only an analytical construct but that many phenomena in reality conform to these laws, we study the goodness of t of two distributions, a continuous and a discrete one, to two non-simulated data sets. vii
2 1. Un recorrido por las distribuciones notables discretas suceso una probabilidad igual a 1/n . Lo verdaderamente llamativo de esta distribución es que los paramétros estadísticos descriptivos que se pueden obtener a partir de una muestra aleatoria simple de datos procedentes de esta distribución (media muestral, varianza muestra...) coinciden con los correspondientes poblacionales. Así, por ejemplo, dada una muestra de datos x1, x2, ..., xn procedente de una distribución Uniforme Discreta se tiene que 1 X=1 n n X i=1 xi=1 n n(n+ 1) 2=n+ 1 2=µ, S2 X=1 n n X i=1 (xi−x)2=1 n n X i=1 x2 i−2x2+x2=(n+ 1)(2n+ 1) 6−(n+ 1)2 4=n2−1 12 =σ2. Por último, en cuanto a las aplicaciones de esta distribución, cabe mencionar que los ámbitos en que se utiliza están ligados a la condición de que los sucesos sean equiprobables, algo que no siempre se puede asumir porque puede haber factores que provoquen que la probabilidad no sea repartida a partes iguales. Sin embargo, hay situaciones donde se acostumbra a hacer esta suposición como, por ejemplo, en el lanzamiento de un dado (supuesto que no esté cargado). 1.2. Bernoulli, Binomial e Hipergeométrica La distribución de Bernoulli es una de las distribuciones discretas más simples que existen y se emplea para modelar un experimento aleatorio binario en el cual se tienen dos sucesos posibles que son mutuamente excluyentes, es decir, solamente hay dos posibilidades: éxito, que numéricamente se asociará con el 1; o fracaso, que se asociará con el 0. Así, una variable aleatoria que siga esta distribución toma únicamente los valores 0 o 1 . Además, se denota por p , con p∈(0,1) , a la probabilidad de que tome el valor 1 (probabilidad de éxito), mientras que la probabilidad de que tome el valor 0 (probabilidad de fracaso) se denota por q , con q= 1 −p 2 . Introducida la distribución de Bernoulli, se van a presentar de manera conjunta las distribuciones Binomial , Hipergeométrica , Binomial Negativa e Hipergeométrica Negativa , en lugar de cada una por separado, atendiendo a los criterios que se emplean en [Miller and Fridell, 2007] y [Chae, 1993] para claricar la relación que existe entre ellas. La motivación de dicho estudio se debe a una serie de relaciones que unen a estas distribuciones y las cuales se reejan en la Figura 2. Para mostrar con claridad estas relaciones se presentará una misma situación en la cual, 1 Para las cuentas posteriores se emplean las fórmulas de la Proposición I.1 y la Proposición I.2. 2 Para profundizar acerca de las principales características de la distribución de Bernoulli véase la Tabla 2.
1.2. Bernoulli, Binomial e Hipergeométrica 3 en función del enfoque que se dé, aparecerá una de las cuatro distribuciones. Supóngase que se dispone de una urna en la cual se tienen a bolas marcadas con la letra A y b bolas marcadas con la letra B . La extracción de las bolas se hará de dos posible maneras: con o sin reemplazamiento. 1) Extracción con reemplazamiento. En este caso se van extrayendo de dicha urna bolas al azar con reemplazamiento, esto es, una vez que se saca una bola esta es devuelta a la misma antes de extraer la siguiente. Se tienen dos situaciones bien diferenciadas. 1.1) Se extraen n bolas y se cuentan cuántas tienen la letra A . Sea X la variable aleatoria que mide el número de bolas A en la muestra de n bolas. En esta situación, se tiene un proceso de Bernoulli o dicotómico (pues la bola tendrá la letra A o la letra B ) repetido n veces de manera independiente (el hecho de devolver la bola a la urna hace que la probabilidad de sacar una bola A no cambie en cada extracción) y se dice que la variable aleatoria X sigue una distribución Binomial de parámetros n , el tamaño de la muestra, y p , la probabilidad de éxito, que en este caso se asociará con ser una bola A 3 . Nótese que en este caso, como se conoce el número de bolas de cada categoría, se tiene que p=a/(a+b) . 1.2) Se extraen bolas hasta obtener k de las que están marcadas con la letra A . Sea Y la variable aleatoria que mide el número de bolas B extraídas hasta que se obtienen las k marcadas con la A . Se tiene pues un experimento de Bernoulli de parámetro p repetido varias veces de manera independiente en el que se persigue, bajo la noción de éxito antes introducida, calcular el número de fracasos antes de obtener el éxito k -ésimo. Tal denición se corresponde con la llamada distribución Binomial Negativa 4 , cuya función masa de probabilidad es sencilla de obtener sin más que tener en cuenta que la bola k -ésima será un éxito y que, por tanto, se deben considerar las diferentes maneras de extraer k−1 bolas A en los y+k−1 intentos anteriores. 2) Extracción sin reemplazamiento. En este tipo de extracción las bolas se sacan al azar y de forma que una vez una bola es extraída esta no es devuelta a la misma. Bajo este supuesto se distinguen dos posibles puntos de vista. 3 El nombre de esta distribución se debe a que en la expresión de su función masa de probabilidad aparecen coecientes binomiales. Para ver más características acerca de esta distribución consúltese la Tabla 3. 4 En algunas obras, como [Balakrishnan and Nevzorov, 2004], también se le conoce bajo el nombre de distribución Geométrica de orden k o distribución de Pascal. Para conocer sus principales características véase la Tabla 4.
4 1. Un recorrido por las distribuciones notables discretas 2.1) Se extraen n bolas y se cuentan cuántas tienen la letra A . Se vuelve a denir la variable aleatoria X de la misma manera que se ha hecho con anterioridad. Sin embargo, el hecho de que ahora no se produzca reemplazamiento implica, por un lado, que ahora el tamaño de la muestra no es arbitrario sino que n∈ {0,1, ... , a +b} y, por otro lado, se produce un importante cambio en la ley de probabilidades que determina a la distribución puesto que ahora la probabilidad de extraer una bola con la letra A cambia con cada extracción. Así pues, en esta situación se dice que la variable aleatoria X sigue una distribución Hipergeométrica 5 de parámetros N , el total de la población; M , el número de éxitos en la población y, n , el tamaño de la muestra . En este caso, pues, se tiene una Hipergeométrica de parámetros N=a+b , M=a y n=n , cuya función masa de probabilidad viene dada por P{X=x}=a x b n−x a+b n,m´ax{0, n −b} ≤ x≤m´ın{a, n}. (1.1) Esta fórmula puede obtenerse de manera sencilla sin más que considerar una extracción aleatoria de n bolas en la que se tienen X=x éxitos, calcular la probabilidad asociada a tal extracción y luego multiplicar por n x , es decir, por las diferentes maneras en que pueden aparecer los x éxitos en la muestra de tamaño n . 2.2) Se extraen bolas hasta obtener k de las que están marcadas con la letra A . De nuevo se dene la variable aleatoria Y de la misma forma que se ha hecho anteriormente. No obstante, de nuevo debido a la falta de reemplazamiento, la función masa de probabilidad de Y no es la de la Binomial Negativa, ya que la probabilidad de éxito varía con cada extracción. En esta situación la variable Y se dice que sigue una distribución Hipergeométrica Negativa 6 de parámetros N , el número total de individuos en la población; M , el número de éxitos en la población y k , el número de éxitos a obtener antes de detener el proceso de extracción. Así pues, en este caso se tiene una Hipergeométrica Negativa de parámetros N=a+b, M =a y k=k y, por tanto, su función masa de probabilidad viene dada por 7 : P{Y=y}=y+k−1 k−1a+b−y−k a−k a+b a, y ∈ {0, ... , b}. (1.2) 5 El nombre se debe a que su función masa de probabilidad (véase la Tabla 5) está en progresión hipergeométrica, esto es, que P{X=k+1} P{X=k} es una función racional de k . 6 En algunas obras, como [Kotz et al., 2005], también se conoce como Hipergeométrica Inversa. 7 Véase la Tabla 6.
1.2. Bernoulli, Binomial e Hipergeométrica 5 Para comprender mejor la naturaleza de esta distribución conviene aclarar lo que representa cada elemento de la fórmula (1.2). Supóngase que se retiran la totalidad de las bolas y que se van etiquetando a medida que se extraen, siendo el 1 la primera de ellas y a+b , la última (véase la Figura 3). Considérese que Y toma el valor y , es decir que, antes de sacar k bolas que tienen la letra A , se han extraído y bolas marcadas con la letra B . Esto es equivalente a que en la posición y+k esté la k -ésima bola A . Por lo tanto, dado que las extracciones son todas equiprobables puede emplearse la Ley de Laplace 8 y, en consecuencia, en el numerador debe aparecer el número de formas de extraer bolas de tal manera que se tengan k−1 bolas A en las primeras y+k−1 posiciones, una bola A en la y+k -ésima posición 9 y las restantes a−k bolas A en las posiciones restantes ( a+b−y−k ). Finalmente, en el denominador deberá aparecer el número de formas diferentes de ordenar las a+b bolas, llegando de este modo a la expresión (1.2). Observación 1.2 . La interpretación que acaba de darse para la Hipergeométrica Negativa puede hacerse de forma similar para la Hipergeométrica sin más que hacer algunas operaciones con los números combinatorios del numerador. Se tendrían x bolas con la letra A en las primeras n posiciones y las a−x restantes estarían en las a+b−n posiciones. Esto es, que la función masa de probabilidad de la Hipergeométrica puede reescribirse como P{X=x}=n xa+b−n a−x a+b a,m´ax{0, n −b} ≤ x≤m´ın{a, n}. Observación 1.3 . Dados a y b , la cantidad de bolas de cada uno de los dos tipos que hay en la urna, para cada k∈ {1, ... , a} se tiene una distribución Hipergeométrica Negativa diferente. Considerando esa familia de distribuciones, esta verica que para todo y∈ {0,1, ... , b} se tiene que 10 a X k=1 Pk{Y=y}=a+b a−1a X k=1 y+k−1 k−1a+b−y−k a−k =a+b a−1a+b a−1 =a b+ 1, lo cual reeja cierta homogeneidad en dicha familia de distribuciones, pues para cualquier valor de Y la suma de las probabilidades en toda la familia es constante. De igual modo, la Binomial 8 En un espacio muestral formado por sucesos equiprobables la probabilidad de un suceso es igual al número de casos posibles dividido entre el número de casos probables. 9 Este término conduce al número combinatorio uno sobre uno, por lo que no es relevante. 10 La segunda de las igualdades se debe a las identidades de la Proposición I.4 y la Proposición I.5.
6 1. Un recorrido por las distribuciones notables discretas Negativa cumple también una propiedad similar: ∞ X k=1 Pk{Y=y}=∞ X k=1 y+k−1 k−1pkqy=pqy∞ X k=0 y+k kpk=pqy(1 −p)−(y+1). Observación 1.4 . El motivo por el cual las distribuciones Binomial Negativa e Hipergeométrica Negativa contienen ese adjetivo negativa en su nombre se debe a una de las identidades que se ha mencionado en la Observación 1.3, puesto que en las funciones masa de probabilidad de ambas distribuciones el factor que tienen en común se puede reescribir de manera que aparezca un coeciente binomial negativo 11 . En efecto, se cumple que y+k−1 k−1= (−1)k−1−(y+ 1) k−1. Así pues, introducidas las cuatro distribuciones, se aclararán a continuación las relaciones que aparecen reejadas en la Figura 2 (véase [Chae, 1993]): De la Hipergeométrica a la Binomial: Considérese la función masa de probabilidad de la Hipergeométrica dada por (1.1). Se tiene que, cuando a y b tienden a +∞ , dicha expresión puede aproximarse por 12 axbn−x x! (n−x)! (a+b)n n! =n! (n−x)! x! axbn−x (a+b)n=n x a a+bxb a+bn−x . En consecuencia, se tiene que, cuando a y b son sucientemente grandes la Hipergeométrica puede aproximarse por una Binomial. Intuitivamente, puede verse que, teniendo en cuenta que una de las diferencias fundamentales entre ambos modelos probabilísticos reside en que al no haber reemplazamiento la probabilidad de sacar una bola A o B cambia en cada extracción, si la cantidad de los dos tipos de bolas es muy grande, puede pensarse en que esas probabilidades son aproximadamente iguales en cada extracción y por tanto cada extracción puede pensarse como un proceso de Bernoulli con p=a/(a+b) . De la Hipergeométrica Negativa a la Binomial Negativa: En primer lugar, tal y como se desprende de las deniciones de sus respectivas funciones masa de probabilidad, ambas tienen un factor común. En cuanto a los términos en los que discrepan la idea fundamental consiste en desarrollar los números combinatorios y emplear 11 Nótese que el factorial puede generalizarse a través de la función Γ (véase la Denición I.6) a cualquier número complejo de parte real positiva. 12 Dicha aproximación se hace a partir de la fórmula de Stirling (véase la Proposición I.3), ya que permite deducir que a! (a−x)! ≃ax .
1.2. Bernoulli, Binomial e Hipergeométrica 7 de nuevo la fórmula de Stirling: a+b−y−k a−k a+b a=(a+b−k−y)! a!b! (a+b)! (a−k)! (b−y)! ≃1 (a+b)k+yakby=a a+bkb a+by . En consecuencia, cuando a y b son sucientemente grandes, la Hipergeométrica Negativa puede ser aproximada por una Binomial Negativa gracias, de nuevo, a que ahora cada extracción puede considerarse como un proceso de Bernoulli. De la Binomial a la Hipergeométrica: En la Figura 2 se dice que la Binomial construye la Binomial Negativa. Para comprender tal armación considérense dos variables aleatorias X1∈ Binomial (n1, p) y X2∈ Binomial (n2, p) independientes. Resulta que la distribución que sigue la variable aleatoria X1 condicionada a que X1+X2=k tiene la función masa de probabilidad siguiente: P{X1=x|X1+X2=k}=n1 xpxqn1−xn2 k−xpk−xqn2−k+x n1+n2 kpkqn1+n2−k =n1 x n2 k−x n1+n2 k, donde m´ax{0, k −n2} ≤ x≤m´ın{n1, k} . Con lo cual, la distribución que sigue X1 condicionada a que X1+X2=k es una Hipergeométrica (n1+n2, n1, k) . De la Binomial Negativa a la Hipergeométrica Negativa: Para comprender la relación entre ambas distribuciones en el contexto de las bolas considérense dos nuevas variables aleatorias. Sea Y1 la cantidad de bolas B sacadas antes de extraer k1 bolas A y sea Y2 el número de bolas B que se quitan hasta extraer otras k2 bolas A , donde k1, k2≥1 y, además, k1+k2≤a . Lo que se pretende es obtener la probabilidad de que Y1 tome el valor y1 condicionado a que Y1+Y2=y1+y2 , con y1, y2≥0 e y1+y2≤b . Para ello, considerando una representación similar a la de la Figura 3, se tiene que en total se extraen y1+y2+k1+k2 bolas, de las cuales k1+k2 son bolas A e y1+y2 son bolas B . Ahora bien, a diferencia de lo que ocurría anteriormente, ahora se tienen dos restricciones: La última bola siempre estará marcada con la letra A . En consecuencia, el denominador de la función masa de probabilidad buscada deberá ser el número de maneras diferentes que hay de colocar las y1+y2+k1+k2−1 bolas restantes. La bola y1+k1 -ésima debe ser A . Por lo tanto, en el numerador deberán aparecer las distintas formas que hay de extraer las k1−1 bolas A y las y1 bolas B en las y1+k1−1 primeras posiciones, así como la forma de extraer las k2−1 bolas A y las y2 bolas B en las y2+k2−1 posiciones nales.
8 1. Un recorrido por las distribuciones notables discretas Teniendo en cuenta estas observaciones se llega a que la función masa de probabilidad buscada es P{Y1=y1|Y1+Y2=y1+y2}=y1+k1−1 k1−1y2+k2−1 k2−1 y1+k1+y2+k2−1 k1+k2−1, y1, y2≥0, y1+y2≤b. Como puede observarse, tal función masa de probabilidad se corresponde con una Hipergeométrica Negativa (y1+ k1+y2+k2−1, k1+k2−1, k1) . Por último, se destacan algunos ejemplos de aplicación de estas distribuciones. En lo relativo a las aplicaciones de la Binomial Negativa 13 esta se emplea, por ejemplo, para modelar los datos obtenidos a través de tomografías por emisión de positrones ( PET ) (véase [Santarelli et al., 2017]) o en bioinformática, donde se emplea en algoritmos de búsqueda de secuencias de genoma donde se produzcan interacciones de ciertas proteínas con el ADN (véase [Spyrou et al., 2009]). En cuanto a la Hipergeométrica Negativa, se emplea en ámbitos muy diversos. Por ejemplo, en [Ridout, 1999] se emplea para modelar una situación relacionada con la memoria de una especie de pájaro que debe encontrar, de entre varios comederos cuyo contenido está oculto, el que tiene comida. Otros ejemplos de aplicación pueden encontrarse en [Cuzick, 2001], donde se propone una modicación del sistema de muestreo empleado en el análisis de los tiempos de los ensayos clínicos basado en la Hipergeométrica Negativa; en [Crooks and Brenner, 2005], donde se propone un modelo alternativo para el proceso de sustitución de aminoácidos durante la evolución de una proteína 14 y tal modelo está basado en la Hipergeométrica Negativa; o en [Autin and Gerstenschlager, 2019], donde, para mostrar a los alumnos de un curso en probabilidades la motivación de esta distribución, se presenta el juego hundir la ota 15 y se estudia la distribución que subyace a tal juego, que atendiendo a cómo se plantea, consiste en seleccionar todas las casillas donde hay un barco para conseguir hundirlos todos, llegando a que tal distribución es la Hipergeométrica Negativa. Por otro lado, en cuanto a aplicaciones de la distribución Hipergeométrica en la vida real puede destacarse la que se menciona en [Stanislevic and VoteTrustUSA, 2006], donde se explica la manera en que esta distribución se emplea para determinar el tamaño de la muestra aleatoria que 13 En la página 1 de [Al-Khasawneh, 2010] se da una formulación alternativa de la función masa de probabilidad de la Binomial Negativa donde aparece un parámetro llamado parámetro de dispersión y que, con la notación que se ha empleado a lo largo de esta sección sería igual a k . 14 Los modelos usuales están basados en procesos de Markov y tienen como inconveniente principal que supone que cada localización en la secuencia de de la proteína tiene asociada la misma distribución de aminoácidos, homogeneidad que no se corresponde con los resultados que se observan. Además, en este artículo los autores señalan de forma explícita que esta distribución es un modelo con varias aplicaciones en bioestadística que está infravalorado. 15 Para poder hablar de variable aleatoria suponen que la elección de dónde lanzar el misil se toma generando un número aleatorio, el cual está asociado con una determinada casilla del tablero.
1.3. Geométrica 9 se toma durante el proceso de auditoría al que se someten los votos americanos para detectar aquellos recintos electorales donde se haya cometido algún tipo de ilegalidad. 1.3. Geométrica Considérese una secuencia de intentos procedentes de una distribución de Bernoulli de parámetro p , con 0< p < 1 . En ese contexto, sea X la variable aleatoria que mide el número de fracasos hasta que ocurre el primer éxito. Se tiene que X sigue una distribución Geométrica y su función masa de probabilidad 16 viene dada por P{X=k}=P{“ Observar k fallos ”}·P{“ Observar un éxito en el (k+ 1) -ésimo intento ”} = (1 −p)kp, k = 0,1, ... Observación 1.5 . Tal y como puede desprenderse de su denición, se trata de un caso particular de distribución Binomial Negativa de parámetro k= 1 . Observación 1.6 . El motivo por el cual esta distribución recibe el nombre de distribución Geométrica es que su función masa de probabilidad está en progresión geométrica, siendo 1−p el factor común. La variable aleatoria X puede denirse de una manera ligeramente diferente, que será denotada como Y , viéndola como el número de intentos que son necesarios hasta obtener el primer éxito, es decir, es un tiempo de espera. Resulta evidente que se tiene la relación Y=X+ 1 y la función masa de probabilidad de la Y viene dada por PY{Y=k}=P{X+ 1 = k}=P{X=k−1}= (1 −p)k−1p, k = 1,2... Una de las propiedades más interesantes de esta distribución es la llamada propiedad de pérdida de memoria , la cual se muestra en la Proposición 1.7. Esta propiedad caracteriza a la distribución, en el sentido de que no hay ninguna otra variable aleatoria discreta 17 que la cumpla. Proposición 1.7 ([Krishnamoorthy, 2016]) . Sea X una variable aleatoria que sigue una distribución geométrica de parámetro p , con p∈(0,1) , y sean dos números naturales n y k . Se tiene que P{X≥k+n|X≥n}=P{X≥k}. 16 Para ver algunas de sus características consúltese la Tabla 7. Otras medidas notables de esta distribución pueden encontrarse en [Balakrishnan and Nevzorov, 2004]. 17 Más adelante se verá que lo mismo ocurre con la distribución Exponencial en el caso continuo.
10 1. Un recorrido por las distribuciones notables discretas Demostración. Se tiene que P{X≥k+n|X≥n}=P{X≥k+n y X≥n} P{X≥n}=P{X≥k+n} P{X≥n}=(1 −p)k+n (1 −p)n= (1 −p)k =P{X≥k}. ■ 1.4. Poisson La distribución de Poisson debe su nombre al matemático francés Siméon Denis Poisson (1781-1840), que en [Poisson, 1837] la introduce por primera vez para estudiar el número de condenas injustas que se producían en un determinado país y para ello empleaba una serie de variables aleatorias que medían el número de veces que se daba este suceso durante un intervalo de tiempo concreto. Así, esta distribución la siguen variables aleatorias que cuentan el número de éxitos en un intervalo continuo (de tiempo, espacio, volumen...) pero tal que la tasa de aparición de sucesos por unidad de tiempo, denotada por λ , permanece constante. La distribución de Poisson puede ser introducida desde diversos puntos de vista. Históricamente, esta distribución ha sido vista como un caso límite de la distribución Binomial cuando el número de intentos es grande y la probabilidad de éxito es pequeña (véase [Ramanathan, 1993]). Más concretamente, sea X una variable aleatoria siguiendo una distribución Binomial de parámetros n y p y sea λ:=np . Se tiene que P{X=x}=n xpx(1 −p)n−x=n(n−1) ···(n−x+ 1) x!λ nx1−λ nn−x =11−1 n1−2 n···1−x−1 n x!λx1−λ nn−x . (1.3) En consecuencia, tomando límites cuando n→ ∞ y p→0 en la ecuación (1.3) se llega a la distribución de Poisson 18 : P{X=x}=e−λλx x!, x = 0,1, ... (1.4) No obstante, esta distribución surge también de manera natural en el contexto de la Teoría de colas 19 a través de la resolución de ecuaciones diferenciales (véanse [Ramanathan, 1993] o [Bowers, 2008]). Más concretamente, hay unos procesos que se pueden enmarcar dentro de la 18 Para conocer más características de esta distribución consúltese la Tabla 8. 19 En concreto, está relacionada con los llamados procesos de colas de espera , que, por ejemplo, en el caso de los supermercados, consiste en la llegada de un nuevo cliente a la caja de pago. Estos procesos modelan situaciones de esa naturaleza tratando de predecir tanto la longitud de la cola como el tiempo de espera de los clientes. Para profundizar más en este campo puede consultarse [Abad, 2002].
1.4. Poisson 11 Teoría de colas conocidos como procesos de Poisson . Estos son un tipo de series de tiempo 20 que se caracterizan por modelar una situación donde se conoce el tiempo medio que hay entre dos sucesos pero no ocurre así con el momento de aparición de cada suceso, que es donde reside la aleatoriedad. Estos procesos poseen tres características fundamentales. Por un lado, los tiempos que hay entre las apariciones de los sucesos son independientes, es decir, el tiempo que hay entre dos eventos consecutivos tiene ausencia de memoria, en el sentido de que el tiempo que tarda en aparecer un suceso no se ve afectado por lo que haya tardado antes. Por otro lado, la cantidad media de sucesos por unidad de tiempo es constante y, por último, dos sucesos no pueden ocurrir a la vez, es decir, dado un subintervalo de tiempo muy pequeño o no ocurre ningún suceso u ocurre exactamente uno 21 . Ejemplo 1.8. Supóngase que un pequeño supermercado es conocedor, en base a un estudio realizado, de que en su establecimiento pasan por la caja de pagos en media unas 50 personas al día. El hecho de que en un determinado momento aparezca un cliente no afecta a la probabilidad del tiempo que tarde en aparecer el siguiente. Por otro lado, está claro que en un determinado lapso de tiempo muy pequeño pueden ocurrir solamente dos cosas: o aparece un cliente para pagar o no aparece ninguno. Observación 1.9 . Tal y como puede verse en [Tse et al., 2014] los procesos de Poisson surgen también en otros ámbitos muy diferentes al del ejemplo del supermercado como la emisión de partículas radioactivas, los archivos que llegan a un servidor, la aparición de desastres naturales... La utilidad de los procesos de Poisson radica en poder resolver cuestiones como hallar la probabilidad de que haya un cierto número de éxitos en un periodo de tiempo determinado. La respuesta a esto se encuentra, precisamente, en la distribución de Poisson, ya que esta permite calcular la probabilidad de observar x eventos en un determinado periodo de tiempo a partir de la longitud de ese intervalo de tiempo ( t ) y la media de aparición de los sucesos por unidad de tiempo ( α ). En efecto, considérese un proceso de Poisson y sea f(x, t) la probabilidad de que ocurran x sucesos en un intervalo de tiempo de longitud t . Lo que se pretende es encontrar una fórmula explícita para tal función de probabilidad y para ello se busca calcular la probabilidad de que ocurran x sucesos en un intervalo ligeramente mayor que [0, t] , esto es, en un intervalo [0, t+∆t] , lo cual viene representado por f(x, t + ∆t) . Dado que se está en un proceso de Poisson, este verica las tres condiciones mencionadas anteriormente, por lo que se tiene lo siguiente: 20 Una serie de tiempo es un conjunto de datos tomados a lo largo de períodos de tiempo consecutivos, normalmente de igual longitud (véase [Encyclopedia Britannica, 2022]). 21 Por lo tanto, la aparición de un suceso en un intervalo de tiempo innitamente pequeño puede verse como un proceso de Bernoulli, lo cual enlaza con lo visto al comienzo de esta sección cuando se establecía la relación entre la distribución de Poisson y la Binomial: la distribución de Poisson permite contar sucesos raros en un intervalo continuo.
18 1. Un recorrido por las distribuciones notables discretas distribución de Pareto. Para describir esto en términos matemáticos, sea X la variable aleatoria que mide la renta per cápita de un individuo y sea X(r) la que mide la del individuo r -ésimo. La Ley de Zipf establece que X(r)=C r−α=x para una cierta constante C , luego la probabilidad de que una renta per cápita sea mayor o igual que x es proporcional a r , es decir: P{X≥x}=C′r=C′x C−1 α. Esta última expresión es precisamente la representación de la situación desde el punto de vista de la distribución de Pareto y en base a ella se puede calcular la probabilidad de que una determinada renta per cápita sea igual a x . En efecto, se tiene que P{X < x}= 1 −P{X≥x}= 1 −C′x C−1 α. Por lo tanto, derivando con respecto a x se llega a que P{X=x}=C′C1 α αx−(1 α+1). (1.10) En consecuencia, se acaba de ver la manera en que una distribución de Zipf de exponente α se corresponde con una función de distribución de Pareto de parámetro 1/α . Además, la función de densidad (1.10) permite ver que se trata de una distribución de potencias de exponente (1/α + 1) . Observación 1.17 . Sea X una variable aleatoria que sigue una distribución de Zipf de parámetros n , el número de elementos, y α . Entonces se tiene que P{X=k}=C k−α, k = 1, ... , n. Para calcular el valor de esa constante C > 0 basta con aplicar que la suma de las probabilidades para todo k∈ {1, ..., n} debe ser igual a 1 . Así, se cumple que n X k=1 C1 kα = 1 ⇐⇒ C="n X k=1 1 kα#−1 . De este modo se llega a la función masa de probabilidad que aparece en la Tabla 10. Observación 1.18 . La expresión que se ha calculado para la C en la Observación 1.17 guarda una evidente relación con la función Zeta de Riemann, ζ(s) , que es tal que ζ(s) = 1 + 1 2s +1 3s +···+1 ks +··· De este modo, en algunas obras como [Ross, 2019] se reeren a la ley o distribución de Zipf como la distribución Zeta .
1.6. Zipf 19 En cuanto a la naturaleza de los datos que verican esta ley, se trata, en general, de conjuntos de datos ordenados en los cuales sucesos grandes son poco frecuentes, mientras que sucesos pequeños ocurren con mucha más frecuencia 28 . Por ejemplo, hay pocas personas que sean millonarias pero hay muchas que tiene una renta humilde; hay pocas ciudades de más de un millón de personas, mientras que hay muchas que tienen miles... En estas situaciones, por una cuestión de escala, los datos suelen representarse en grácas de tipo log-log , lo que además ayuda a entender lo que trata de explicar la Ley de Zipf. Para claricar lo anterior, tómense logaritmos en la expresión (1.9) y considérese la variable aleatoria X(r) . Se tiene que Elog(X(r))= log(C)−αlog(r) = β0+β1log(r). Con lo cual, la ley de Zipf en una gráca de tipo log-log lo que pretende es ajustar un modelo de regresión en media de log(X(r)) sobre log(r) y la pendiente asociada a dicho modelo está cercana a −1 . Por consiguiente, para dar una interpretación del parámetro α hay que irse al nivel de log(r) , ya que representa la tasa de variación de la variable log(X(r)) con respecto al logaritmo de la posición que ocupa el dato. Dicho de otra manera, el parámetro α establece la rapidez con la que va disminuyendo en escala logarítmica la probabilidad asociada a la posición r -ésima. Por último, en cuanto a las aplicaciones de esta distribución, esta se emplea en una gran diversidad de campos (véase [Li, 2002]) como son estudiar la frecuencia de uso de las palabras en un idioma, cuanticar el prestigio de un investigador ordenando los investigadores según el número de citas que tienen sus artículos, establecer un ranking de tamaño de empresas en base a su número de empleados o a sus benecios anuales, estudiar la intensidad de las erupciones solares, estudiar las visitas a portales web o incluso se emplea para describir las interacciones entre pacientes en un centro psiquiátrico (véase [Piqueira et al., 1999]). Sin embargo, algo que conviene destacar también es que no cualquier conjunto de datos verica necesariamente esta ley, ya que por ejemplo, cuando se tienen pocos datos resulta complicado pensar en que una distribución de potencias como es la ley de Zipf modele la situación adecuadamente. Así, hay una gran cantidad de ejemplos en los cuales en una escala log-log los datos ordenados se desvían signicativamente de la teórica recta de regresión que deberían seguir bajo la ley de Zipf. 28 Grácamente esto se corresponde con nubes de puntos que presentan una cola en el extremo positivo del eje de abscisas.
Capítulo 2 Un recorrido por las distribuciones notables continuas Siempre me recuerdo a mí mismo que lo que se observa es como mucho una combinación de probabilidades y resultados, no sólo resultados. Nassim Nicholas Taleb A lo largo de este capítulo se estudiarán algunas de las distribuciones continuas más importantes. Cabe destacar que en el caso de este tipo de variables aleatorias sus distribuciones de probabilidad quedan perfectamente determinadas si se determina su función de densidad. Denición 2.1 (Función de densidad) . Sea (Ω,A,P) un espacio de probabilidad y sea X: Ω −→ R una variable aleatoria continua. La función de densidad de X es una función fX:R−→ [0,1] tal que fX(x) = l´ım h→0 P{x−h≤X≤x+h} 2h, x ∈R, donde P es una probabilidad con la axiomática de Kolmogorov y, además, verica que Z+∞ −∞ f(x)dx = 1. Además, se denomina soporte de X al conjunto de los x∈R tales que fX(x)>0 . 21
22 2. Un recorrido por las distribuciones notables continuas 2.1. Weibull La distribución de Weibull debe su nombre al ingeniero y matemático sueco Ernst Hjalmar Waloddi Weibull 1 (1887-1979) que la introdujo en su trabajo [Weibull, 1951] para comprobar que esta ajustaba mejor que muchas otras distribuciones conocidas una serie de datos procedentes de varios ámbitos diferentes (ingeniería, agricultura, biometría ...). En él, parte de la idea de que cualquier función de distribución F(x) puede ser expresada de la forma F(x) = 1 −e−φ(x), (2.1) para una determinada función φ(x) , a la que luego se le exige que sea positiva, no decreciente y que se anule en un determinado xu . Bajo esas condiciones Weibull establece que la función más simple es 2 φ(x) = (x−xu)m x0 . Sustituyendo en (2.1) se tiene que F(x)=1−e−(x−xu)m x0. Para dilucidar la motivación de esta distribución, Weibull introduce el ejemplo de una cadena formada por n eslabones, para cada uno de los cuales se ha comprobado que la probabilidad de rotura para una carga de como mucho x es P 3 . Además, se entenderá que la cadena está rota cuando lo esté alguno de sus eslabones, es decir, que será tan fuerte como el más débil de sus eslabones. En este escenario, se desea calcular la probabilidad de que la cadena no se rompa (en adelante, 1−Pn ), la cual, según las hipótesis que se han considerado, será igual a la probabilidad de que no se rompa ninguno de los eslabones, es decir: 1−Pn= (1 −P)n. (2.2) Ahora bien, si la función de distribución de cada uno de los eslabones es de la forma (2.1), se obtiene una expresión mucho más sencilla: (1 −P)n=he−φ(x)in=e−n φ(x). 1 Si bien fue él quien popularizó su uso, anteriormente varios matemáticos como Fréchet [Fréchet, 1927], von Bortkiewicz [von Bortkiewicz, 1922], von Mises [von Mises, 1923] o Fisher y Tippet [Fisher and Tippett, 1928] ya habían publicado artículos donde se trataban aspectos de los valores extremos de una muestra que de alguna manera sirvieron como primer paso hacia esta distribución y poco tiempo después comenzaron a estudiarse algunas de sus aplicaciones (véase [Rosin et al., 1933]). 2 Esto seguramente se trate de una errata en el artículo, ya que los datos que aparecen en él se corresponden con φ(x) = x−xu x0m . Cuestiones como estas son tratadas con mayor detalle en [Tsu et al., 1951]. 3 Por respeto al autor y a su obra se ha decidido preservar la notación que emplea él pero verdaderamente tanto P como Pn son valores que, en principio, dependen tanto de la x como de cada uno de los eslabones.
2.1. Weibull 23 Con lo cual, considerando la ecuación (2.2) se llega a una función de distribución que es, en una forma muy general, la de la distribución de Weibull 4 : Pn= 1 −e−n φ(x). Uno de los principales ámbitos de aplicación de esta distribución es el análisis de supervivencia y los estudios de abilidad. En estos campos cobran mucha importancia dos funciones vinculadas a la función de distribución y que caracterizan a las distribuciones: la función de supervivencia y la tasa o razón de fallo. 1) Función de supervivencia : la función de supervivencia de una variable aleatoria, denotada por s(x) , se dene de la manera siguiente: s(x) = 1 −F(x). En cuanto a su interpretación, representa la probabilidad de supervivencia del individuo (o lo análogo en el contexto en el que se esté) pasado el tiempo x . 2) Razón de fallo : la razón de fallo, también llamada función de riesgo, de una variable aleatoria, se denota por τ(x) y se dene como sigue: τ(x) = f(x) s(x). En lo relativo a su signicado, establece la probabilidad de que un superviviente de x años fallezca inmediatamente tras alcanzar dicha edad, es decir, es una medida del riesgo de fallecimiento que tiene asociado el individuo que ha sobrevivido hasta el instante x . Un sencillo cálculo permite llegar a que la razón de fallo de la Weibull de parámetros α y λ viene dada por τ(x) = λα(λ x)α−1, x > 0. Es decir, que atendiendo a este nuevo enfoque, la distribución de Weibull trata de modelar la distribución de los fallos cuando la función de riesgo es, salvo constantes, una potencia del tiempo. En la Figura 4 puede observarse la forma de la razón de fallo para diferentes valores del parámetro forma α : α < 1 : la razón de fallo decrece con el tiempo. Este fenómeno acostumbra a darse cuando hay una alta mortalidad infantil , que se interpreta como un fallo prematuro debido a alguna tara o error de diseño, y luego la razón de fallo disminuye con el tiempo debido a que se detecta el fallo y se subsana el error. 4 Para ver algunas de las características de esta distribución véase la Tabla 12 .
24 2. Un recorrido por las distribuciones notables continuas α= 1 : la razón de fallo es constante y representa la vida útil del individuo. α > 1 : la razón de fallo aumenta con el tiempo y ocurre cuando se ha producido un desgaste en el individuo. Hay algunas situaciones donde se suponen ciertas consideraciones físicas acerca del mecanismo de fallo y estas conducen a una determinada distribución pero, con mucha mayor frecuencia, a la hora de seleccionar la distribución prima el criterio de selección basado en la bondad de ajuste de la distribución a los tiempos de fallo observados. En ese contexto, en el ámbito de los estudios de abilidad o en el biológico hay multitud de ocasiones, como podría ser el caso de un ser humano, en las cuales los datos llevan asociada una función de riesgo que no es monótona, sino que frecuentemente tiene forma de bañera (véase la Figura 4), por lo que la distribución de Weibull no ajusta adecuadamente los datos. Por esa razón, a lo largo de los años se han ido estudiando distribuciones de probabilidad cuya función de riesgo fuese más exible y permitiese modelar situaciones donde tuviese una forma similar a una bañera y que, además, extendiese a la distribución de Weibull. Algunas de tales distribuciones se presentan a continuación. 2.1.1. Exponencial La distribución Exponencial se trata de un caso particular de distribución de Weibull de parámetro α= 1 , tal y como puede comprobarse a partir de la Tabla 13. Esta distribución se emplea en situaciones donde los datos tienen una razón de fallo aproximadamente constante ya que esta cumple que τ(t) = λ e−λ t 1−(1 −e−λ t)=λ, t > 0. No obstante, esta característica hace que existan numerosos ámbitos donde este modelo no sería aplicable como, por ejemplo, para modelar el tiempo de vida de una persona, pues se estarían considerando individuos que no envejecen, en tanto que la probabilidad de fallecimiento alcanzada cualquier edad es siempre la misma, λ . Por último, cabe destacar que esta distribución se dene como el tiempo transcurrido entre dos sucesos consecutivos en un proceso de Poisson de parámetro λ y es, además, la generalización de la distribución Geométrica al caso continuo. De este modo, de una manera análoga a como se procedió en la demostración de la Proposición 1.7, se puede demostrar que la Exponencial posee la propiedad de ausencia de memoria. Además, es la única distribución continua que verica esta propiedad. Proposición 2.2 ([Ross, 2019]) . Sea X una variable aleatoria continua. Entonces X verica la propiedad de ausencia de memoria si, y solo si, sigue una distribución exponencial.
2.1. Weibull 25 Demostración. Supóngase que X es una variable aleatoria no negativa que verica la propiedad de ausencia de memoria, esto es: P{X > s +t, X > t} P{X > t}=P{X > s} ⇐⇒ P{X > s +t}=P{X > s}P{X > t}. (2.3) Si se escribe F(x) = P{X > x} entonces la ecuación (2.3) implica que F satisface la ecuación diferencial f(x+t) = f(x)f(t) . Esta ecuación diferencial tiene como única solución continua por la derecha la función f(x) = e−λx , por lo que F(x) = 1 −F(x) = 1 −e−λx . En efecto, tal solución única se puede obtener de la manera siguiente. Sea f(x) solución de la ecuación diferencial. Entonces se cumple que fm n=f1 n+1 n+ m ···+1 n=f1 nf1 nm ···f1 n=f1 nm . (2.4) De igual modo también se verica lo siguiente: f(1) = f1 n+1 n+ n ···+1 n=f1 nn =⇒f1 n= [f(1)] 1 n. (2.5) En consecuencia, juntando las ecuaciones (2.4) y (2.5) se tiene que fm n= [f(1)]m n. Por lo tanto, como se busca una solución continua por la derecha se tiene que la solución f(x) será tal que f(x) = [f(1)]x y como f(1) = [f(1/2)]2≥0 , se tiene que f(x) = e−λx, λ =−ln (f(1)) . En cuanto a la otra implicación, resulta inmediato comprobar que la función de distribución de la Exponencial verica tal ecuación diferencial. ■ 2.1.2. Weibull-Geométrica En [Barreto-Souza et al., 2011] se introduce una nueva familia de distribuciones de probabilidad conocida como Weibull-Geométrica . Dicha familia surge de manera natural a partir de la Exponencial-Geométrica 5 , ya que, como se ha visto con anterioridad, la Weibull generaliza la distribución Exponencial. La utilidad principal de esta familia es que, al contrario de lo que ocurría con la Weibull, permite ajustar datos cuya razón de fallo es no necesariamente monótona. 5 Por una cuestión de brevedad, tanto esta distribución como la Exponencial-Geométrica extendida no son tratadas en este trabajo pero pueden ser consultadas en [Adamidis and Loukas, 1998] y [Adamidis et al., 2005], respectivamente.
26 2. Un recorrido por las distribuciones notables continuas El contexto en el que surge esta familia de distribuciones es en el de los llamados dentro del análisis de supervivencia problemas de riesgos competitivos . Sean X1, X2, ... , Xk variables aleatorias independientes cuyas funciones de distribución vienen dadas por F1(x), ... , Fk(x) , respectivamente, y donde k es el valor de una variable aleatoria discreta Z que toma valores en N . Sea, además U= m´ın{X1, ... , Xn} . Se denomina problema de riesgos competitivos al problema de obtener las distribuciones de las variables Xi, i = 1, ... , n , cuyos valores no son observables, a partir de la distribución de U , cuyos valores sí lo son. Intuitivamente, estos problemas se pueden interpretar como un individuo que se encuentra sometido a k posibles causas de muerte o fallo cuyos tiempos de vida (en el sentido de tiempo hasta ser observado) vienen determinados por las variables Xi antes introducidas. En esa circunstancia está claro que el único valor que podrá ser observado será el mínimo de los tiempos de vida. Tal y como se muestra en [Basu and Ghosh, 1980], este tipo de problemas son muy comunes y surgen en ámbitos muy diversos, que pasan desde un sistema eléctrico cuyas componentes están conectadas en serie hasta un humano que está sometido a una serie de factores que pueden causarle la muerte. En ese mismo artículo, se lleva a cabo un análisis acerca de cómo abordar este tipo de problemas cuando las variables no son independientes, cuando sí lo son pero no están idénticamente distribuidas y algún caso más que involucra el concepto de identicabilidad y que no se tratará en este trabajo (véase [Basu and Klein, 1982]). Sin embargo, en el caso de la Weibull-Geométrica el caso que se considera es el más sencillo, puesto que es el caso de variables aleatorias independientes e idénticamente distribuidas. Más concretamente, se tiene que siguen una Weibull (λ, α) y la variable Z , cuyo valor determina el número de fallos, sigue una distribución Geométrica de parámetro p 6 . En cuanto al cálculo de la densidad marginal de U , previamente es necesario realizar algunos cálculos: fU(x|z) = z z Q i=1 1−FXi(x)fXi(x) = z α λ (λ x)α−1e−(λ x)αz Q i=1 1−1 + e−(λ x)α =z α λαxα−1e−(λ x)α(1+z), fU(x) = ∞ P z=1 f(x, z) = ∞ P z=1 fX(x|z)·fZ(z) = ∞ P z=1 z α λαxα−1e−(λ x)α(1+z)(1 −p)pz−1 =α λα(1 −p)xα−1∞ P z=1 z pz−1e−(λ x)α(1+z) =α λα(1 −p)xα−1e−2(λ x)α1−p e−(λ x)α−2, (2.6) 6 En este caso la geométrica se ve como el número de componentes defectuosas (o lo análogo en el problema concreto) que son necesarias hasta obtener el primer fracaso (no defecto), donde p es la probabilidad de que una componente sea defectuosa, es decir: fZ(z) = (1 −p)pz−1, z ∈N . No obstante, aunque en el artículo no se haga mención, esta visión es algo limitada porque supone que una vez detectada una pieza que no sea defectuosa, ya no habrá ninguna más que lo sea.
2.1. Weibull 27 donde el último de los pasos se debe a que ∞ P n=1 n kn−1= (1−k)−2 . Mediante unos sencillos cálculos, a partir de esta función de densidad puede obtenerse la función de riesgo de esta distribución: τ(x) = α λαxα−1 1−p e−(λ x)α. (2.7) Observación 2.3 . La función de densidad (2.6) lo es para todo p≤1 , aunque en caso de tomar p≤0 no podría interpretarse como el parámetro de una distribución Geométrica. No obstante, si se toma p= 0 la función de densidad que se obtiene es la de una Weibull (λ, α) . Además, tal y como puede observarse en la Figura 5, cuando este parámetro tiende a 1 , la distribución tiende a degenerar en el cero, por lo que p puede verse como un parámetro de escala. Otras dos distribuciones que se pueden obtener a partir de la Weibull-Geométrica son, por un lado, la Exponencial-Geométrica, que se obtiene cuando α= 1 y p∈(0,1) , y, por otro, la Exponencial- Geométrica extendida, que se obtiene cuando α= 1 y p < 1 . Lo interesante de esta distribución es su función de riesgo y la amplia gama de posibilidades que ofrece en función de los tres parámetros, tal y como muestra la Figura 4. Por ejemplo, un sencillo análisis de la fórmula (2.7) permite ver que si p∈[−1,1) y α∈(0,1] esta es decreciente e incluso, como se puede apreciar en la Figura 6, hay casos donde esta no es monótona. En denitiva, se ha obtenido una familia de distribuciones cuya función de riesgo es su- cientemente exible como para ajustar datos cuya razón de fallo sea más general que la de una Weibull. En la práctica lo que se hará es estimar (p, λ, α) , que en [Barreto-Souza et al., 2011] se hace por el método de máxima verosimilitud empleando el algoritmo de esperanza-maximización de tal manera que la distribución se ajuste lo mejor posible a los datos. 2.1.3. Weibull-Geométrica complementaria En [Tojeiro et al., 2014] se introduce una nueva familia de distribuciones de tiempos de vida triparamétrica llamada Weibull-Geométrica complementaria . Dicha familia se caracteriza por tener una función de riesgo unimodal, creciente y decreciente y ser una distribución complementaria a la introducida en [Barreto-Souza et al., 2011], en el sentido de que la contiene como un caso particular. El contexto en el que surge esta familia de distribuciones es dentro de los problemas de riesgos complementarios , que son un tipo problemas que surgen en el mismo escenario que los de riesgos competitivos con la salvedad de que ahora los únicos valores que son observables son los de V= m´ax{X1, ... , Xn} . De nuevo en [Basu and Ghosh, 1980] pueden encontrarse multitud de
34 2. Un recorrido por las distribuciones notables continuas Observación 2.17 . Además de haber sido desarrollados en la misma época, a la vista de las dos propiedades que se han presentado en esta sección queda patente la estrecha relación que existe entre la teoría de valores extremos y el Teorema Central del Límite. Este vínculo se pone de maniesto con claridad a través de un sencillo ejemplo en la página 3 de [De Haan et al., 2006]: El neumático de un coche puede fallar de dos maneras diferentes. Cada día de conducción se va desgastando un poco el neumático y, tras mucho tiempo, el deterioro acumulado acaba provocando la rotura del neumático (es decir, las sumas parciales de los deterioros superan un determinado umbral). Pero también durante la conducción uno podría tropezar con un bache o chocar contra la acera. Accidentes como estos pueden o bien no tener efecto alguno sobre el neumático o bien provocar un pinchazo. En el segundo caso, es tan solo un gran valor observado el que provoca el fallo, lo cual signica que el máximo de los deterioros superó algún umbral . No obstante, hay una diferencia bastante notoria entre el Teorema 2.9 y el Teorema 2.12 y es que, así como el Teorema Central del Límite establece un criterio de convergencia para la distribución de la media muestral, el Teorema de Fisher-Tippett-Gnedenko tan solo establece que en caso de que la distribución del máximo normalizado converja a una distribución no degenerada entonces lo hará a uno de esas tres clases de distribuciones pero no establece que esta deba converger. Para nalizar esta sección se van a comentar algunas de las numerosas aplicaciones que tienen las distribuciones de Valores Extremos. Por ejemplo, en [Koutsoyiannis, 2003] se puede ver lo extendido que está el uso de esta familia de distribuciones en el ámbito de la hidrología. En la página 306 de la citada referencia se argumenta desde cuatro puntos de vista diferentes el por qué es la distribución de Gumbel la que más se utiliza a la hora de estudiar la pluviosidad máxima, si bien más adelante se pone en tela de juicio su uso y se propone como alternativa la distribución de Fréchet. Cabe destacar que obtener un buen modelo de las lluvias extremas es algo de vital importancia en el mundo de la ingeniería, ya que las pruebas de inundaciones a las que son sometidas las estructuras y la maquinaria que se desarrollan están diseñadas a partir de modelos de tormenta cticios, por lo que conseguir un modelo que se adapte lo mejor posible a las tormentas reales es estratégico para que puedan soportar esas condiciones extremas una vez que empiecen a ser utilizadas 17 . Persiguiendo ese objetivo, en [Nadarajah, 2006], inspirados en la generalización de la distribución Exponencial (llamada Exponencial Exponenciada ) que se propone en [Gupta et al., 1998] y [Gupta and Kundu, 2001], se introduce la Gumbel Exponenciada y se muestra que este puede ser un modelo mejor que la Gumbel a través de una serie de datos pluviométricos correspondientes a las precipitaciones diarias máximas anuales de Orlando, Florida, entre los años 1901 y 2001. No obstante, más allá del ámbito climatológico, las distribuciones de Valores Extremos tienen muchas aplicaciones tal y como puede observarse en la página 179 de [Kotz and Nadarajah, 2000], donde aparecen recogidos más de 50 campos 17 Para profundizar en los métodos que se emplean para modelar las lluvias máximas pueden consultarse obras como [Katz et al., 2002], [Sutclie, 1978] o [Yevjevich, 1972].
2.3. Normal 35 en los que estas distribuciones se emplean como modelo: estudio de la resistencia de sólidos a rotura por fatiga (véanse [Weibull, 1939a] y [Weibull, 1939b]), modelización direccional de las velocidades del viento (véase [Coles and Walshaw, 1994]), modelado de reuniones de multitudes de personas en el estudio de problemas de cargas vivas extremas (véase [Kanda, 1994]), diseño de redes (véase [Greis and Wood, 1981])... Por último, cabe mencionar que en los últimos años se han empleado también dentro del ámbito nanciero en el estudio de los riesgos de cola 18 en algunos de los principales índices bursátiles (véase [Gilli et al., 2006]), así como para asignar activos bajo restricciones de tipo safety-rst , calcular pérdidas esperadas o estudiar dependencias entre mercados en condiciones de estrés (véase [Rocco, 2014]). 2.3. Normal La distribución Normal o distribución Gaussiana está estrechamente relacionada con la gura de Carl Friedrich Gauss (1777-1855) debido a que la introduce en su obra [Gauss, 1823]. Sin embargo, cabe mencionar que esta distribución ya aparece, aunque no de manera directa, en 1738 en una obra del matemático francés Abraham de Moivre (1667-1754) (véase [De Moivre, 1738]), en la cual al estudiar los coecientes binomiales de (a+b)n acaba por llegar a una expresión que algunos matemáticos (véase [Pearson, 1926]) consideran la primera aparición de la Normal. Ahora bien, por aquel entonces, tal y como se arma en la página 76 de [Stigler, 1986], de Moivre no conocía el concepto de función de densidad pero, aún así, veía esa expresión como una curva: Si se consideran los términos del binomio de manera que dada una línea recta estos estén colocados en posición vertical, igualmente espaciados, y formando un ángulo recto por encima de la línea recta, los extremos de los términos siguen una curva. La curva así descrita tiene dos puntos de inexión, uno a cada lado del término máximo 19 . La situación que de Moivre describía era como la representada en la Figura 9 y esa intuición le llevó a desarrollar el siguiente teorema, cuya demostración puede ser encontrada en las páginas 156 y 157 de [Papoulis and Pillai, 2002]. Teorema 2.18 (de Moivre-Laplace, [Papoulis and Pillai, 2002] p. 105) . Sea p∈(0,1) jo y sea q= 1 −p . Entonces, para k∈(np −√npq, np +√npq) se tiene que l´ım n→∞ n kpkqn−k 1 √2πnpq e−(k−np)2 2npq = 1. (2.11) 18 Es un tipo de riesgo de cartera que surge en situaciones que tienen una pequeña probabilidad de ocurrir. Concretamente, los riesgos de cola aparecen cuando la probabilidad de que una inversión se aleje más de tres desviaciones típicas de la media es mayor que lo que ocurre en el caso de la distribución Normal (véase [Hayes, 2022]). 19 Para una introducción histórica más detallada de la distribución Normal puede consultarse [Johnson et al., 1994] de la página 85 a la 88.
36 2. Un recorrido por las distribuciones notables continuas Observación 2.19 . El Teorema de de Moivre-Laplace no es más que un caso particular del Teorema 2.9, ya que el numerador y el denominador de 2.11 se corresponden con la función masa de probabilidad de una Binomial (n, p) (es decir, una suma de n distribuciones de Bernoulli de parámetro p independientes) y la función de densidad de una Normal (np, npq) , respectivamente 20 . Ahora bien, como se está aproximando una variable discreta por una continua es necesario aplicar una corrección por continuidad , esto es, P{X discreta =k} ≃ P{k−0.5< X continua < k + 0.5}. Observación 2.20 . Sea X∈ Poisson (λ) . Entonces 21 Xd =X1+X2+···+Xλ , con Xi∈ Poisson (1) , i= 1, ... , λ . Por lo tanto, aplicando el Teorema 2.9 se tiene que n P i=1 Xi−λ √λ d −−−−−−−→ λ→∞ Normal (0,1). Es decir, que una Poisson (λ) puede aproximarse (de nuevo aplicando la corrección por continuidad) mediante una Normal (λ, λ) para valores de λ grandes. Un ejemplo de esta aproximación puede verse en la Figura 10. Así pues, el Teorema Central del Límite es uno de los principales motivos por los cuales el uso de la distribución Normal está tan extendido y por los que tiene una importancia tan grande dentro de la Estadística, ya que permite encontrar la distribución asintótica de la suma de variables aleatorias independientes 22 bajo condiciones muy generales como son la existencia de media y varianza. No obstante, no todas las aplicaciones de la distribución Normal son directamente a través de este teorema. Por ejemplo, en multitud de ocasiones cuando no se tiene ningún tipo de información acerca de la distribución que siguen unas determinadas observaciones se suele suponer que es una Normal 23 . El motivo por el cual se acostumbra hacer esto está vinculado con la llamada cota de Fréchet-Cramer-Rao 24 , que en la estimación de un parámetro determinista es una cota inferior de la varianza de un estimador insesgado. Un análisis detallado acerca de la inuencia que tiene la suposición de Normalidad en dicha cota puede ser encontrado en los 20 Tanto la función de densidad como otras características de la distribución Normal pueden ser consultadas en la Tabla 16. Cabe destacar que su función de distribución no tiene primitivas elementales pero puede calcularse mediante métodos numéricos y obtener valores concretos que luego son recogidos en tablas (véase [Abramowitz and Stegun, 1972] p. 966). 21 A partir de la función característica resulta sencillo demostrar que si X1∈ Poisson (λ1) y X2∈ Poisson (λ2) entonces X1+X2∈ Poisson (λ1+λ2) . 22 Hay versiones del Teorema Central del Límite en las cuales no es necesario suponer la independencia de las variables aleatorias. Para profundizar en ello y en otras generalizaciones de dicho teorema puede consultarse el Capítulo 1 de [Prokhorov and Statulevi£ius, 2000]. 23 Esta suposición puede ser contrastada mediante test estadísticos como el Test de Shapiro-Wilk (véase [Shapiro and Wilk, 1965]). 24 Para profundizar más acerca de este concepto puede leerse la Sección 6.3 de [Vélez and García, 2012].
2.3. Normal 37 artículos [Park et al., 2013] y [Stoica and Babu, 2011]. En la época en que fue introducida esta distribución, Gauss dedujo que en el modelo lineal que había planteado la distribución que debían seguir los errores era una Normal (véase [Seal, 1967], páginas de la 3 a la 6). Años más tarde, el físico y matemático escocés James Clerk Maxwell (1831-1879) en uno de sus trabajos sobre teoría cinética de gases se encontró con la distribución Normal como la distribución de las componentes ortogonales de la velocidad de partículas que se mueven libremente en el vacío (véase [Maxwell, 1860], páginas 22 y 23) llegando a decir: Las velocidades están distribuidas entre las partículas de acuerdo a la misma ley con la que están distribuidos los errores entre las observaciones en la teoría del método de mínimos cuadrados . De este modo, se convertía en el primer paso hacia describir el movimiento de los gases a través de una función estadística en lugar de una determinista. Trabajos como los que se han citado sirvieron de inspiración a numerosos cientícos a lo largo de los años y con el tiempo el uso de la distribución Normal fue extendiéndose cada vez a ámbitos más diversos. Así, en nuestros días esta distribución está presente en ámbitos como la biología y la medicina, donde, por ejemplo, la distribución de la glucosa en los perros (véase [Kaneko et al., 2008], p. 5), las dimensiones de la tiroides de los recién nacidos medidas mediante técnicas de sonografía (véase [Degroot et al., 2016], p. 1412) o algunos modelos de codicación de la información en la teoría de detección de señales 25 (véase [Gabbiani and Cox, 2010], p.) siguen una distribución Normal; en la psicología, donde la puntuación que recibe una persona como co- eciente intelectual en la mayor parte de los test de inteligencia sigue una distribución Normal 26 ; en la química, donde en el campo de la espectroscopía electrónica se emplea la distribución Normal para estudiar la absorción continua de rayos ultravioleta de determinadas sustancias (véase [Salman, 1999], p. 252)... No obstante, en contra de lo que pudiera parecer, el universo de las distribuciones de probabilidad no se reduce a la distribución Normal, ya que hay variables aleatorias cuya distribución puede presentar, por ejemplo, algún tipo de asimetría que haga que esté lejos de ser Normal, como ocurre con la distribución del colesterol y los triglicéridos en los humanos (véase [Dasgupta and Wahed, 2014], p. 50). En ocasiones, mediante algún tipo de transformación (lo- 25 Es una teoría matemática que se emplea para detectar de manera óptima las señales que se encuentran embebidas en ruido y que en el campo de la neurociencia se emplea para medir la información acerca de estímulos sensoriales transmitida por las neuronas o conjuntos de las mismas. 26 En un test WAIS ( Wechsler Adult Intelligence Scale ) la variable aleatoria que mide la puntuación de la persona que hace el test sigue una distribución Normal de media 100 y desviación típica 15 (véase [Moore et al., 2015], página 95). Así, una persona con un IQ de 120 estará en aproximadamente el 9% de personas con mayores habilidades intelectuales.
38 2. Un recorrido por las distribuciones notables continuas garitmos, transformaciones Box-Cox ...) es posible hacer que los datos adquieran una distribución que se asemeje más a la Normal pero hay veces en las que tales transformaciones no son efectivas. Como puede verse en la Figura 12 hay una gran variedad de distribuciones y es fundamental conocer un buen número de ellas para poder encontrar aquella que mejor se adapte a las observaciones. En esa misma gura pueden observarse bastantes distribuciones que están estrechamente ligadas a la distribución Normal, algunas de las cuales serán presentadas a continuación. 2.3.1. Ji-cuadrado La distribución χ2 ( ji-cuadrado ) de k grados de libertad surge de manera natural en el contexto de las distancias cuadráticas, ya que es la distribución que presenta la varianza muestral. Denición 2.21 (Ji-cuadrado) . Sean Z1, Z2, ... , Zk variables aleatorias siguiendo una distribución Normal estándar (esto es, de media 0 y varianza 1) e independientes. En tal caso se dice que la variable aleatoria X=Z2 1+Z2 2+···+Z2 k sigue una distribución χ2 27 con k grados de libertad y se denota por X∈χ2(k) . Observación 2.22 . A través del Teorema 2.9, para un k grande puede aproximarse χ2(k) por una Normal (k, 2k) . En efecto, pues en base a ese teorema se cumple que k P i=1 X2 i−k √2k d −−−−−−−→ k→∞ Normal (0,1). El motivo por el cual se dice que esta distribución es la distribución de las distancias cuadráticas es por el llamado Teorema de Fisher, que a continuación se presenta. Una demostración del mismo puede encontrarse en las páginas 65 y 66 de [Vélez and García, 2012]. Teorema 2.23 (Teorema de Fisher, [Vélez and García, 2012], p. 66) . Sean X1, X2, ... , Xn variables aleatorias independientes e idénticamente distribuidas como una Normal (µ, σ2) . Entonces la media muestral, X , y la varianza muestral. S2 X son independientes y siguen las siguientes distribuciones: X∈ Normal µ, σ2 n,n S 2 X σ2∈χ2(n−1). Observación 2.24 . En el Teorema de Fisher, a partir de la distribución de la varianza muestral puede obtenerse la de la cuasivarianza muestral S2 cX=1 n−1Pn i=1 Xi−X2 ya que se cumple que n S 2 X σ2=(n−1) S2 cX σ2 . 27 En la Tabla 17 pueden verse algunas de las principales características de esta distribución y en la Figura 11 pueden verse varias funciones de densidad para diferentes grados de libertad.
2.3. Normal 39 Observación 2.25 . Nótese que así como el Teorema Central del Límite proporciona una distribución asintótica, el Teorema de Fisher da una distribución exacta para todo n∈N . Por último, cabe destacar que la distribución χ2 no es más que una distribución Γk 2,1 2 , cuyas principales características pueden verse en la Tabla 18 pero en la que, por una cuestión de extensión del trabajo, no se profundizará. Para un mayor conocimiento de esta distribución puede consultarse el Capítulo 17 de [Johnson et al., 1994]. 2.3.2. t de Student La distribución t de Student fue desarrollada por el matemático inglés William Sealy Gosset (1876-1937) que trabajaba en la cervecera Guiness en Dublín y estaba interesado en estudiar cuestiones como las propiedades químicas de la cebada pero tomando muestras de tamaño muy pequeño (véase [Lehmann, 2012]). Gosset desarrolló esta distribución bajo el pseudónimo de Student en [Student, 1908] y por eso Ronald Aylmer Fisher (1880-1962) en su obra [Fisher, 1925], que fue la que provocó que esta distribución alcanzase una gran fama, se rerió a ella como distribución de Student. La distribución t de Student surge en el contexto de medir la distancia que hay de la media muestral X a la media poblacional µ . Por el Teorema 2.23 X∈ Normal µ, σ2 n o, equivalentemente, X−µ σ/√n∈ Normal (0,1) . Sin embargo, si se desconoce la varianza poblacional σ2 esa distribución no puede ser empleada para analizar la diferencia de medias. La idea que tuvo Gosset fue que, como para muestras grandes σ2 y S2 cX tendrán valores parecidos, podría considerarse el estadístico 28 : t=X−µ ScX/√n. A partir de como está denido t y aplicando el Teorema de Fisher resulta sencillo ver que se trata de un cociente entre una Normal (0,1) y la raíz cuadrada de una χ2 dividida entre sus grados de libertad. En efecto, se cumple que t= X−µ σ/√n q(n−1) S2 cX σ2: (n−1) . En las páginas 67 y 68 de [Vélez and García, 2012] puede verse cómo, suponiendo la independencia entre ambas distribuciones, a partir de las densidades de la Normal y la Ji-cuadrado se obtiene la función de densidad del estadístico, que será la función de densidad de la t de Student. Así, esta distribución se dene como sigue. 28 Véase la Denición I.11.
40 2. Un recorrido por las distribuciones notables continuas Denición 2.26 ( t de Student, [Vélez and García, 2012], p. 68) . Sean dos variables aleatorias Z∈ Normal (0,1) y X∈χ2(k) independientes. En tal caso se tiene que la variable aleatoria X=Z pX/k sigue una distribución t de Student 29 con k grados de libertad y se denota por X∈t(k) . Observación 2.27 . A la vista de la Figura 13 puede verse que la forma que tiene la función de densidad de la t de Student es muy similar a la de la Normal estándar. Sin embargo, existe una diferencia notable entre ambas grácas y es que las colas de la t de Student son más pesadas que las de la Normal estándar, es decir, en la t de Student los valores alejados de la media tienen una frecuencia de aparición mayor que la que tendrían en una Normal estándar. Esto puede comprobarse de manera analítica sin más que observar que, así como las funciones de densidad de ambas tienden a cero cuando x tiende a ±∞ , en la Normal hay un término de orden exponencial, mientras que en la t hay uno de orden potencial. A pesar de la diferencia que se ha comentado en la Observación 2.27 la distribución Normal estándar y la t de Student guardan una relación bastante estrecha que va más allá de tener una forma similar y es que una t de Student de n grados de libertad, cuando n tiende a ∞ , converge en distribución a una Normal estándar. No obstante, para demostrar esto serán necesarios algunos resultados previos cuyas demostraciones no se incluyen en este trabajo pero pueden ser consultadas en [Vélez, 2019]. Teorema 2.28 (Teorema de la aplicación continua, [Vélez, 2019], pp. 289 y 297) . Sea {Xn} una sucesión de variables aleatorias y sea g:R−→ R una aplicación continua. Se verican los siguientes apartados 1) Si Xnc. s. −→ X , entonces g(Xn)c. s. −→ g(X) . 2) Si Xn p −→ X , entonces g(Xn)p −→ g(X) . 3) Si Xnd −→ X , entonces g(Xn)d −→ g(X) . Teorema 2.29 (Teorema de Slutsky, [Vélez, 2019], p. 298) . Sean {Xn} e {Yn} dos sucesiones denidas en el mismo espacio de probabilidad. Si Xnd −→ X e Yn p −→ c∈R , entonces se tiene que Xn+Ynd −→ X+c, XnYnd −→ cX. Teorema 2.30 (Teorema de Khintchine, [Vélez, 2019], p. 328) . Sea {Xn} una sucesión de variables aleatorias independientes, idénticamente distribuidas y con media µ < ∞ . Entonces se 29 En la Tabla 19 pueden verse algunas de las principales características de esta distribución y en la Figura 13 pueden verse varias funciones de densidad para diferentes grados de libertad.
2.3. Normal 41 tiene que 1 n n X i=1 Xi p −−−−→ n→∞ µ. Con estos resultados ya se está en condiciones de demostrar la convergencia de la t de Student a la Normal estándar. Proposición 2.31. Sea {Xn} una sucesión de variables aleatoria tal que Xn∈t(n) . Entonces Xnd −→ Z , donde Z es una variable aleatoria que sigue una distribución Normal estándar. Demostración. Por hipótesis Xn∈t(n) , por lo que se tiene que Xn=Z pYn/n, donde Yn∈χ2(n) . Considérese una sucesión {Zj} de variables aleatorias independientes e idénticamente distribuidas como una Normal estándar y sea {Yn} la sucesión de variables aleatorias tal que Yn= n P i=1 Z2 i∈χ2(n) . Dado que EZ2 i= 1 , por el Teorema 2.30 se cumple que 1 n n X i=1 Z2 i=Yn n p −−−−−−−→ n→∞ 1. En consecuencia, por el Teorema 2.28 se verica que rYn n p −−−−−−−→ n→∞ 1. Por lo tanto, aplicando el Teorema 2.29 se llega a que Xn d −−−−→ n→∞ Z . ■ Para nalizar este apartado cabe destacar que esta distribución se utiliza en una gran variedad de procesos de inferencia estadística como puede ser, por ejemplo, obtener un intervalo de conanza para la media poblacional cuando la varianza es desconocida o en contrastes de hipótesis acerca de la igualdad de medias en dos poblaciones cuyas varianzas son desconocidas 30 . 2.3.3. F de Snédecor La distribución F de Fisher o F de Snédecor debe su nombre a los matemáticos Ronald Aylmer Fisher (1880-1962) y George Waddel Snedecor (1881-1974) que la introdujeron en sus 30 En este trabajo no se profundizará demasiado en estas cuestiones relacionadas con el contraste de hipótesis porque estas son objeto de estudio en la asignatura del tercer curso del Grado en Matemáticas Probabilidad y Estadística .
42 2. Un recorrido por las distribuciones notables continuas obras [Fisher, 1924] y [Snedecor, 1934]. Así como la t de Student surgía de la necesidad de un estadístico que permitiese medir la distancia entre medias, el contexto en el que surge la F de Snédecor es en el de obtener un estadístico que permita comparar las varianzas de dos poblaciones y cuya distribución en el muestreo pueda ser obtenida de manera explícita. El estadístico que se propuso para tal motivo fue F=S2 c1 S2 c2 , donde S2 c1 representa la cuasivarianza muestral de la población 1 y S2 c2 representa la de la población 2. A partir de esa expresión de F , mediante el Teorema de Fisher puede verse que, bajo la hipótesis de que σ2 1=σ2 2 31 , se trata del cociente entre dos distribuciones χ2 divididas entre sus grados de libertad ya que F=(n1−1) S2 c1 σ2 1(n1−1) :(n2−1) S2 c2 σ2 2(n2−1) . Así pues, se trata de encontrar la función de densidad de una variable aleatoria que se obtiene como cociente de dos distribuciones χ2 independientes y tal función de densidad será, por denición, la de la F de Snédecor. En las páginas 75 y 76 de [Vélez and García, 2012] se pueden encontrar los cálculos que permiten llegar a tal función de densidad. Denición 2.32 ( F de Snédecor, [Vélez and García, 2012], p. 76) . Sean dos variables aleatorias X1∈χ2(k1) y X2∈χ2(k2) independientes. En tal caso se tiene que la variable aleatoria X=X1/k1 X2/k2 sigue una distribución F de Snédecor 32 con k1 y k2 grados de libertad y se denota por X∈F(k1, k2). Observación 2.33 . Razonando de manera similar a como se hizo en la Proposición 2.31 se puede ver que la distribución asintótica de F(k1, k2) cuando k2 tiende a ∞ es χ2(k1)/k1 . 31 Esta es la hipótesis nula del contraste en el que surge este estadístico aunque, por la misma razón que en el caso de la t de Student, no se profundizará más en esa cuestión. 32 Algunas de las características principales de esta distribución pueden ser consultadas en la Tabla 20. Además, en la Figura 14 pueden verse diferentes densidades para distintos grados de libertad.
Capítulo 3 Simulación articial de distribuciones La estadística es el único tribunal de apelación para juzgar el nuevo conocimiento Prasanta Chandra Mahalanobis En este capítulo se llevará a cabo una breve introducción a la simulación de variables aleatorias y se presentarán algunos de los procedimientos más habituales ilustrándolos mediante algunos ejemplos 1 . Se introduci¯an algunos métodos generales de simulación de variables, tanto continuas como discretas, y se darán algoritmos especícos para algunas distribuciones notables. Las principales referencias para la elaboración de este capítulo y que pueden ser consultadas para profundizar en las cuestiones que se tratarán a lo largo de él han sido [Devroye, 1986] y [Abad, 2002]. 3.1. Introducción a la simulación A la hora de realizar estudios hay situaciones en las cuales no resulta posible llevar a cabo experimentos directamente sobre la realidad por diversos motivos (lentitud, coste elevado, falta de medios, gran complejidad...) y en las cuales disponer de un modelo que se asemeje al sistema real resulta de vital importancia. En este punto es donde surge la simulación , que es un procedimiento que se basa en experimentar sobre ese modelo del sistema en lugar de sobre la realidad. Tal modelo articial consistirá en un conjunto de variables sometidas a ciertas restricciones que se encuentran relacionadas mediante una serie de ecuaciones matemáticas. No obstante, la obtención de un modelo que se ajuste adecuadamente a la realidad puede ser una 1 No se llevará a cabo un estudio muy pormenorizado debido a que este capítulo se correspondería con la asignatura Simulación Estadística del segundo cuatrimestre del Máster Universitario en Técnicas Estadísticas. 43
50 3. Simulación articial de distribuciones Algoritmo 7 Método de transformación cuantil con búsqueda secuencial 1: Generar U∈ Uniforme (0,1) . 2: Tomar I= 1 . 3: Tomar S=p1 . 4: Mientras U > S hacer 5: I=I+ 1 . 6: S=S+pI . 7: Devolver X=xI . Ejemplo 3.10 ([Abad, 2002], p. 57) . Sea una variable aleatoria X∈ Uniforme Discreta (n) cuya masa de probabilidad viene dada por pi=1 n , para i= 1,2, ... , n . Se tiene que m X i=1 pi≥U > m−1 X i=1 pi⇐⇒ m n≥U > m−1 n⇐⇒ m≥nU > m −1. Esto equivale a que m=⌈nU⌉ . Por lo tanto, el método de simulación que se obtiene es el siguiente: Algoritmo 8 Método de simulación de la Uniforme Discreta (n) 1: Generar U∈ Uniforme (0,1) . 2: Devolver X=⌈nU⌉ . 3.3.2. Métodos de truncamiento Los conocidos como métodos de truncamiento son procedimientos de simulación de variables aleatorias discretas que están basados en la utilización de una distribución continua auxiliar de función de distribución similar a la de la variable que se pretende simular. Supóngase que se desea simular una variable aleatoria X que toma valores x1, x2, ..., xn y considérese la variable I asociada, cuya función de distribución es la función F(x) = P i≤x pi . Además, considérese otra variable aleatoria continua de distribución G(x) y n valores a0=−∞ < a1< a2<··· < an−1< an=∞, tales que para i= 1,2, ... , n se cumple que F(i)−F(i−) = pi=G(ai)−G(ai−1).
3.3. Simulación de distribuciones discretas 51 Así pues, la probabilidad de que I tome el valor i coincide con la probabilidad de que la variable aleatoria continua tome valores en [ai−1, ai) . De esta manera se establece una correspondencia entre los valores de ambas variables, por lo que si la distribución continua se puede simular de manera sencilla (por ejemplo, mediante el método de inversión) tan solo hay que generar valores de ella y luego convertirlos en valores de I . Este razonamiento da lugar al Algoritmo 9, conocido como método de simulación por truncamiento . Algoritmo 9 Método de simulación por truncamiento 1: Generar T con distribución G . 2: Encontrar i tal que ai−1≤T < ai . 3: Devolver I=i . Observación 3.11 . Hay casos especiales donde el Algoritmo 9 especialmente rápido como, por ejemplo, en el caso en que G(0) = 0 y para i= 0,1, ... , n se tiene que ai=i . En esa situación se llega al Algoritmo 10, conocido como método de simulación por truncamiento a la parte entera . Algoritmo 10 Método de simulación por truncamiento a la parte entera 1: Generar T con distribución G . 2: Hacer i=⌈T⌉ . 3: Devolver I=i . Ejemplo 3.12 ([Devroye, 1986], p. 88) . Sea X una variable aleatoria con masa de probabilidad pi=1 ib−1 (i+ 1)b, b > 0, i = 1,2, ... Considérese como variable aleatoria continua auxiliar una cuya función de distribución viene dada por G(x)=1−x−b , para x≥1 y tal que G(1) = 0 . Se verica que G(i+ 1) −G(i) = i−b−(1 + i)−b. Por lo tanto, la variable aleatoria X puede ser simulada mediante ⌈U−1/b⌉ . 3.3.3. Algoritmos para distribuciones especícas En este apartado se especican algunos algoritmos de simulación de dos distribuciones discretas notables que han sido estudiadas a lo largo del trabajo: la distribución Binomial y la Geométrica. En cuanto a la Binomial, considérese una variable aleatoria X∈ Binomial (n, p) . A partir de su denición esta distribución se puede simular mediante el siguiente algoritmo:
52 3. Simulación articial de distribuciones Algoritmo 11 Método de simulación de la Binomial (n, p) 1: Tomar S= 0 . 2: Repetir n veces 3: Generar U∈ Uniforme (0,1) . 4: Si U≤p 5: Tomar S=S+ 1 . 6: Devolver X=S . En cuanto a la distribución Geométrica, también podría simularse directamente a partir de su denición como el número de fracasos hasta que ocurre el primer éxito. Sin embargo, en este caso se optará por presentar un algoritmo de simulación que consiste en un método de truncamiento. Así, sea X∈ Geométrica (p) . Se cumple que P{X=i}=p(1 −p)i, i = 0,1... Considérese como variable continua auxiliar la distribución Exponencial, cuya función de distribución viene dada por la función G(x) = 1 −e−λx , para x≥0 . Entonces se tiene que G(i)−G(i−1) = e−λ(i−1) −e−λi =e−λi−11−e−λ. Por lo tanto, tomando p= 1 −e−λ se cumple que G(i)−G(i−1) = p(1 −p)i−1 . En cuanto a la simulación de la distribución G se tiene que λ=−ln (1 −p) y, además u=G(x) = 1 −e−λx ⇐⇒ −λx = ln (1 −u)⇐⇒ x=−1 λln (1 −u). De este modo, para simular una muestra aleatoria de la distribución Geométrica (p) se puede emplear el siguiente algoritmo: Algoritmo 12 Método de simulación de la Geométrica (p) 1: Tomar L=1 ln(1−p) . 2: Generar U∈ Uniforme (0,1) . 3: Tomar T=L·ln(U) . 4: Devolver X=⌈T⌉ .
Capítulo 4 Aplicaciones con datos reales Es importante comprender lo que puedes hacer antes de aprender a valorar lo bien que parece que lo has hecho. John Wilder Tukey En este capítulo se presentan dos ejemplos reales en los cuales se contrasta el ajuste de dos distribuciones de probabilidad a dos conjuntos de datos. Con ello se pretende ilustrar una idea que se ha ido mostrando a lo largo de todo este trabajo y es que los modelos probabilísticos no son elaboraciones articiales que quedan relegadas al ámbito teórico sino que, aun no siendo modelos exactos, en la vida real tienen un nivel de ajuste elevado y por ello son muy empleados. El código de que se ha empleado para el análisis de los conjuntos de datos puede encontrarse en la Sección IV.1. Se consideran seis conjuntos de datos relativos a las exportaciones e importaciones realizadas en Japón en los años 1995, 2006 y 2017 medidas en yenes 1 . El objetivo es ver si en esos tres conjuntos se verica la Ley de Benford para el primer dígito y para ello se emplearán un test χ2 y un test de Kolmogorov-Smirnov de bondad de ajuste 2 . De este modo, se obtienen un total de doce p -valores, los cuales pueden ser consultados en la Tabla 22. En ella puede verse que únicamente el valor crítico asociado al test de Kolmogorov-Smirnov en las exportaciones del año 2006 y los asociados a las importaciones del año 1995 están por debajo del 5%. Sin embargo, como 1 Esos tres conjuntos de datos han sido extraídos de [Ministry of Finance Japan, 2022] la web del Ministerio de Finanzas del gobierno japonés. 2 Esos dos test estadísticos son objeto de estudio de las materias Probabilidad y Estadística e Inferencia Estadística del tercer curso del grado, por lo que no se profundizará en sus características. No obstante, para consultar algunos acerca de ellos puede consultarse el Capítulo 10 de [Vélez and García, 2012]. Además, los test concretos que se han empleado en el caso de la Ley de Benford han sido extraídos de [Joenssen, 2015]. 53
54 4. Aplicaciones con datos reales el test χ2 de las exportaciones del año 2006 es grande, el único conjunto de datos en el que podría sospecharse que no se verica la Ley de Benford es el relativo a las importaciones del año 1995. Si esos valores fuesen todavía más bajos podría pensarse que hay algún tipo de anomalía en ellos y podría llevarse a cabo alguna investigación más en profundidad para tratar de averiguar el porqué de la discrepancia. No obstante, como tampoco son unos valores excesivamente bajos, se puede armar que no hay evidencias estadísticamente signcativas a un nivel de conanza menor del 1% en contra de suponer que la totalidad de esos conjuntos de datos verican la Ley de Benford. En cuanto a la distribución continua, se considera un conjunto de datos relativos a todos los meteoritos cuyo aterrizaje ha sido detectado en nuestro planeta 3 . El objetivo es tratar de averiguar si los 100 meteoritos de mayor masa (medida en gramos) siguen la distribución de Pareto. En primer lugar, la Figura 15 muestra una representación gráca de la masa de los 100 mayores meteoritos frente a la posición que ocupan en dicho ranking. A la vista de dicha gráca, parece razonable sospechar que la función de densidad de una distribución de Pareto se ajustaría bastante bien a los datos, por lo que se aplica un test de bondad de ajuste para la distribución de Pareto 4 . El resultado de dicho test es un p -valor de 0.9328233, por lo que no hay evidencias estadísticamente signicativas a ningún nivel de conanza razonable en contra de que los 100 meteoritos de mayor tamaño siguen una distribución de Pareto. 3 Este conjunto de datos ha sido extraído de [NASA Open Data Portal, 2018]. 4 El software no tiene ningún contraste para la distribución de Pareto por defecto pero en el paquete ptsuite puede encontrarse uno. Para conocer más detalles acerca de dicho test puede consultarse [Munasinghe et al., 2019].
Anexos
Anexo I Resultados auxiliares Proposición I.1 ([Apostol, 1984]) . Se verica que n X k=1 k=n(n+ 1) 2. Proposición I.2 ([Apostol, 1984]) . Se tiene que: n X k=1 k2=n(n+ 1)(2n+ 1) 6. Proposición I.3 (Fórmula de Stirling, [García et al., 1994]) . Se cumple que l´ım n→∞ n! √2π n n en= 1. Esto es, que cuando n es grande n!≃√2π n n en . Proposición I.4 ([Abramowitz and Stegun, 1972]) . Sean n, k ∈N tales que k≤n . Se tiene que n k=n n−k. Proposición I.5 ([Spivey, 2019]) . Sean m, n, r, s ∈Z+ tales que n≥s . Entonces se cumple que r X k=0 r−k ms+k n=r+s+ 1 m+n+ 1. Denición I.6 (Función gamma, [Abramowitz and Stegun, 1972]) . Dado z∈C tal que ℜ(z)>0 se dene la función gamma de la manera siguiente: Γ(z) = Z∞ 0 tz−1e−tdt. 57
58 I. Resultados auxiliares Denición I.7 (Convergencia en distribución, [Petrov and Mordecki, 2008] p. 119) . Sea {Xn} una sucesión de variables aleatorias y sea X otra variable aleatoria, denidas todas sobre el mismo espacio de probabilidad (Ω,A,P) . Se dice que {Xn} converge en distribución a X , y se escribe Xnd −→ X , si l´ım n→∞Fn(x) = F(x),∀x∈CF, donde Fn es la función de distribución de Xn , F es la de X y CF es el conjunto de puntos de continuidad de F . Denición I.8 (Convergencia en probabilidad, [Petrov and Mordecki, 2008] p. 117) . Sea {Xn} una sucesión de variables aleatorias y sea X otra variable aleatoria, denidas todas sobre el mismo espacio de probabilidad (Ω,A,P) . Se dice que {Xn} converge en probabilidad a X , y se escribe Xn p −→ X , si ∀ε > 0,l´ım n→∞(P{ω∈Ω| |Xn(ω)−X(ω)|< ε})=1. Denición I.9 (Convergencia casi segura, [Petrov and Mordecki, 2008] p.118) . Sea {Xn} una sucesión de variables aleatorias y sea X otra variable aleatoria, denidas todas sobre el mismo espacio de probabilidad (Ω,A,P) . Se dice que {Xn} converge de forma casi segura a X , y se escribe Xnc. s. −→ X , si Pnω∈Ω|l´ım n→∞Xn(ω) = X(ω)|o= 1. Proposición I.10 ([Gut, 2013] p. 213) . Sea {Xn} una sucesión monótona de variables aleatorias denidas sobre el mismo espacio de probabilidad (Ω,A,P) y supóngase que Xn p −→ X , con X una variable aleatoria denida sobre ese mismo espacio. Entonces se tiene que Xnc.s. −→ X . Denición I.11 (Estadístico, [Vélez and García, 2012], p. 15) . Se denomina estadístico a cualquier función T del espacio medible (Ω,A) en un espacio euclídeo (Rk,Bk) que sea medible respecto a las σ -álgebras A y Bk . La dimensión k del espacio euclídeo imagen se denomina dimensión del estadístico . Teorema I.12 (Teorema de Weierstrass, [Bartle and Sherbert, 2011], p. 136) . Sea I= [a, b] un intervalo compacto y sea f:I−→ R una función continua en I . Entonces f alcanza un máximo absoluto y un mínimo absoluto en I .
Anexo II Figuras Sucesos (Ω,A,P) Números (R,B,P∗) X Figura 1: De lo abstracto a lo concreto. Muestreo Con reemplazamiento Sin reemplazamiento k éxitos en n intentos jados Binomial Hipergeométrica n intentos hasta un número de éxitos jados Binomial Negativa Hipergeométrica Negativa B (n, p) H (N, M, n) BN (n, k) HN (N, M, k) Construye Límite Construye Límite Figura 2: Relaciones entre las cuatro distribuciones. Basada en [Miller and Fridell, 2007]. 59
66 II. Figuras Figura 12: Una muestra del universo de distribuciones univariantes. Extraída de [Leemis and McQueston, 2008].
67 Figura 13: Distribuciones t2(k) para diferentes valores de k . Figura 14: Distribuciones F(k1, k2) para diferentes valores de k1 y k2 .
68 II. Figuras Figura 15: Diagrama de dispersión de los 100 mayores meteoritos ordenados por tamaño.
Anexo III Tablas X∈ Uniforme Discreta (n) Parámetros n∈N Soporte {1,2, ... , n} Función masa de probabilidad P{X=k}=1 n Función de distribución P{X≤k}=k n Función característica φ(t) = eit 1−eint n(1 −eit) Media n+ 1 2 Varianza n2−1 12 Tabla 1: La distribución Uniforme Discreta. 69
70 III. Tablas X∈Bernoulli(p) Parámetros p∈(0,1) Soporte {0,1} Función masa de probabilidad P{X=k}=pk(1 −p)1−k Función de distribución P{X≤x}= 1 −p Función característica φ(t) = 1 −p+peit Media p Varianza p(1 −p) Tabla 2: La distribución de Bernoulli.
71 X∈Binomial(n, p) Parámetros n∈N, p ∈(0,1) Soporte {0,1, ... , n} Función masa de probabilidad P{X=k}=n kpk(1 −p)n−k Función de distribución P{X≤x}=⌊x⌋ X i=0 n ipi(1 −p)n−i Función característica φ(t) = (1 −p+peit)n Media np Varianza n p (1 −p) Tabla 3: La distribución Binomial.
72 III. Tablas X∈Binomial Negativa(k, p) Parámetros k∈N, p ∈(0,1) Soporte {0,1,2, ... } Función masa de probabilidad P{X=x}=x+k−1 k−1pk(1 −p)x Función de distribución P{X≤x}= k+x X i=kk+x ipi(1 −p)k+x−i Función característica φ(t) = p 1−(1 −p)eit k , t ∈R Media k(1 −p) p Varianza k(1 −p) p2 Tabla 4: La distribución Binomial Negativa.
73 X∈ Hipergeométrica (N, M, n) Parámetros N∈N, M ∈N∩(0, N), n ∈N∩(0, N) Soporte {L= m´ax{0, n −(N−M)}, L + 1... , m´ın{M, n}} Función masa de probabilidad P{X=x}=M xN−M n−x N n Función de distribución P{X≤x}= x X i=LM iN−M n−i N n Media n p, p =M N Varianza np(1 −p)(N−n) N−1 Tabla 5: La distribución Hipergeométrica.
74 III. Tablas X∈ Hipergeométrica Negativa (N, M, n) Parámetros N∈N, M ∈N∩(0, N), n ∈N∩(0, N) Soporte {0,1... , N −M} Función masa de probabilidad P{X=x}=x+k−1 k−1N−x−k M−k N M Media kM N−M+1 Varianza k(N+1)M (N−M+1)(N−M+2) 1−k N−M+1 Tabla 6: La distribución Hipergeométrica Negativa.
75 X∈ Geométrica (p) Parámetros p∈(0,1) Soporte {0,1,2, ...} Función masa de probabilidad P{X=x}= (1 −p)xp Función de distribución P{X≤x}=p x X k=0 (1 −p)k= 1 −(1 −p)x+1 Función característica φ(t) = p 1−(1 −p)eit , t ∈R Media 1−p p Varianza 1−p p2 Tabla 7: La distribución Geométrica.
82 III. Tablas X∈ Valores Extremos Generalizada (µ, σ, γ) Parámetros µ∈R, σ > 0, γ ∈R\{0} Soporte [µ−σ/γ, +∞) si γ > 0 (−∞, µ −σ/γ ] si γ < 0 Función de densidad f(x) = 1 σ1 + γx−µ σ−(1+1/γ)exp n−1 + γx−µ σ−1/γo Función de distribución F(x) = exp (−1 + γx−µ σ−1/γ) Media µ+σ[ Γ(1 −γ)−1] γ si γ < 1 Varianza σ2hΓ(1 −2γ)−(Γ(1 −γ))2i γ2 si γ < 1 2 Tabla 14: La distribución de Valores Extremos Generalizada con γ= 0 .
83 X∈ Valores Extremos Generalizada (µ, σ, 0) Parámetros µ∈R, σ > 0 Soporte R Función de densidad f(x) = 1 σe−(x−µ σ)exp n−e−(x−µ σ)o Función de distribución F(x) = exp n−e−(x−µ σ)o Función característica φ(t) = eiµt Γ (1 −iσt) Media µ+γbγ∗ Varianza σ2π2 6 ∗bγ es la constante de Euler-Mascheroni (véase [Havil, 2003]). Tabla 15: La distribución de Valores Extremos Generalizada con γ= 0 .
84 III. Tablas X∈ Normal (µ, σ2) Parámetros µ∈R, σ2>0 Soporte R Función de densidad f(x) = 1 σ√2πe−1 2(x−µ σ)2 Función de distribución F(x) = 1 σ√2πZx −∞ e−1 2(x−µ σ)2 dx Función característica φ(t) = eitµ−σ2t2/2 Media µ Varianza σ2 Tabla 16: La distribución Normal.
85 X∈χ2(k) Parámetros k∈N\{0} Soporte R+ Función de densidad f(x) = 1 2k/2Γ(k/2) xk/2−1e−x/2 Función de distribución F(x) = 1 2k/2Γ(k/2) Zx 0 tk/2e−t/2dt Función característica φ(t) = (1 −2it)−k/2 Media k Varianza 2k Tabla 17: La distribución χ2 .
86 III. Tablas X∈Γ(α, β) Parámetros α > 0, β > 0 Soporte R+ Función de densidad f(x) = βα Γ (α)xα−1e−βx Función de distribución F(x) = βα Γ (α)Zx 0 tα−1e−βt dt Media α β Varianza α β2 Tabla 18: La distribución Gamma .
87 X∈t(k) Parámetros k∈N\/{0} Soporte R Función de densidad f(x) = Γn+1 2 √nπ Γn 21 + x2 nn+1 2 Media 0, si k > 1 Varianza k k−2, si k > 2 Tabla 19: La distribución t de Student.
88 III. Tablas X∈F(k1, k2) Parámetros k1∈N\{0}, k2∈N\{0} Soporte R+ Función de densidad f(x) = kk1/2 1kk2/2 2Γk1+k2 2xk1/2−1 Γk1 2Γk2 2(k1x+k2)k1+k2 2 Media k2 k2−2, si k2>2 Varianza 2k2 2(k1+k2−2) k1(k2−2)2(k2−4), si k2>4 Tabla 20: La distribución F de Snédecor.
89 X∈ Uniforme (a, b) Parámetros a, b ∈R Soporte (a, b) Función de densidad f(x) = 1 b−a Función de distribución f(x) = x−a b−a Función característica f(x) = eitb −eita it(b−a) Media a+b 2 Varianza (b−a)2 12 Tabla 21: La distribución Uniforme (a, b) .
90 III. Tablas Exportaciones Importaciones χ2K−Sχ2K−S Año 2017 0.1376467 0.2135000 0.5912239 0.4191000 Año 2006 0.0972103 0.0288000 0.3304694 0.5023000 Año 1995 0.7214202 0.5661000 0.0326048 0.0207000 Tabla 22: Contrastes de bondad de ajuste para la Ley de Benford.
Anexo IV Código de En el presente an se recoge el código que se ha utilizado a lo largo del trabajo. El lenguaje de programación empleado ha sido ([R Core Team, 2022]) y este ha sido escrito mediante el programa Rstudio ([RStudio Team, 2022]). Asimismo, los paquetes que se han utilizado han sido [Xie, 2022], [Joenssen, 2015], [Munasinghe et al., 2019], [Chang, 2022], [Csárdi et al., 2021], [Wolodzko, 2020], [Wickham, 2016] y [Yee, 2022]. IV.1. Código relativo a los test del Capítulo 4 # Se establece el directorio de trabajo setwd("C:/Users/pablo/OneDrive/Escritorio/UNIVERSIDAD/4 º AÑO/TFG/Simulacion") # Se cargan las librerías necesarias library(BenfordTests) library(ptsuite) library(extrafont) library(remotes) library(ggplot2) library(extraDistr) library(VGAM) # Se carga la fuente LM Roman 10 (solo una vez y ya descargado el archivo .ttf) remotes::install_version("Rttf2pt1",version ="1.3.8") extrafont::font_import(pattern ="lm.*") loadfonts(device ="win") 91
98 IV. Código de labels =c("0.00","","0.02","","0.04", "","0.06","","0.08"), limits =c(0.0,0.08), expand =c(0,0)) + scale_x_continuous(breaks =seq(from = xmin+1,to = xmax-1,by =3), labels =seq(from = xmin+1,to = xmax-1,by =3), limits =c(xmin, xmax)) + scale_fill_manual(name =NULL,values =c("bin" ="grey"), labels ="Binomial(100, 0.5)")+ scale_colour_manual(name ="Distribución\n", values =c("normal" ="red"), labels ="Normal(50, 25)") bin_norm; ggsave("binomial_to_normal.jpeg",width =1560,height =1532, units ="px",dpi =320) # Aproximación de la Poisson a la Normal lambda =50 xmin =round(lambda -3*sqrt(lambda)) xmax =round(lambda +3*sqrt(lambda)) df_poi =data.frame(k= xmin:xmax, masa =dpois(xmin:xmax, lambda = lambda)) poi_norm =ggplot(data = df_poi , aes(x=k,y= masa)) + geom_bar(width =0.7,stat ="identity", color ="black",aes(fill ="pois")) + geom_function(fun = dnorm, args =list(mean = lambda, sd =sqrt(lambda)), size =0.8,aes(colour ="normal")) + ggtitle("Aproximación de la Poisson a la Normal\n")+ xlab("\nx")+ylab("Masa de probabilidad o densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =13,hjust =0.5,vjust =1,
IV.2. Código relativo a guras del Anexo II 99 margin =margin(0,0,15,0)), axis.title.x =element_blank(), axis.title.y =element_blank(), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5), axis.ticks.x =element_blank(), legend.title =element_text(family ="LM Roman 10",size =9), legend.text =element_text(family ="LM Roman 10",size =9), legend.key =element_rect(fill ="white"), legend.key.size =unit(0.7,"line"), legend.position =c(0.197,0.89), legend.direction ="vertical", legend.spacing.y =unit(-1.8, ' pt ' ), legend.text.align =0, legend.box.background =element_rect(size =0.3), legend.box.margin =margin(1,1,1,1,unit ="mm")) + guides(fill =guide_legend(byrow =TRUE)) + scale_y_continuous(breaks =seq(from =0.0,to =0.06,0.005), labels =c("0.00","","0.01","","0.02", "","0.03","","0.04","","0.05", "","0.06"), limits =c(0.0,0.06), expand =c(0,0)) + scale_x_continuous(breaks =seq(from = xmin, to = xmax, by =3), labels =seq(from = xmin, to = xmax, by =3), limits =c(xmin-1, xmax+1)) + scale_fill_manual(name =NULL,values =c("pois" ="grey"), labels ="Poisson(50)")+ scale_colour_manual(name ="Distribución\n", values =c("normal" ="red"), labels ="Normal(50, 50) ") poi_norm ggsave("poisson_to_normal.jpeg",width =1560,height =1532, units ="px",dpi =320)
100 IV. Código de # Densidades Ji-cuadrado chi =ggplot() +geom_function(fun = dchisq, args =list(df =1), size =0.8,aes(colour ="1")) + geom_function(fun = dchisq, args =list(df =3), size =0.8,aes(colour ="3")) + geom_function(fun = dchisq, args =list(df =5), size =0.8,aes(colour ="5")) + geom_function(fun = dchisq, args =list(df =9), size =0.8,aes(colour ="9")) + ggtitle("Distribuciones Ji-cuadrado")+xlab("\nx")+ ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =13,hjust =0.5,vjust =1, margin =margin(0,0,15,0)), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5), legend.title =element_text(family ="LM Roman 10",size =10), legend.text =element_text(family ="LM Roman 10",size =10), legend.key =element_rect(fill ="white"), legend.position =c(0.845,0.755), legend.direction ="vertical", legend.spacing.x =unit(-0.5, ' pt ' ), legend.text.align =0, legend.box.background =element_rect(size =0.3),
IV.2. Código relativo a guras del Anexo II 101 legend.box.margin =margin(1,1,1,1,unit ="mm")) + scale_y_continuous(breaks =seq(from =0.0,to =0.30,0.05), labels =c("0.00","","0.10","","0.20","","0.30"), limits =c(0.0,0.30), expand =c(0.02,0.01)) + scale_x_continuous(breaks =seq(from =0,to =8,by =1), labels =seq(from =0,to =8,by =1), limits =c(-0.1,8.1), expand =c(0.03,0.03)) + scale_colour_manual(name ="Distribución",values =c("1" ="orange", "3" ="blue", "5" ="purple", "9" ="green"), labels =expression(" k = 1"," k = 3","k=5", " k = 9")) chi; ggsave("chi_dens.jpeg",width =1560,height =1532, units ="px",dpi =320) # Densidades t de Student y Normal(0, 1) t_dens =ggplot() +geom_function(fun = dt, args =list(df =3), size =0.8,aes(colour ="3")) + geom_function(fun = dt, args =list(df =5), size =0.8, aes(colour ="5")) + geom_function(fun = dt, args =list(df =20), size =0.8,aes(colour ="20")) + geom_function(fun = dnorm, args =list(mean =0,sd =1), size =0.8,aes(colour ="normal")) + ggtitle("Distribuciones t de Student")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =13,hjust =0.5,vjust =1,
102 IV. Código de margin =margin(0,0,15,0)), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5), legend.title =element_text(family ="LM Roman 10",size =10), legend.text =element_text(family ="LM Roman 10",size =10), legend.key =element_rect(fill ="white"), legend.position =c(0.81,0.775), legend.direction ="vertical", legend.spacing.x =unit(-0.5, ' pt ' ), legend.text.align =0, legend.box.background =element_rect(size =0.3), legend.box.margin =margin(1,1,1,1,unit ="mm")) + scale_y_continuous(breaks =seq(from =0.0,to =0.40,0.05), labels =c("0.00","","0.10","","0.20","", "0.30","","0.40"), limits =c(0.0,0.40), expand =c(0.02,0.01)) + scale_x_continuous(breaks =seq(from =-4,to =4,by =1), labels =seq(from =-4,to =4,by =1), limits =c(-4.1,4.1), expand =c(0.03,0.03)) + scale_colour_manual(name ="Distribución",values =c("3" ="orange", "5" ="blue", "20" ="green", "normal" ="red"), labels =expression(" k = 3"," k = 5"," k = 20", " Normal(0, 1)")) t_dens; ggsave("student_dens.jpeg",width =1560,height =1532, units ="px",dpi =320)
IV.2. Código relativo a guras del Anexo II 103 # Densidades F de Snédecor F_dens =ggplot() +geom_function(fun = df, args =list(df1 =3,df2 =1), size =0.8,aes(colour ="31")) + geom_function(fun = df, args =list(df1 =7,df2 =5), size =0.8,aes(colour ="75")) + geom_function(fun = df, args =list(df1 =20,df2 =12), size =0.8,aes(colour ="2012")) + ggtitle("Distribuciones F de Snédecor")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =13,hjust =0.5,vjust =1, margin =margin(0,0,20,0)), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5), legend.title =element_text(family ="LM Roman 10",size =10), legend.text =element_text(family ="LM Roman 10",size =10), legend.key =element_rect(fill ="white"), legend.position =c(0.825,0.815), legend.direction ="vertical", legend.spacing.x =unit(-0.5, ' pt ' ), legend.text.align =0, legend.box.background =element_rect(size =0.3), legend.box.margin =margin(1,1,1,1,unit ="mm")) + scale_y_continuous(breaks =seq(from =0.0,to =1,0.125),
104 IV. Código de labels =c("0.00","","0.25","","0.50", "","0.75","","1.00"), limits =c(0.0,1.00), expand =c(0.02,0.01)) + scale_x_continuous(breaks =seq(from =0,to =4,by =1), labels =seq(from =0,to =4,by =1), limits =c(-0.1,4.1), expand =c(0.03,0.03)) + scale_colour_manual(name ="Distribución",values =c("31" ="orange", "75" ="blue", "2012" ="green"), labels =expression(" (3, 1)"," (7, 5)", " (20, 12)")) F_dens; ggsave("snedecor_dens.jpeg",width =1560,height =1532, units ="px",dpi =320) IV.3. Código relativo a las tablas del Apéndice III # 1. Distribuciones Discretas # Uniforme Discreta df_ud =data.frame(k=1:10,masa =ddunif(1:10,min =1,max =10)) ud =ggplot(data = df_ud, aes(x= k, y= masa)) + geom_bar(width =0.8,stat ="identity",color ="black",fill ="grey")+ ggtitle("Uniforme Discreta(10)")+ xlab("\nk")+ylab("Función masa de probabilidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =14,hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5),
IV.3. Código relativo a las tablas del Apéndice III 105 axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5), axis.ticks.x =element_blank()) + scale_y_continuous(breaks =seq(from =0.0,to =0.12,0.01), labels =c("0.00","","0.02","","0.04","","0.06", "","0.08","","0.10","","0.12"), limits =c(0.0,0.12), expand =c(0,0)) + scale_x_continuous(breaks =1:10,labels =1:10) ud ggsave("uniforme_discreta.jpeg",width =1130,height =987, units ="px",dpi =320) # Bernoulli df_bern =data.frame(k=c(0,1), masa =dbern(c(0,1), 0.4)) bern =ggplot(data = df_bern , aes(x=k,y= masa)) + geom_bar(width =0.8,stat ="identity",color ="black",fill ="grey")+ ggtitle("Bernoulli(0.4)\n")+ xlab("\nÉxitos")+ylab("Función masa de probabilidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =14,hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1),
106 IV. Código de axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5), axis.ticks.x =element_blank()) + scale_y_continuous(breaks =seq(from =0.0,to =0.8,0.1), labels =c("0.00","","0.02","","0.04","","0.06", "","0.08"), limits =c(0.0,0.8), expand =c(0,0)) + scale_x_continuous(breaks =c(0,1), labels =c(0,1)) bern; ggsave("bernoulli.jpeg",width =1130,height =987, units ="px",dpi =320) # Binomial df_bin =data.frame(k=0:30,masa =dbinom(0:30,size =30,prob =0.5)) bin =ggplot(data = df_bin , aes(x=k,y= masa)) + geom_bar(width =0.7,stat ="identity",color ="black",fill ="grey")+ ggtitle("Binomial(30, 0.5)\n")+ xlab("\nÉxitos en 30 intentos")+ylab("Función masa de probabilidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =14,hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5), axis.ticks.x =element_blank()) +
IV.3. Código relativo a las tablas del Apéndice III 107 scale_y_continuous(breaks =seq(from =0.0,to =0.175,0.025), labels =c("0.00","","0.05","","0.10","","0.15", ""), limits =c(0.0,0.16), expand =c(0,0)) + scale_x_continuous(breaks =seq(from =0,to =30,by =3), labels =seq(from =0,to =30,by =3)) bin ggsave("binomial.jpeg",width =1130,height =987, units ="px",dpi =320) # Binomial Negativa df_nbin =data.frame(k=0:30,masa =dnbinom(0:30,size =4,prob =.25)) nbin =ggplot(data = df_nbin , aes(x=k,y= masa)) + geom_bar(width =0.7,stat ="identity",color ="black",fill ="grey")+ ggtitle("Binomial Negativa(4, 0.25)\n")+ xlab("\nFallos hasta el cuarto éxito")+ ylab("Función masa de probabilidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =14,hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5), axis.ticks.x =element_blank()) +
114 IV. Código de zipf; ggsave("zipf.jpeg",width =1130,height =987, units ="px",dpi =320) # 2. Distribuciones continuas # Pareto par =ggplot() + geom_function(fun = dpareto, args =list(location =1,shape =2), size =0.8)+ ggtitle("Pareto(2, 1)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10", face ="bold",size =13, hjust =0.5,vjust =1, margin =margin(0,0,15,0)), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =2.5,0.25), labels =c("0.00","","0.50","","1.00","","1.50", "","2.00","","2.50"), limits =c(0.0,2.5), expand =c(0.005,0.005)) + scale_x_continuous(breaks =seq(from =0,to =5,by =1), labels =seq(from =0,to =5,by =1), limits =c(1,5), expand =c(0.03,0.03))
IV.3. Código relativo a las tablas del Apéndice III 115 par ggsave("pareto.jpeg",width =1130,height =987,units ="px",dpi =320) # Weibull weib =ggplot() + geom_function(fun = dweibull, args =list(shape =2,scale =0.5), size =0.8)+ ggtitle("Weibull(2, 2)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =14,hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =1.75,0.25), labels =c("0.0","","0.5","","1.0","","1.5", ""), limits =c(0.0,1.75), expand =c(0.02,0.01)) + scale_x_continuous(breaks =seq(from =0,to =3,by =0.5), labels =seq(from =0,to =3,by =0.5), limits =c(0.0,3.0), expand =c(0.03,0.03)) weib; ggsave("weibull.jpeg",width =1130,height =987, units ="px",dpi =320)
116 IV. Código de # Exponencial exp =ggplot() + geom_function(fun = dexp, args =list(rate =1.5), size =0.8)+ ggtitle("Exponencial(1.5)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =14,hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =1.5,0.5), labels =c("0.0","","0.5","","1.0","","1.5"), limits =c(0.0,1.5), expand =c(0.02,0.01)) + scale_x_continuous(breaks =seq(from =0,to =6,by =1), labels =seq(from =0,to =6,by =1), limits =c(0.0,6.0), expand =c(0.03,0.03)) exp; ggsave("exponencial.jpeg",width =1130,height =987, units ="px",dpi =320) # Valores Extremos dev =function(mu,sigma,gamma){ if(gamma!=0){ function(x){
IV.3. Código relativo a las tablas del Apéndice III 117 1/sigma*(1+gamma*(x-mu)/sigma)** (-(1+1/gamma))* exp(-(1+gamma*(x-mu)/sigma)**(-1/gamma)) } } else{ function(x){ exp(-(x-mu)/sigma)*exp(-exp(-(x-mu)/sigma)) } } } evg =ggplot() + geom_function(fun =dev(1,2,1/4), size =0.8)+ ggtitle("Valores Extremos(1, 2, 1/4)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =13,hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =0.2,0.025), labels =c("0.00","","0.05","","0.01","","0.15", "","0.20"), limits =c(0.0,0.2), expand =c(0.005,0.005)) + scale_x_continuous(breaks =seq(from =-7,to =11,by =2),
118 IV. Código de labels =seq(from =-7,to =11,by =2), limits =c(-6.9,11), expand =c(0.03,0.03)) evg; ggsave("valores_extremos_generalizada.jpeg",width =1130,height =987, units ="px",dpi =320) ev_0 =ggplot() + geom_function(fun =dev(1,1,0), size =0.8)+ ggtitle("Valores Extremos(1, 1, 0)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10",face ="bold", size =13,hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =0.4,0.05), labels =c("0.00","","0.10","","0.20","","0.30", "","0.40"), limits =c(0.0,0.4), expand =c(0.005,0.005)) + scale_x_continuous(breaks =seq(from =-2,to =4,by =1), labels =seq(from =-2,to =4,by =1), limits =c(-2,4), expand =c(0.03,0.03)) ev_0; ggsave("valores_extremos_gumbel.jpeg",width =1130,height =987, units ="px",dpi =320)
IV.3. Código relativo a las tablas del Apéndice III 119 # Normal norm =ggplot() + geom_function(fun = dnorm, args =list(mean =0,sd =1), size =0.8)+ ggtitle("Normal(0, 1)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10", face ="bold",size =13, hjust =0.5,vjust =1, margin =margin(0,0,15,0)), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =0.4,0.05), labels =c("0.00","","0.10","","0.20", "","0.30","","0.40"), limits =c(0.0,0.4), expand =c(0.005,0.005)) + scale_x_continuous(breaks =seq(from =-3,to =3,by =1), labels =seq(from =-3,to =3,by =1), limits =c(-3.1,3.1), expand =c(0.03,0.03)) norm; ggsave("norm.jpeg",width =1130,height =987, units ="px",dpi =320)
120 IV. Código de # Chi cuadrado chisq =ggplot() + geom_function(fun = dchisq, args =list(df =7), size =0.8)+ ggtitle("Ji-cuadrado(7)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10", face ="bold",size =13, hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =0.15,0.025), labels =c("0.00","","0.05","","0.10", "","0.15"), limits =c(0.0,0.15), expand =c(0.005,0.005)) + scale_x_continuous(breaks =seq(from =0,to =20,by =2), labels =seq(from =0,to =20,by =2), limits =c(0,20), expand =c(0.03,0.03)) chisq; ggsave("chisquared.jpeg",width =1130,height =987, units ="px",dpi =320) # Gamma
IV.3. Código relativo a las tablas del Apéndice III 121 gam =ggplot() + geom_function(fun = dgamma, args =list(shape =6,scale =1.25), size =0.8)+ ggtitle("Gamma(6, 1.25)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10", face ="bold",size =13, hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =0.15,0.025), labels =c("0.00","","0.05","","0.10", "","0.15"), limits =c(0.0,0.15), expand =c(0.005,0.005)) + scale_x_continuous(breaks =seq(from =0,to =20,by =2), labels =seq(from =0,to =20,by =2), limits =c(0,20), expand =c(0.03,0.03)) gam; ggsave("gamma.jpeg",width =1130,height =987, units ="px",dpi =320) # t de Student t=ggplot() +
122 IV. Código de geom_function(fun = dt, args =list(df =5), size =0.8)+ ggtitle("t(5)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10", face ="bold",size =13, hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =0.4,0.05), labels =c("0.00","","0.10","","0.20", "","0.30","","0.40"), limits =c(0.0,0.40), expand =c(0.005,0.005)) + scale_x_continuous(breaks =seq(from =-5,to =5,by =1), labels =seq(from =-5,to =5,by =1), limits =c(-5,5), expand =c(0.03,0.03)) t; ggsave("student.jpeg",width =1130,height =987, units ="px",dpi =320) # F de Fisher f=ggplot() + geom_function(fun = df, args =list(df1 =5,df2 =2), size =0.8)+
IV.3. Código relativo a las tablas del Apéndice III 123 ggtitle("F(5, 2)")+ xlab("\nx")+ylab("Función de densidad\n")+ theme(panel.background =element_blank(), panel.grid.major =element_blank(), panel.grid.minor =element_blank(), axis.line.y =element_line(color ="black"), axis.line.x =element_line(color ="black"), plot.title =element_text(family ="LM Roman 10", face ="bold",size =13, hjust =0.5,vjust =1), axis.title.x =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =0.5), axis.title.y =element_text(family ="LM Roman 10", face ="bold",size =11, hjust =0.5,vjust =1), axis.text.x =element_text(family ="LM Roman 10",size =11), axis.text.y =element_text(family ="LM Roman 10",size =11, angle =90,hjust =0.5)) + scale_y_continuous(breaks =seq(from =0.0,to =0.6,0.1), labels =c("0.00","","0.20","","0.40", "","0.60"), limits =c(0.0,0.6), expand =c(0.005,0.005)) + scale_x_continuous(breaks =seq(from =0,to =5,by =1), labels =seq(from =0,to =5,by =1), limits =c(0,5), expand =c(0.03,0.03)) f;ggsave("snedecor.jpeg",width =1130,height =987, units ="px",dpi =320) # Uniforme continua uni_cont =ggplot() + geom_function(fun = dunif, args =list(min =0,max =1), size =0.8)+ ggtitle("Uniforme Continua(0, 1)")+ xlab("\nx")+ylab("Función de densidad\n")+
130 BIBLIOGRAFÍA [Krishnamoorthy, 2016] Krishnamoorthy, K. (2016). Handbook of statistical distributions with applications . Chapman and Hall/CRC. [Leemis and McQueston, 2008] Leemis, L. M. and McQueston, J. T. (2008). Univariate distribution relationships. The American Statistician , 62(1):4553. [Lehmann, 2012] Lehmann, E. L. (2012). student and small-sample theory. In Selected works of EL Lehmann , pages 9971004. Springer. [Li, 2002] Li, W. (2002). Zipf's law everywhere. Glottometrics , 5:1421. [Lu et al., 2016] Lu, Y., Miller, A. A., Homann, R., and Johnson, C. W. (2016). Towards the automated verication of Weibull distributions for system failure rates. In Critical Systems: Formal Methods and Automated Verication , pages 8196. Springer. [Maxwell, 1860] Maxwell, J. C. (1860). Illustrations of the dynamical theory of gases. part i. on the motions and collisions of perfectly elastic spheres. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science , 19(124):1932. [Mebane et al., 2008] Mebane, W., Álvarez, R. M., Hall, T. E., and Hyde, S. D. (2008). Election forensics: The second-digit Benford's law test and recent american presidential elections. Election fraud: detecting and deterring electoral manipulation , pages 162181. [Miller and Fridell, 2007] Miller, G. K. and Fridell, S. L. (2007). A forgotten discrete distribution? Reviving the Negative Hypergeometric model. The American Statistician , 61(4):347350. [Ministry of Finance Japan, 2022] Ministry of Finance Japan (2022). Trade Statistics of Japan: Values by Country. Recuperado el 3 de julio de 2022 de https://www.customs.go.jp/toukei/ srch/indexe.htm?M=23&P=0. [Moore et al., 2015] Moore, D., Notz, W., and Fligner, M. (2015). The Basic Practice of Statistics . W. H. Freeman. [Munasinghe et al., 2019] Munasinghe, R., Kossinna, P., Jayasinghe, D., and Wijeratne, D. (2019). ptsuite: Tail Index Estimation for Power Law Distributions . R package version 1.0.0. [Muraleedharan et al., 2011] Muraleedharan, G., Soares, C. G., and Lucas, C. (2011). Characteristic and moment generating functions of generalised extreme value distribution (GEV). Sea Level Rise, Coastal Engineering, Shorelines and Tides , pages 269276. [Nadarajah, 2006] Nadarajah, S. (2006). The Exponentiated Gumbel distribution with climate application. Environmetrics: The ocial journal of the International Environmetrics Society , 17(1):1323.
BIBLIOGRAFÍA 131 [Nadarajah and Pogány, 2013] Nadarajah, S. and Pogány, T. K. (2013). On the characteristic functions for extreme value distributions. Extremes , 16(1):2738. [NASA Open Data Portal, 2018] NASA Open Data Portal (2018). Meteorite Landings . Recuperado el 3 de julio de 2022 de https://data.nasa.gov/Space-Science/Meteorite-Landings/ gh4g-9sfh. [Newcomb, 1881] Newcomb, S. (1881). Note on the frequency of use of the dierent digits in natural numbers. American Journal of mathematics , 4(1):3940. [Nigrini, 1996] Nigrini, M. J. (1996). A taxpayer compliance application of Benford's law. The Journal of the American Taxation Association , 18(1):72. [Papoulis and Pillai, 2002] Papoulis, A. and Pillai, S. U. (2002). Probability, random variables, and stochastic processes . Tata McGraw-Hill Education. [Park et al., 2013] Park, S., Serpedin, E., and Qaraqe, K. (2013). Gaussian assumption: The least favorable but the most useful [lecture notes]. IEEE Signal Processing Magazine , 30(3):183186. [Pearson, 1926] Pearson, K. (1926). Abraham de moivre. Nature , 117:551552. [Petrov and Mordecki, 2008] Petrov, V. and Mordecki, E. (2008). Teoría de la probabilidad . DIRAC. [Piqueira et al., 1999] Piqueira, J., Monteiro, L., De Magalhães, T., Ramos, R., Sassi, R., and Cruz, E. (1999). Zipf's law organizes a psychiatric ward. Journal of theoretical biology , 198(3):439443. [Poisson, 1837] Poisson, S. (1837). Recherches sur la probabilité des jugements en matière criminelle et en matière civile . Bachelier. [Prokhorov and Statulevi£ius, 2000] Prokhorov, Y. V. and Statulevi£ius, V. (2000). Limit theorems of probability theory . Springer. [R Core Team, 2022] R Core Team (2022). R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing, Vienna, Austria. [Ramanathan, 1993] Ramanathan, R. (1993). Statistical Methods in Econometrics . Emerald Group Publishing Limited. [Ridout, 1999] Ridout, M. S. (1999). Memory in coal tits: An alternative model. Biometrics , 55(2):660662. [Rocco, 2014] Rocco, M. (2014). Extreme value theory in nance: A survey. Journal of Economic Surveys , 28(1):82108.
132 BIBLIOGRAFÍA [Rosin et al., 1933] Rosin, P., Rammler, E., and Sperling, K. (1933). Size of powdered coal and its meaning to grinding (in german). Bericht C 52 des Reichskohlenrats . [Ross, 2019] Ross, S. M. (2019). A rst course in probability . Pearson Boston. [RStudio Team, 2022] RStudio Team (2022). RStudio: Integrated Development Environment for R . RStudio, PBC, Boston, MA. [Salman, 1999] Salman, S. R. (1999). Chemical Reactions Studied by Electronic Spectroscopy. [Santarelli et al., 2017] Santarelli, M. F., Positano, V., and Landini, L. (2017). Measured PET data characterization with the Negative Binomial distribution model. Journal of Medical and Biological Engineering , 37(3):299312. [Seal, 1967] Seal, H. L. (1967). Studies in the History of Probability and Statistics. XV The historical development of the Gauss linear model. Biometrika , 54(1-2):124. [Shapiro and Wilk, 1965] Shapiro, S. S. and Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika , 52(3/4):591611. [Snedecor, 1934] Snedecor, G. W. (1934). Calculation and interpretation of analysis of variance and covariance. [Spivey, 2019] Spivey, M. Z. (2019). The art of proving binomial identities . CRC Press. [Spyrou et al., 2009] Spyrou, C., Stark, R., Lynch, A. G., and Tavaré, S. (2009). BayesPeak: Bayesian analysis of ChIP-seq data. BMC bioinformatics , 10(1):117. [Stanislevic and VoteTrustUSA, 2006] Stanislevic, H. and VoteTrustUSA, E. (2006). Random auditing of e-voting systems: How much is enough. National Coalition for Election Integrity, E-Voter Education Project . [Stigler, 1986] Stigler, S. M. (1986). The history of statistics: The measurement of uncertainty before 1900 . Harvard University Press. [Stoica and Babu, 2011] Stoica, P. and Babu, P. (2011). The gaussian data assumption leads to the largest Cramér-Rao bound [lecture notes]. IEEE Signal Processing Magazine , 28(3):132 133. [Student, 1908] Student (1908). The probable error of a mean. Biometrika , pages 125. [Sutclie, 1978] Sutclie, J. V. (1978). Methods of ood estimation: A guide to the ood studies report. [Tam Cho and Gaines, 2007] Tam Cho, W. K. and Gaines, B. J. (2007). Breaking the (Benford) law: Statistical fraud detection in campaign nance. The american statistician , 61(3):218223.
BIBLIOGRAFÍA 133 [Tojeiro et al., 2014] Tojeiro, C., Louzada, F., Roman, M., and Borges, P. (2014). The complementary Weibull geometric distribution. Journal of Statistical Computation and Simulation , 84(6):13451362. [Tse et al., 2014] Tse, K.-K. et al. (2014). Some applications of the Poisson process. Applied Mathematics , 5(19):3011. [Tsu et al., 1951] Tsu, T. C., Mugele, R. A., and McClintock, F. A. (1951). Discussion: A statistical distribution function of wide applicability. Journal of Applied Mechanics , 18(2):293297. [Van Slambrouck et al., 2014] Van Slambrouck, K., Stute, S., Comtat, C., Sibomana, M., van Velden, F. H., Boellaard, R., and Nuyts, J. (2014). Bias reduction for low-statistics PET: maximum likelihood reconstruction with a modied Poisson distribution. IEEE Transactions on medical imaging , 34(1):126136. [Vélez, 2019] Vélez, R. (2019). Cálculo de Probabilidades 2. UNED. [von Bortkiewicz, 1898] von Bortkiewicz, L. (1898). Das Gesetz der kleinen Zahlen . B.G. Teubner. [von Bortkiewicz, 1922] von Bortkiewicz, L. (1922). Range and mean error (in german). Sitzungsberichte der BerlinerMathematischen Gesellschaft , 27:333. [von Mises, 1923] von Mises, R. (1923). On the range in a series of observations (in german). Sitzungsberichte derBerliner Mathematischen Gesellschaft , 222:38. [Vélez and García, 2012] Vélez, R. and García, A. (2012). Principios de Inferencia Estadística . UNED. [Vélez and Hernández, 1995] Vélez, R. and Hernández, V. (1995). Dados, monedas y urnas . UNED. [Weibull, 1939a] Weibull, W. (1939a). The phenomenon of rupture in solids. IVA Handlingar , 153. [Weibull, 1939b] Weibull, W. (1939b). A Statistical Theory of the Strength of Materials . Number 151 in Handlingar / Ingeniörsvetenskapsakademien. Generalstabens litograska anstalts förlag. [Weibull, 1951] Weibull, W. (1951). A statistical distribution function of wide applicability. Journal of applied mechanics , (18):293297. [Westermeier and Michaelis, 1995] Westermeier, T. and Michaelis, J. (1995). Applicability of the Poisson distribution to model the data of the German children's cancer registry. Radiation and Environmental Biophysics , 34(1):711.
134 BIBLIOGRAFÍA [Westfall, 2014] Westfall, P. H. (2014). Kurtosis as peakedness, 19052014. RIP. The American Statistician , 68(3):191195. [Wickham, 2016] Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis . Springer- Verlag New York. [Wolodzko, 2020] Wolodzko, T. (2020). extraDistr: Additional Univariate and Multivariate Distributions . R package version 1.9.1. [Xie, 2022] Xie, Y. (2022). knitr: A General-Purpose Package for Dynamic Report Generation in R . R package version 1.39. [Yee, 2022] Yee, T. W. (2022). VGAM: Vector Generalized Linear and Additive Models . R package version 1.1-6. [Yevjevich, 1972] Yevjevich, V. M. (1972). Probability and statistics in hydrology. [Zhu et al., 2018] Zhu, Y., Zhang, B., Wang, Q. A., Li, W., and Cai, X. (2018). The principle of least eort and Zipf distribution. Journal of Physics: Conference Series , 1113(1):012007. [Zipf, 1949] Zipf, G. K. (1949). Human behavior and the principle of least eort . Addison Wesley Press.