Inferencia bayesiana

TRABAJO FIN DE GRADO Inferencia Bayesiana Realizado por: Íñigo Ortiz Padilla Tutores: DR. RAFAEL BLANQUERO BRAVO DR. EMILIO CARRIZOSA PRIEGO FACULTAD DE MATEMÁTICAS Departamento de Estadística e Investigación Operativa Sevilla, 2018 Índice general Abstract 5 1. Conceptos previos 7 1.1. Propiedades de la inferencia bayesiana . . . . . . . . . . . . . . . . . 7 1.2. FórmuladeBayes............................ 7 1.2.1. Formulación del teorema . . . . . . . . . . . . . . . . . . . . 7 1.2.2. Fórmula de Bayes y función de verosimilitud . . . . . . . . . 8 1.2.3. Naturaleza secuencial de la fórmula de Bayes . . . . . . . . . 9 1.3. Probabilidad subjetiva . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4. Procedimiento a seguir . . . . . . . . . . . . . . . . . . . . . . . . . 10 2. Herramientas de la Inferencia Bayesiana 13 2.1. Estimación bayesiana . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.1. Estimación puntual . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.2. Intervalos............................ 15 2.2. Elección de la distribución a priori . . . . . . . . . . . . . . . . . . . 18 2.2.1. Distribuciones a priori conjugadas . . . . . . . . . . . . . . . 18 2.2.2. Distribuciones a priori impropias . . . . . . . . . . . . . . . . 19 2.2.3. Distribuciones a priori no informativas: Regla de Jeffreys . . . 19 2.3. Distribución a posteriori predictiva . . . . . . . . . . . . . . . . . . . 21 3. Inferencia bayesiana sobre una proporción 23 3.1. Cálculo de la distribución a posteriori . . . . . . . . . . . . . . . . . 23 3.1.1. Distribución a priori discreta para ˜p.............. 23 3.1.2. Distribución a priori continua para ˜p.............. 24 3.1.2.1. Distribución continua a priori para ˜p........ 24 3.1.2.2. Distribución continua uniforme para ˜p....... 25 3.1.3. Distribución Beta a priori . . . . . . . . . . . . . . . . . . . . 27 3.2. RegladeJeffreys ............................ 28 3.3. Distribución predictiva . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4. Contraste de hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3 4 Índice general 4. Distribución de Poisson 31 4.1. Función de verosimilitud . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2. Cálculo de la distribución a posteriori . . . . . . . . . . . . . . . . . 32 4.2.1. Distribución a priori no informativa . . . . . . . . . . . . . . 33 4.2.2. Distribución Gamma a priori para Poisson ˜ λ.......... 33 4.3. Distribución predictiva . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.4. Contraste de hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5. Distribución normal 37 5.1. Función de verosimilitud para (µ, σ2)................. 37 5.2. RegladeJeffreys ............................ 39 5.2.1. Caso ˜µdesconocido y σconocido............... 40 5.2.2. Caso µconocido, ˜σdesconocido................ 40 5.3. Cálculo de la distribución a posteriori . . . . . . . . . . . . . . . . . 41 5.3.1. Caso ˜µdesconocido y σ2conocido............... 41 5.3.1.1. Distribución a priori no informativa . . . . . . . . . 41 5.3.1.2. Distribución a priori conjugada para ˜µ....... 42 5.3.2. Caso µconocido y ˜σ2desconocido . . . . . . . . . . . . . . . 42 5.3.2.1. Distribución a priori no informativa . . . . . . . . . 42 5.3.2.2. Distribución a priori conjugada . . . . . . . . . . . 47 5.3.3. Caso ˜µ, ˜σdesconocidos .................... 47 5.3.3.1. Distribución a priori no informativa . . . . . . . . . 47 5.3.3.2. Distribución conjugada . . . . . . . . . . . . . . . 51 5.4. Distribución predictiva . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.5. Testdehipótesis............................. 53 5.5.1. Media de una normal con varianza desconocida . . . . . . . . 53 5.5.2. Varianza............................. 54 6. Inferencia Bayesiana en R 55 6.1. Inferencia sobre una proporción . . . . . . . . . . . . . . . . . . . . 55 6.2. Distribución de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.3. Distribución Normal . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.4. Ejemplopráctico ............................ 60 Abstract The subject of this document is Bayesian Inference, an inference system based on Bayes’ Formula. In the first chapter we will state this formula and will discuss how to use it. It will be shown that, according to the formula, posterior distributions are proportional to the likelihood function times the prior distribution. This is the essential notion of Bayesian inference. In the second chapter, definitions and results related to Bayesian Analysis will be given in order to accomplish our inference. These mathematical tools will be of great use in the following chapters, where inference on proportions, Poisson distribution and Normal distribution will be studied. Both conjugate and noninformative prior distributions are considered. This documents ends with the transcription of an R code, allowing us to compute posterior distributions. 5 Capítulo 1 Conceptos previos 1.1. Propiedades de la inferencia bayesiana La inferencia bayesiana es una estrategia de inferencia estadística que tiene su base en el teorema de Bayes. Algunas de sus propiedades, enumeradas a continuación, la hacen particularmente útil para ser usada en la investigación científica. Dado un modelo, la inferencia bayesiana hace uso de los datos obtenidos empíricamente. Inferencias que resulten, según nuestros propósitos, inaceptables han de ser consecuencia de afirmaciones inapropiadas, no de un inadecuado sistema de inferencia. Por tanto, todas los factores del modelo, incluidas las distribuiones a priori, son susceptibles de ser modificadas. Ya que este sistema de inferencia puede ser aplicado a casi cualquier modelo de probabilidad, no se requiere prestar tanta atención a la conveniencia del modelo considerado como al valor científico de éste. 1.2. Fórmula de Bayes Nos detendremos ahora en analizar un importante teorema que da lugar a la materia de estudio de este trabajo. Aparecido en un artículo del 1764 firmado por Thomas Bayes, este resultado es designado de manera usual como fórmula de Bayes. 1.2.1. Formulación del teorema Sea x= (x1, ..., xn)un vector de nobservaciones cuya distribución de probabilidad P(x|θ)depende de los valores de kparámetros θ0= (θ1, ..., θk). Supongamos además que θtiene distribución de probalibilidad P(θ). Entonces se tiene: 7 8 1.2. Fórmula de Bayes P(x|θ)P(θ) = P(x, θ) = P(θ|x)P(x) Dado el vector de datos observados x, la distribución condicionada de θaxes la siguiente: P(θ|x) = P(x|θ)P(θ) P(x). Además podemos escribir: P(x) = E[P(x|θ)] = c−1. Para los casos continuo y dicreto de θtenemos entonces, respectivamente: P(x) = RP(x|θ)P(θ)dθ P(x) = PP(x|θ)P(θ) donde la suma y la integral están definidas en el dominio de θ, y la esperanza está tomada con respecto a la distribución P(θ). Podemos escribir en consecuencia P(θ|x) = cP(x|θ)P(θ). En la expresión precedente, P(θ), conocida como distribución a priori de θ, nos dice qué sabemos de θsin conocimiento de los datos muestrales, mientras que P(θ|x) nos dice qué sabemos de θuna vez conocidos los datos muestrales, y es por ello llamada distribución a posteriori de θdado x. Por otra parte, la constante ctiene como utilidad normalizar la suma o la integral para que sea igual a uno. 1.2.2. Fórmula de Bayes y función de verosimilitud Supongamos ahora que conocemos los datos muestrales x; entonces P(x|θ)puede ser considerada una función de θ. Si es éste el caso, entonces dicha función de probabilidad es la función de máxima verosimilitud de θdado x y se escribe L(θ|x). Podemos por tanto escribir la fórmula de Bayes como P(θ|x) = cL(θ|x)P(θ). Esto quiere decir que la fórmula de Bayes nos asegura que la distribución de probabilidad de θconocidos los datos x es proporcional al producto de la distribución de θa priori y la función de verosimilitud de θdado x: distribución a posteriori ∝verosimilitud ×distribución a priori. La función de verosimilitud L(θ|x)juega un importante papel en la fórmula de Bayes, pues es la función que determina cómo los datos x modifican a θ. Puede ser, por tanto, considerada como una representación de la información sobre θobtenida a partir de los datos observados. La función de verosimilitud está definida salvo constante multiplicativa. Esto es coherente con el lugar que ocupa en la fórmula de Bayes, ya que multiplicar la función de verosimilitud por una constante arbitraria no afecta a la distribución a posteriori de θ. Capítulo 1. Conceptos previos 9 1.2.3. Naturaleza secuencial de la fórmula de Bayes La fórmula de Bayes resulta especialmente útil porque traduce a formulación matemática cómo los conocimientos previos de un suceso o experimento pueden combinarse con nuevos datos. La fórmula permite, en efecto, actualizar de forma continuada la información disponible acerca de los parámetros conforme se van tomando observaciones muestrales. Supongamos que tenemos una muestra inicial de observaciones x1; entonces la fórmula de Bayes se escribe P(θ|x1)∝P(θ)L(θ|x1). Supongamos ahora que disponemos de una segunda muestra de observaciones x2, independientes de la primera muestra. Entonces P(θ|x2, x1)∝P(θ)L(θ|x1)L(θ|x2) ∝P(θ|x1)L(θ|x2). En la expresión precedente, la distribución a posteriori de θdado x1,P(θ|x1), hace de distribución a priori de la segunda muestra. Es obvio que este proceso puede repetirse tantas veces como se quiera. En particular, si tenemos nobservaciones independientes, la distribución a posteriori puede ser recalculada con cada nueva observación, por lo que en el paso mla verosimilitud asociada con las mobservaciones aparece combinada con la distribución a posteriori de θdespués de m−1observacio- nes, dando como resultado la distribución a posteriori siguiente: P(θ|x1, .., xm)∝P(θ|x1, .., xm−1)L(θ|xm), m = 2, ..., n donde P(θ|x1)∝P(θ)L(θ|x1). En consecuencia, la fórmula de Bayes describe el proceso de aprendizaje a partir de la experiencia, y muestra cómo el conocimiento del estado del mundo físico representado por los parámetros está continuamente en cambio según aparecen nuevos datos. 1.3. Probabilidad subjetiva Definición Probabilidad subjetiva: Una probabilidad de que un evento ocurra es un número entre cero y uno que cuantifica la opinión particular (subjetiva) de una persona acerca de cuán probable es que este evento ocurra (o haya ocurrido). 16 2.1. Estimación bayesiana La cantidad γes el nivel creíble del intervalo creíble (tl, tu). Esta definición implica que la variable aleatoria θ|xestá contenida en un γ·100 % intervalo creíble con probabilidad γ. La manera más sencilla de construir intervalos creíbles es la de elegir tlcomo el cuantil (1 −γ)/2ytucomo el cuantil (1 + γ)/2de la distribución a posteriori. Definición intervalos de Máxima Densidad de Probabilidad: Sea γ∈(0,1) un nivel creíble fijo. Un intervalo creíble I= (tl, tu)se llama intervalo de Máxima Densidad de Probabilidad, intervalo M.D.P. (Highest Probability Density en inglés), si verifica f(θ|x)≥f(˜ θ|x) para todo θ∈Iy todo ˜ θ6∈ I. Un intervalo M.D.P. contiene todos los valores del paráetro que tienen mayor densidad de probabilidad a posteriori que aquellos valores del parámetro no contenidos en el intervalo. Bajo ciertas condiciones, la función de densidad a posteriori es igual para los límites del intervalo M.D.P, es decir, f(tl|x) = f(tu|x). Existen contraejemplos en que esto último no se cumple; por ejemplo, cuando la distribución a posteriori es exponencial. El límite menor tlde un intervalo M.D.P. a cualquier nivel será cero puesto que la densidad de una exponencial es monótona decreciente. Para espacios de parámetros discretos Θ, la definición de intervalo creíble debe ser modificada ya que no siempre será posible encontrar un intervalo con el exacto nivel creíble γbuscado. En consecuencia, un γ·100 % intervalo creíble I= (tl, tu)para el parámetro θestará definido como Pθ∈I∩Θf(θ|x)≥γ. Surge ahora la cuestión de si, al igual que en la estimación puntual, los intervalos creíbles son óptimos con respecto a una cierta función de pérdida. Por simplicidad, asumiremos que el parámetro θ∈Θes un escalar con densidad a posteriori f(θ|x). Primero introducimos el concepto de región creíble, una generalización del intervalo creíble. De forma similar puede ser definida una región M.D.P. Definición Región creíble: Un subconjunto C ⊆Θcon RCf(θ|x)dθ =γ. Capítulo 2. Herramientas de la Inferencia Bayesiana 17 se llama γ·100 % región creíble para θcon respecto a f(θ|x). Si C es un intervalo real, entonces C también se llama intervalo creíble. No existe una única γ·100 % región creíble para cada γfijo. Es lógico pensar que la unicidad se obtiene al especificar una cierta función de pérdida £(C, θ). La siguiente función de pérdida tiene en cuenta (para cada γ) la medida de Lebesgue |C|de la región creíble, de forma que regiones pequeñas serán preferibles. Además, la función de pérdida también da resultado desfavorable en el caso de que el parámetro real θno esté contenido en la región creíble. Restringiremos nuestra atención a distribuciones a posteriori continuas. Proposición: Sea f(θ|x)la función de densidad a posteriori de θy sea, para cada γ∈(0,1) fijo, A={C|P(θ∈C|x) = γ} el conjunto de todas las γ·100 % regiones creíbles para θ. Consideremos ahora la función de pérdida £(C, θ) = |C|−IC(θ)para C∈A, θ ∈Θ. Entonces, C es óptimo con respecto a £(C, θ)si y sólo si, para todo θ1∈Cyθ26∈ C, f(θ1|x)≥f(θ2|x), es decir, si y sólo si C es una región M.D.P. Prueba: consideremos un conjunto C∈Acon pérdida esperada RΘ£(C, θ)f(θ|x)dθ =|C|−RΘIC(θ)f(θ|x)dθ =|C|−P(θ∈C|x) = |C|−γ. Para un γfijo, la γ·100 % región creíble C con medida más pequeña |C|minimizará por tanto la pérdida esperada. Sea ahora C∈Averificando f(θ1|x)≥f(θ2|x)para todo θ1∈Cyθ26∈ C, y sea D∈Aotro elemento de A. Tenemos que probar que |C|≤|D|. Sea A˙ ∪Bla unión disjunta de A y B (A˙ ∪B=A∪ByA∩B=∅). Entonces: C=C∩Θ = C∩(D˙ ∪Dc) = (C∩D)˙ ∪(C∩Dc) y análogamente D= (C∩D)˙ ∪(Cc∩D). Queda por tanto probar que |C∩Dc| ≤ |Cc∩D|. Como C es una región M.D.P., se tiene que supCc∩Df(θ|x)≤´ınfC∩Dcf(θ|x), y ya que C, D ∈Ase obtiene 18 2.2. Elección de la distribución a priori P(θ∈C|x) = P(θ∈D|x) = γ. Usando la precedente descomposición de C y D, obtenemos también lo siguiente: P(θ∈C∩D|x) + P(θ∈C∩Dc|x) = P(θ∈C|x) =P(θ∈D|x) = P(θ∈C∩D|x) + P(θ∈Cc∩D|x) y en consecuencia P(θ∈C∩Dc|x) = P(θ∈Cc∩D|x). En definitiva, obtenemos ´ınfC∩Dcf(θ|x)·|C∩Dc| ≤ RC∩Dcf(θ|x)dθ =RCc∩Df(θ|x)dθ ≤supCc∩Df(θ|x)·|Cc∩D| ≤´ınfC∩Dcf(θ|x)·|Cc∩D|. Y por tanto, |C∩Dc|≤|Cc∩D|y se tiene |C|≤|D|. La prueba en la otra dirección es similar. 2.2. Elección de la distribución a priori En la inferencia Bayesiana es fundamental elegir correctamente la distribución a priori, puesto que a través de ella expresamos matemáticamente las ideas previas respecto al modelo. Por tanto, si se dispone de información sobre los datos o el modelo, ésta debe usarse a la hora de construir la distribución a priori. En ciertas situaciones, sin embargo, la elección es algo delicada debido a la ausencia de información previa al modelo. Por otra parte, a menudo es útil restringir el rango de distribuciones a priori a una familia de uno o dos parámetros. 2.2.1. Distribuciones a priori conjugadas Una solución pragmática al problema de elección de la distribución a priori consiste en seleccionar un miembro de una familia de distribuciones tal que la distribución a posteriori pertenezca a la misma familia. Definición distribución a priori conjugada: Sea L(θ) = f(x|θ)una función de verosimilitud basada en la observación X=x. Una clase Gde distribuciones se llama conjugada con respecto a L(θ)si la distribución a posteriori f(θ|x)está en Gpara todo xsiempre que la distribución a priori f(θ)esté en G. La familia G={todas las distribuciones}es trivialmente conjugada con respecto a cualquier función de verosimilitud. En la práctica trataremos de encontrar conjuntos Gmás pequeños específicos para la verosimilitud Lx(θ). Capítulo 2. Herramientas de la Inferencia Bayesiana 19 Notemos que es suficiente estudiar si es conjugado un miembro Xide la muestra aleatoria X1:nEn efecto, si la distribución a priori es conjugada, la distribución a posteriori después de haber observado la primera observación es, por definición, del mismo tipo y sirve como la nueva distribución a priori para la próxima observación. La nueva distribución a posteriori, incorporando ahora la segunda observación, pertenece de nuevo a la clase de conjugación. Este proceso de asimilación de datos continúa secuencialmente y sólo cambia el valor de los parámetros de la distribución. 2.2.2. Distribuciones a priori impropias La distribución a priori tiene influencia sobre la distribución a posteriori. Si se desea minimizar dicha influencia, entonces se puede especificar una distribución a priori imprecisa, por ejemplo una con varianza muy grande. En un caso límite, esto puede conducir a una distribución a priori impropia, con función de densidad cuya integral no es igual a la unidad. Usando estas distribuciones siempre es necesario comprobar que al menos la distribución a posteriori es propia. Definición distribución a priori impropia: Una distribución a priori con función de densidad f(θ)≥0se llama impropia si RΘf(θ)dθ =∞Pθ∈Θf(θ) = ∞ para parámetros continuos y discrretos, respectivamente. 2.2.3. Distribuciones a priori no informativas: Regla de Jeffreys Las distribuciones a priori no informativas son útiles cuando deseamos que la inferencia no se vea afectada por información que no provenga de los datos presentes. Estas distribuciones son también apropiadas cuando tenemos muy poco conocimiento previo en comparación con la información contenida en los nuevos datos (en estos casos, es de esperar que la verosimilitud domine a la distribución a priori). Para establecer una distribución a priori no informativa, la opción más obvia es escoger una distribución a priori localmente uniforme tal que fθ(θ)∝1, en cuyo caso la distribución a posteriori es proporcional a la función de verosimilitud. Si el espacio de parámetros no está acotado, corremos el riesgo de que la distribución localmente uniforme sea impropia. Ahora bien, surgen problemas en este planteamiento. Supongamos que h(θ)es una transformación inyectiva y diferenciable del parámetro θ, que tiene distribución a priori localmente uniforme con densidad fθ(θ)∝1. Usando un cambio de variables, obtenemos la correspondiente distribución a priori para φ 20 2.2. Elección de la distribución a priori fφ(φ) = fθ{h−1(φ)}·|dh−1(φ) dφ | ∝ |dh−1(φ) dφ |. Este término no es necesariamente constante. En efecto, fφ(φ)será independiente de φúnicamente si hes lineal. Si hno lo es, la función de densidad a priori fφ(φ)dependerá de θy no será, por tanto, localmente uniforme. Aun así, si hemos elegido una parametrización en función de φdesde el comienzo, podemos también haber elegido una distribución a priori localmente uniforme fφ(φ)∝1. Que no exista invarianza bajo reparametrización es un inconveniente de este tipo de distribuciones a priori. Hagamos notar que asumimos implícitamente que podemos aplicar el cambio de variables a funciones de densidad impropias de la misma manera que a las propias. A pesar de todo lo dicho, existe un caso particular de elección de la distribución a priori que es invariante bajo reparametrización: se trata de la distribución a priori de Jeffreys. Definición distribución a priori de Jeffreys: Sea X una variable aleatoria con función de verosimilitud f(x|θ)donde θes un parámetro escalar desconocido. Entonces la distribución a priori de Jeffreys está definida como f(θ)∝pJ(θ) donde J(θ)es la información de Fisher esperada de θ. La ecuación enunciada también se conoce como Regla de Jeffreys. La distribución de Jeffreys es proporcional a la raíz cuadrada de la información de Fisher esperada, lo cual puede proporcionar una distribución a priori impropia. A primera vista resulta sorprendente que esta elección sea invariante bajo reparametrización, pero el siguiente resultado nos confirma que, en efecto, así es. Proposición: Invarianza de la distribución a priori de Jeffreys La distribución a priori de Jeffreys es invariante bajo parametrizaciones inyectivas de θ: Si fθ(θ)∝pJθ(θ) entonces la función de densidad φ=h(θ)es fφ(φ)∝pJφ(φ) Para demostrarlo, usaremos el cambio de variables y también el siguiente resultado: Información de Fisher esperada de una transformación Capítulo 2. Herramientas de la Inferencia Bayesiana 21 Sea Jθ(θ)la información de Fisher esperada de un parámetro escalar θ, y φ= h(θ)una transformación inyectiva de θ. La información de Fisher esperada Jφ(φ)de φ puede entonces calcularse como sigue: Jφ(φ) = Jθ(θ){dh−1(φ) dφ }2=Jθ(θ){dh(θ) dθ }−2. Con esto ya podemos probar el resultado de invarianza de la distribución a priori de Jeffreys. Prueba: De f(θ)∝pJ(θ)se deduce, usando el cambio de variables y el resultado precedente, que fφ(φ) = fθ(θ)·|dh−1(φ) dφ | ∝pJθ(θ)·|dh−1(φ) dφ |=qJθ(θ)·{dh−1(φ) dφ }2=pJφ(φ) y con esto queda probada la proposición. Si hubiésemos elegido la parametrización φ=h(θ)desde el principio, al aplicar la regla de Jeffreys habríamos obtenido la misma distribución a priori: fφ(φ)∝pJφ(φ). 2.3. Distribución a posteriori predictiva Consideremos una muestra aleatoria X1:nde una distribución con función de densidad f(x|θ). Queremos predecir una observación (independiente) que se produzca en el futuro, Y=Xn+1 de f(x|θ). La distribución a posteriori predictiva de Y|x1:n, que puede calcularse fácilmente usando leyes de probabilidad, es la siguiente: f(y|x1:n) = Rf(y, θ|x1:n)dθ =Rf(y|θ, x1:n)f(θ|x1:n)dθ =Rf(y|θ)f(θ|x1:n)dθ. Observemos que la última línea se deduce de la hipótesis de independencia condicional de Y y X1:n, dado θ. Por tanto, para calcular f(y|x1:n), basta con integrar el producto de la verosimilitud f(y|θ)y la función de densidad a posteriori f(θ|x1:n) con respecto de θ. El resultado f(y|x1:n)se llama distibución a posteriori predictiva, en contraposición a la distribución a priori predictiva: f(y) = Rf(y|θ)f(θ)dθ. donde f(y|x1:n)se ha sustituido por la densidad a priori f(θ). Conviene señalar que la distribución a posteriori predictiva f(y|x1:n)no es igual a la distribución f(y|ˆ θ). De hecho, f(y|ˆ θ)se obtiene al sustituir la distribución a posteriori f(y|x1:n)por una medida puntual en ˆ θ, normalmente el estimador de máxima 22 2.3. Distribución a posteriori predictiva verosimilitud, en la expresión para la distribución predictiva. Por su parte, la predicción Bayesiana automáticamente incorpora la falta de certeza en la estimación del parámetro. Cálculo de la distribución a posteriori predictiva Ya que la distribución a priori predictiva es sencillamente el denominador en el teorema de Bayes, puede calcularse si se conocen la verosimilitud y las distribuciones a priori y a posteriori: f(y) = f(y|θ)f(θ) f(θ|y) lo cual está bien definido para todo θ. Una fórmula similar también estará entonces bien definida para la distribución a posteriori predictiva: f(y|x) = f(y|x,θ)f(θ|x) f(θ|x,y)=f(y|θ)f(θ|x) f(θ|x,y) donde la última igualdad es consecuencia de la independencia condicionada de X y Y. Si f(θ)es conjugada con respecto a f(x|θ), entonces f(θ|x)yf(θ|x, y)pertenecen a la misma familia de distribuciones, y la distribución a posteriori predictiva puede obtenerse sin integrar explícitamente. Capítulo 3 Inferencia bayesiana sobre una proporción En este capítulo se estudiarán métodos Bayesianos para inferir sobre el valor desconocido pde una distribución Bernoulli asumiendo conocida la distribución a priori de py habiendo observado datos de un proceso que siga esta distribución. 3.1. Cálculo de la distribución a posteriori Verermos ahora distintos modos de establecer la distribución a priori, tanto para caso continuo como para discreto, y obtendremos como consecuencia distribuciones a posteriori diferentes en cada caso. 3.1.1. Distribución a priori discreta para ˜p Enunciamos de nuevo la fórmula de Bayes para mvalores de una distribución Bernoulli de parámetro p. Si la probabilidad a priori de que la distribución tome el valor pies P(˜p=pi), donde i= 1,2, ..., m y se han observado unos datos D, la fórmula de Bayes queda como sigue: P(˜p=pi|D) = P(˜p=pi)L(pi|D) P(D)= ( 1 P(D))P(˜p=pi)L(pi|D) donde 1 P(D)no depende de p; podemos por tanto llamarlo ky tenemos: P(˜p=pi|D) = kP(˜p=pi)L(pi|D) y por proporcionalidad, deducimos, como hicimos en el capítulo anterior, que a posteriori ∝verosimilitud ×distribución a priori. 23 24 3.1. Cálculo de la distribución a posteriori Función verosimilitud de la distribución Bernoulli Si se observan nexperimentos Bernoulli (distribución binomial), los datos observados pueden verse como una secuencia de éxitos y fracasos en un orden específico. Si hay réxitos y n−rfracasos, entonces L(pi|D) = pr iqn−r i. Si el orden exacto de éxitos y fracasos no queda registrado en los datos entonces L(pi|D) = n rpr iqn−r i aunque esta verosimilitud es equivalente a la anterior puesto que n rno depende de p. Caso distribución uniforme Un caso importante es aquél en que la distribución a priori de ˜pes discreta y uniforme. Usando el resultado de proporcionalidad expuesto más arriba, tenemos lo siguiente: a posteriori ∝verosimilitud donde las probabilidades uniformes a priori han sido englobadas por la constante de proporcionalidad. Para obtener la igualdad basta encontrar una constante normalizadora para la verosimilitud. 3.1.2. Distribución a priori continua para ˜p En la sección anterior, los valores del parámetro ˜pde la distribución Bernoulli estaban restringidos a un conjunto finito de valores; de esta forma se dedujo la relación de proporcionalidad entre distribución a priori y a posteriori. Esta idea se puede extender para el caso más general en el que ppuede tomar cualquier valor entre 0 y 1. Dado que el intervalo entre 0 y 1 es continuo, es necesario considerar funciones de densidad para el valor desconocido de p, tanto para la distribución a priori como para la distribución a posteriori. 3.1.2.1. Distribución continua a priori para ˜p Sea p∗un valor fijo de py sea dp la longitud de un intervalo pequeño contenido en los valores de p. Consideremos también un valor rque denote el número de éxitos en nexperimentos Bernoulli. Para la función de densidad a priori f(p)se tiene: f(p∗)dp ≈P(p∗−dp 2<˜p<p∗+dp 2). Capítulo 3. Inferencia bayesiana sobre una proporción 25 Es decir, f(p∗)dp es la probabilidad a priori aproximada de que ˜ptome un valor contenido en el intervalo de centro p∗y longitud 2∗dp. Por la fórmula de Bayes tenemos: P(p∗−dp 2<˜p<p∗+dp 2|r)≈P(p∗−dp 2<˜p<p∗+dp 2)L(p∗|r,n) P(r,n). Si dejamos que P(r|n)sea parte de la constante de proporcionalidad, tenemos, de manera aproximada, lo siguiente: P(p∗−dp 2<˜p<p∗+dp 2|r)∝P(p∗−dp 2<˜p<p∗+dp 2)L(p∗|r, n). El miembro izquierdo de la ecuación es aproximadamente la probabilidad a posteriori f(p∗|r, n)dp por tanto se obtiene: f(p∗|r, n)dp ≈f(p∗)dpL(p∗|r,n) P(r|n). Cuando dp tiende a cero, esta aproximación tiende a la igualdad y se pueden cancelar los dp a ambos lados para obtener la fómula de Bayes para el caso continuo: f(p∗|r, n) = f(p∗)L(p∗|r,n) P(r|n). y usando proporcionalidad f(p∗|r, n)∝f(p∗)L(p∗|r, n). 3.1.2.2. Distribución continua uniforme para ˜p Al igual que en el caso discreto, cobra importancia el caso uniforme para la distribución continua a priori de ˜p. Distribución a priori uniforme para ˜pquiere decir que la función de densidad de probabilidad f(p)es constante para todo p, por tanto debe ser un rectángulo con base en el intervalo del cero al uno; además, como debe tener área igual a uno, el rectángulo debe ser un cuadrado, y f(p)debe ser igual a uno para todos los valores de pentre cero y uno. Usando la ecuación de proporcionalidad del apartado anterior tenemos, dado r,: f(p|r, n)∝1×L(p|r, n). Cuando estudiamos antes el caso uniforme para distribución discreta a priori, vimos que la función de verosimilitud era una función discreta en un conjunto discreto de valores de p, y en un pparticular tomaba un valor proporcional a prqn−r. La verosimilitud es ahora la misma salvo por ser esta vez una función continua definida en todo el intervalo de cero a uno: 32 4.2. Cálculo de la distribución a posteriori t0 1=t1, t0 2=t2−t1, t0 3=t3−t2, ..., t0 r=tr−tr−1, entonces las variables aleatorias t0 1, ..., t0 rson independientes y están distribuidas según una exponencial de parámetro λ. Por tanto, ya que f(t1, ..., tr|λ, r) = f(t0 1, ..., t0 r|λ, r) y los t0 ison independientes, tenemos: f(t0 1, ..., t0 r|λ, r) = f(t0 1|λ)×... ×f(t0 r|λ). Puesto que cada una de estas densidades es de la forma λe−λt0 i, su producto es como sigue: λe−λt0 1×... ×λe−λt0 r=λre−λPt0 i. Y es trivial que Pr 1t0 i=tr. En conclusión, tenemos: L(λ|t1, ..., tr) = f(t1, ..., tr|λ) = λre−λtr La verosimilitud depende sólo de tr; deducimos por tanto que el estadístico tres suficiente. Verosimilitud de la distribución de Poisson Para una muestra de Poisson, con Tfijo, la verosimilitud de λdados los datos t1, ..., tryTes una modificación de la verosimilitud obtenida para el caso de la distribución Poisson inversa. Ahora debemos establecer la verosimilitud de forma que no ocurran más sucesos desde trhasta T. La probabilidad de cero sucesos entre tryTes la probabilidad de observar cero sucesos en un intervalo de tiempo de longitud (T−tr): Pp(˜r= 0|λ×(T−tr)) = e−λ(T−tr)(λ(T−tr))0 0! =e−λ(T−tr). Multiplicando la verosimilitud de la distribución inversa por esta probabilidad obtenemos la verosimilitud para el caso de la distribución de Poisson: L(λ|t1, ..., tr, T) = λre−λtr×e−λ(T−tr) =λre−λ(tr+T−tr)=λre−λT . El estadístico suficiente es ahora r: basta conocer rpara deducir la verosimilitud. 4.2. Cálculo de la distribución a posteriori En este apartado estudiaremos cómo obtener la distribución a posteriori de una distribución de Poisson, partiendo de una distribución a priori no informativa y también estudiando el caso conjugado. Capítulo 4. Distribución de Poisson 33 4.2.1. Distribución a priori no informativa La verosimilitud de la distribución de Poisson es de la forma L(λ|t1, ..., tr)∝λre−λT . Por tanto la log-verosimilitud es: log(L(λ|t1, ..., tr)) = constante +rlog λ −λT. Ahora deduciremos la expresión de una distribución a priori no informativa mediante la regla de Jeffreys. Primero calculamos derivadas sucesivas respecto al parámetro λ: ∂log(L(λ|t1,...,tr)) ∂λ =r λ−T, ∂2log(L(λ|t1,...,tr)) ∂λ2=−r λ2. Notemos que r=Pti, luego r/T =¯ t. Para ¯ tdistinto de cero, el estimador de máxima verosimilitud de λes ˆ λ=¯ t, y tenemos J(ˆ λ) = (−1 n ∂2log(L(λ|t1,...,tr)) ∂λ2)ˆ λ=1 ˆ λ. La regla de Jeffreys nos dice entonces que la distribución a priori no informativa para ˜ λes f(˜ λ)∝qJ(˜ λ)∝λ−1/2. Como consecuencia, la función de densidad a posteriori cumple la siguiente relación de proporcionalidad: f(˜ λ|t1, ..., tr)∝˜ λr−1 2e−λT e integrando se obtiene la constante normalizadora, que es la siguiente: constante =T−(r+1 2)[Γ(r+1 2)]−1. 4.2.2. Distribución Gamma a priori para Poisson ˜ λ De manera similar a la distribución Beta a priori para una Bernoulli, estudiamos ahora el caso conjugado para la distribución de Poisson. La familia paramétrica que estudiaremos es la conocida como familia Gamma de distribuciones, llamada así por su relación con la función Gamma, definida como sigue: Γ(a) = R+∞ 0xa−1e−xdx. Por su parte, la función de densidad de probabilidad que usaremos es de la siguiente forma: 34 4.2. Cálculo de la distribución a posteriori f(x) = xa−1e−x/b Γ(a)ba. La familia Gamma depende de los parámetros ayby, al ser función de densidad, su integral en todo el intervalo real es uno, o expresado de otro modo: R+∞ 0xa−1e−xbdx = Γ(a)ba. Vemos ahora la relación de esta función de densidad con la función Gamma (para b= 1 es, de hecho, la misma expresión integral). Para adaptar estas consideraciones a la distribución de Poisson, reemplazamos la xpor ˜ λ, y los parámetros aybpor r0 y1/T0respectivamente. Obtenemos, como consecuencia, la densidad a priori para el parámetro λde una distribución de Poisson: fγ(˜ λ|r0, T0) = ˜ λr0−1e−˜ λT 0 Γ(r0)(1/T0)r0 para ˜ λ > 0yr0, T0>0. Como eventualmente estudiaremos relaciones de proporcionalidad, las constante pueden obviarse hasta calcular las constante normalizadora de la distribución a posteriori. Los valores r0yT0de la distribución a priori no tienen por qué ser iguales a los de la función de verosimilitud, de aquí que estén distinguidos con la prima. Escribamos la media y la varianza de la distribución Gamma con respecto a los parámetros de una Poisson. Media: E(˜ λ|r0, T0) = r0 T0. Varianza V(˜ λ|r0, T0) = r0 T02. Para finalizar, apliquemos el teorema de Bayes y obtengamos la distribución a posteriori para el parámetro de la distribución de Poisson. Por una parte, la distribución a priori es proporcional a ˜ λr0−1e−˜ λT0 y por otra, la verosimilitud de ˜ λdados los datos (resumidos en los valores ryT) es proporcional a ˜ λre−˜ λT . En consecuencia, la distribución a posteriori es proporcional a ˜ λr0−1e−˜ λT0ט λre−˜ λT . Definimos ahora r00 =r0+ryT00 =T0+Ty obtenemos la siguiente expresión para la distribución a posteriori: ˜ λr00 −1e−˜ λT00 . Capítulo 4. Distribución de Poisson 35 Esta expresión se corresponde con una función de densidad Gamma de parámetros r00 yT00. La constante normalizadora se deduce de la expresión para la función de densidad Gamma que se dio anteriormente, y por tanto la densidad a posteriori (conjugada) que buscamos es f(˜ λ|r00, T00) = ˜ λr00 −1e−˜ λT 00 Γr00 (1/T00 )r00 . 4.3. Distribución predictiva Para ˜ λdesconocido, las predicciones acerca de (˜r|T, distribuc´on a priori)vienen expresadas por la distribución predictiva. Distribución a priori discreta para ˜ λ P(r|T, dist. a priori) = PiPP0(r|λiT)P(λi), r = 0,1, ... Distribución a priori continua para ˜ λ P(r|T, dist. a priori) = R+∞ 0PP0(r|λT)f(λ)dλ, r = 0,1, ... Distribución a priori Gamma (r0, T0)para ˜ λ P(r|T, r0, T0) = R+∞ 0PP0(r|λT)fγ(λ|r0, T0)dλ, r = 0,1, ... que corresponde a la función de distribución de probabilidad de una Gamma-Poisson. Escribamos su media y varianza: E(˜r|T, r0, T0) = ( r0 T0)T V(˜r|T, r0, T)=(r0 T0)T(T0+T T0) 4.4. Contraste de hipótesis Consideremos un test de hipótesis del tipo H0:λ≤λ0 H1:λ > λ0 donde se observan rsucesos en Texperimentos. El p-valor en sentido clásico es PP0(˜r≥r|λ0×T). En el análisis Bayesiano, con una distribución a priori no informativa, se tiene a posteriori para ˜ λuna distribución Gamma con parámetros r00 =r+ 1 yT00 =T. Se calcula lo siguiente: 36 4.4. Contraste de hipótesis P(H0cierta|datos) = Pγ(˜ λ≤λ0|r00 =r+ 1, T00 =T). Cuando los parámetros de la gamma r00 yT00 son enteros, la distribución Gamma y la distribución de Poisson son equivalentes, en el sentido siguiente: Pγ(˜ λ≤λ0|r00 =r+ 1, T00 =T)≡PP0(˜r≥r00|λ0×T00). Y para λ0×Tgrande, se tiene que PP0(˜r≥r|λ0×T) no es muy diferente a PP0(˜r≥r00 =r+ 1|λ0×T) y el p-valor clásico es en gran medida similar a la probabilidad a posteriori de que ˜ λ≤λ. Capítulo 5 Distribución normal 5.1. Función de verosimilitud para (µ, σ2) Sean x1, x2, ..., xnuna muestra de nobservaciones (independientes) de un proceso que sigue una distribución normal, es decir, cada xies un valor observado de ˜xique se distribuye según una normal N(µ, σ2)con función de densidad como sigue: fN(xi|µ, σ2) = 1 σ√2πe−1 2(x−µ σ)2 Ya hemos visto que la función de la verosimilitud es de la fórma L(par0ametros|datos) = f(datos|par0ametros). Ahora tenemos dos parámetros, por lo que la verosimilitud queda L(µ, σ2|x1, x2, ..., xn) = f(x1, x2, ..., xn|n, µ, σ2) donde fdenota la función de densidad conjunta. Como los ˜xison independientes, f(x1, x2, ..., xn|µ, σ2) = fN(x1|µ, σ2)×fN(x2|µ, σ2)×... ×fN(xn|µ, σ2). Reemplazando cada función de densidad por su expresión llegamos a L(µ, σ2|x1, x2, ..., xn) = 1 σ√2πe−1 2(x1−µ σ)2×1 σ√2πe−1 2(x2−µ σ)2×... ×1 σ√2πe−1 2(xn−µ σ)2. Agrupando términos tenemos L(µ, σ2|x1, x2, ..., xn)=(σ√2π)−ne−(1/2σ2)[(x1−µ)2+(x2−µ)2+...+(xn−µ)2]. O de manera equivalente: L(µ, σ2|x1, x2, ..., xn)=(σ√2π)−ne−(1/2σ2)Σn 1(xi−µ)2. Podemos simplificar esta expresión sustrayendo (√2π)−ny cambiando el igual por una relación de proporcionalidad: 37 38 5.1. Función de verosimilitud para (µ, σ2) L(µ, σ2|x1, x2, ..., xn)∝(σ)−n×e−(1/2σ2)Σn 1(xi−µ)2. Los dos términos que quedan no pueden quitarse porque dependen de los parámetros. Ahora haremos algunas transformaciones en lo obtenido para llegar a una expresión que será de utilidad más adelante. Sea ¯xla media muestral de las observaciones x1, x2, ..., xn. Consideremos el término Pn 1(xi−µ)2 que podemos escribir como Pn 1(xi−¯x+ ¯x−µ)2. Desarrollando obtenemos Pn 1(xi−¯x+ ¯x−µ)2=Pn 1[(xi−¯x) + (¯x−µ)]2= Pn 1(xi−¯x)2+ 2 Pn 1(xi−¯x)(¯x−µ) + Pn 1(¯x−µ)2. Estudiemos los tres término separadamente: El primer término está relacionado con la varianza muestral, denotada s2y definida como s2=Pn 1(xi−¯x¯ )2 n. Por tanto Pn 1(xi−¯x)2=ns2. Para el segundo término tenemos 2Pn 1(xi−¯x)(¯x−µ) = 2(¯x−µ)Pn 1(xi−¯x). Pero Pn 1(xi−¯x) = 0, luego el segundo término se anula. En cuanto al tercero, se observa que el sumando es igual para todo n, luego Pn 1(¯x−µ)2=n(¯x−µ)2. Uniendo estos resultados e implementándolos en la función de verosimilitud se obtiene: L(µ, σ2|x1, x2, ..., xn)∝(σ)−ne−(1/2σ2)[ns2+n(¯x−µ)2]. Capítulo 5. Distribución normal 39 Lo que conduce a la expresión para la verosimilitud: L(µ, σ2|x1, x2, ..., xn)∝(σ)−ne−ns2/2σ2×e−n(µ−¯x)2/2σ2. Obsevación: antes de proseguir con nuestro estudio de la distribución Normal, conviene aclarar cierto detalle que estará presente a lo largo del capítulo. Dependiendo de lo que resulte más convenente para los cálculos, a veces manejaremos expresiones en función de la varianza, σ2, y en otras ocasiones en función de la desviación estándar, σ. Ahora bien, es muy sencillo pasar de tener una expresión en función de la primera a tenerla en función de la segunda, y viceversa. Pongamos como ejemplo que tenemos la función de densidad a priori para σ, pero nos interesa tener dicha densidad en función de σ2. Consideremos el siguiente cambio de variable: h:σ7→ σ2 cuya función inversa es h−1:σ27→ √σ2. Si ahora calculamos el Jacobiano de h−1 J(h−1) = |dh−1 dσ2|=1 2σ2 tenemos que la distribución a priori de σ2no es más que f(σ2) = f(σ)J(h−1) = f(σ)1 2σ2. 5.2. Regla de Jeffreys Aplicaremos ahora la regla de Jeffreys para obtener distribuciones a piori no informativas para una normal de parámetros µyσ2. Primero, escribamos lo siguiente: f(x|µ, σ2)∝σ−1h(x−µ σ) donde xno depende de los parámetros y htiene ciertas condiciones de regularidad. Supongamos también que tenemos una muestra de nobservaciones x1, ..., xnde la distribución. 40 5.2. Regla de Jeffreys 5.2.1. Caso ˜µdesconocido y σconocido La verosimilitud es L(˜µ|σ, x)∝Qn i=1 h(xi−˜µ σ). Ahora tenemos ∂log f(x|˜µ,σ) ∂˜µ=−1 σ h0 x h(x), donde y= (x−˜µ σ). Y ahora calculamos la información de Fisher esperada: J(˜µ) = Ex|˜µ[∂log f(x|˜µ,σ) ∂˜µ]2=1 σ2Ey[h0 x h(x)]2 donde la esperanza en el último término se toma sobre la distribución f(y) = h(y). Ya que la expresión para la esperanza no incluye a ˜µyσ2es conocido, concluimos que J(˜µ) = constante. Esto conduce a tomar a ˜µlocalmente uniforme a priori: f(˜µ|σ)∝J(˜µ)1/2=constante y en consecuencia la correspondiente distribución a posteriori para ˜µes f(˜µ|σ, x)∝Qn i=1 h(xi−˜µ σ),−∞ <˜µ < ∞. Nota: conviene señalar que la distribución a priori obtenida para la media es impropia, pues la integral de su función de densidad, R∞ −∞ constanted˜µ, no es finita. 5.2.2. Caso µconocido, ˜σdesconocido La verosimilitud es ahora: L(˜σ|µ, x)∝σ−nQn i=1 h(xi−µ ˜σ). Ya que ∂log f(x|µ,˜σ) ∂˜σ=−1 ˜σ[1 + yh0 x h(x)], donde y= (x−µ ˜σ). se deduce que J(˜σ) = 1 ˜σ2Ex[1 + yh0 x h(x)]2∝1 σ2. En consecuencia, la regla de Jeffreys nos lleva a la distribución a priori siguiente: f(˜σ|µ)∝1 ˜σóf(log ˜σ)∝constante. y por tanto la correspondiente distribución a posteriori para ˜σes f(˜σ|µ, x)∝σ−(n+1) Qn i=1 h(xi−˜µ σ), σ > 0. De nuevo se ha obtenido una distribución a priori impropia (pues la integral de 1/˜σ no es finita). Capítulo 5. Distribución normal 41 5.3. Cálculo de la distribución a posteriori 5.3.1. Caso ˜µdesconocido y σ2conocido Sean x1, x2, ..., xnuna muestra de nobservaciones de una normal N(˜µ, σ2).Asumimos ahora que la varianza σ2es conocida. Haremos inferencia sobre la media desconocida ˜µ. 5.3.1.1. Distribución a priori no informativa En cuanto a la función de verosimilitud, todos los factores, salvo aquellos que dependen de ˜µ, pueden ser absorbidos por la constante de proporcionalidad; por tanto, la función de verosimilitud se reduce a L(˜µ|n, ¯x, σ2)∝e−n(˜µ−¯x)2/2σ2. Los x1, x2, ..., xnhan sido sustituidos por ny¯xya que estos dos últimos términos bastan para determinar la verosimilitud. Como distribución a priori no informativa para ˜µusaremos la obtenida mediante la regla de Jeffreys, es decir: f(˜µ|σ)∝c. La distribución a posteriori para ˜µverifica: f(˜µ|¯x)∝f(˜µ)L(˜µ|¯x)∝f(˜µ)f(¯x|˜µ, σ2), donde f(˜µ)es la distribución a priori. Por tanto tenemos: f(˜µ|n, ¯x, σ2)∝e−n(˜µ−¯x)2/2σ2. La constante normalizadora que hace que la integral del miembro a la derecha sea uno es c−1=R∞ −∞ e−n(˜µ−¯x)2/2σ2d˜µ= (2πσ2 n)1/2 y, en consecuencia, la función de densidad a posteriori es f(˜µ|n, ¯x, σ2)=(2πσ2 n)−1/2e−n(˜µ−¯x)2/2σ2,−∞ <˜µ < ∞. Concluimos que la distribución que sigue ˜µ(a posteriori y de manera aproximada) es una normal N(¯x, σ2/n). 48 5.3. Cálculo de la distribución a posteriori f(¯x|˜µ, ˜σ2) = √n √2π˜σexp[−n 2˜σ2(¯x−˜µ)2] f(s2|˜σ2) = 1 Γ[1 2(n−1)] (n−1 2˜σ2)1 2(n−1)(s2)1 2(n−1)−1×exp[−(n−1)s2 2˜σ2], s > 0. Por tanto, dados los datos x, la distribución conjunta a posteriori de (˜µ, ˜σ)es: f(˜µ, ˜σ|x)∝f(˜µ, ˜σ)f(¯x|˜µ, ˜σ2)f(s2|˜σ) donde f(˜µ, ˜σ)es la distribución a priori. Es conveniente reparar ahora en el hecho de que no es inadecuado considerar un parámetro de situación independientemente distribuido de un parámetro de escala. En efecto, cualquier idea preconcebida que se tenga sobre el parámetro de situación de una distribución no se verá influenciada, de forma general, por la idea que se tenga acerca del valor del parámetro de escala. Tomando nuestro caso de distribución normal, donde µyσson parámetros de situación y escala respectivamente, tenemos que f(˜µ|˜σ). = f(˜µ)(donde hemos asumido que a priori ˜µy˜σson aproximadamente independientes) y, por tanto, concluimos que f(˜µ, ˜σ) = f(˜µ|˜σ)f(˜σ). =f(˜µ)f(˜σ). Adoptamos como distribuciones a priori no informativas para los parámetros las siguientes: f(˜µ)∝c f(log˜σ)∝c, f(˜σ)∝˜ σ−1. Usando lo anterior llegamos a que la distribución a posteiori es: f(˜µ, ˜σ|x) = k˜σ−(n+1)exp−1 2˜σ2[νs2+n(˜µ−¯x)2],−∞ <˜µ < ∞, σ > 0 donde ν=n−1ykes la constante de normalización. Usaremos las siguientes fórmulas de integración: R∞ 0x−(p+1)e−ax−2dx =1 2a−p 2Γ(p/2), a > 0, p > 0 R∞ ∞exp[−1 2(x−η c)2]dx =√2πc, −∞ < η < ∞, c > 0 para llegar a k=pn 2π[1 2Γ(ν 2)]−1(νs2 2)ν/2. En conclusión: f(˜µ, ˜σ|x) = pn 2π[1 2Γ(ν 2)]−1(νs2 2)ν/2˜σ−(n+1)exp−1 2σσ2[νs2+n(˜µ−¯x)2], −∞ <˜µ < ∞,˜σ > 0. Capítulo 5. Distribución normal 49 La distribución conjunta a posteriori de (˜µ, ˜σ)obtenida puede escribirse como f(˜µ, ˜σ|x) = f(˜µ|˜σ, x)f(˜σ|x) donde el primer factor es la distribución condicionada a posteriori de ˜µdado ˜σ, y el segundo factor es la distribución marginal a posteriori de ˜σ. Estudiemos ahora las distribuciones componentes de f(˜µ, ˜σ|x). Distribución de ˜µcondicionada a ˜σ Si tratamos a ˜σcomo una constante conocida, la distribución condicional a posteriori de ˜µdado ˜σes f(˜µ|˜σ, x) = √n √2π˜σexp[−n 2˜σ2(˜µ−¯x)2],−∞ <˜µ < ∞. Es decir, dado ˜σ, ˜µestá distribuida a posteriori según una normal N(¯x, ˜σ2/n. Distribución marginal de ˜σ La distribución marginal a posteriori de ˜σes f(˜σ|¯x) = f(˜µ,˜σ|¯x) f(˜µ|˜σ,¯x). De la distribución a posteriori y la distribución de ˜µcondicionada a ˜σllegamos a f(˜σ|x) = k0˜σ−(ν+1)exp(−νs2 2˜σ2),˜σ > 0. donde k0= [1 2Γ(ν 2)]−1(νs2 2)ν/2. Dicho de otro modo, ˜σestá distribuida a posteriori como √νsχ−1 ν. Esta distribución es la misma que obtuvimos cuando la media era conocida y la varianza desconocida, pero con la diferencia de que cuando la media es conocida la distribución de ˜σ/(√ns), con ns2=Pn i=1(xi−µ)2, tiene ngrados de libertad, mientras que cuando la media es desconocida, la distribución de ˜σ/(√νs), con νs2=Pn i=1(xi−¯x)2, tiene ν=n−1 grados de libertad. Las inferencias sobre ˜σ, con ˜µdesconocido, pueden realizarse tal como hicimos en el caso µconocido y ˜σdesconocido. Distribución marginal a posteriori de ˜µ La distribución marginal a posteriori de ˜µpuede ser obtenida integrando con respecto a ˜σen la distribución conjunta a posteriori de ˜µy˜σ: f(˜µ|x) = R∞ 0f(˜µ, σ|x)d˜σ=R∞ 0f(˜µ|˜σ, x)f(˜σ|x)d˜σ. 50 5.3. Cálculo de la distribución a posteriori Usaremos la siguiente fórmula integral R∞ 0x−(p+1)e−ax−2dx =1 2a−p 2Γ(p/2), a > 0, p > 0 para llegar a f(˜µ|x) = (s/√n)−1 B(1 2ν, 1 2)√ν[1 + n(˜µ−¯x)2 νs2]−1 2(ν+1),−∞ <˜µ < ∞ donde B(p, q)es la función beta B(p, q) = Γ(p)Γ(q)/Γ(p+q). Si escribimos t=˜µ−¯x s/√n entonces tenemos f(t=˜µ−¯x s/√n/x) = 1 B(1 2ν, 1 2)√ν(1 + t2 ν)−1 2(ν+1),−∞ <t<∞ que será una distribución tde Student con ν=n−1grados de libertad. Conviene recordar que ˜µes la variable aleatoria en la expresión y (¯x, s2)son valores muestrales conocidos. Notaremos a la distribución obtenida como t(¯x, s2/n, ν). Agrupamos ahora los resultados y resumimos lo que hemos obtenido: Consideremos los valores muestrales ¯xys2distribuidos independientemente según distribuciones N(˜µ, ˜σ2/n)y(˜σ2/ν)χ2 νrespectivamente. Supongamos que a priori ˜µylog ˜σestán distribuidas (aproximadamente) de manera independiente y localmente uniforme. Entonces, dados (¯x, s2), se tiene: ˜σsigue la distribución √νsχ−1 ν. ˜µcondicionada a ˜σestá distribuida según una normal N(¯x, σ2/n). ˜µtiene distribución marginal t(¯x, s2/n, ν). Intervalos a posteriori para ˜µ La distribución a posteriori obtenida para ˜µes una distribución simétrica centrada en ¯xcon factor de escala s/√n. Para obtener intervalos M.D.P. para la media debemos usar una tabla de valores de la distribución tde Student. Si denotamos t(α/2)(τ)como el valor para el cual P{t > t(α/2)(τ)}=α/2 Capítulo 5. Distribución normal 51 entonces, por la simetría de la distribución tde Student, P{t < −t(α/2)(τ)}=α/2 y por tanto P{|t|< t(α/2)(τ)}= 1 −α/2. Se deduce que los límites de un (1 −α)intervalo M.D.P. de t(¯x, s2/n, τ)vienen dados por ¯x±(s √n)t(α/2)(τ). 5.3.3.2. Distribución conjugada Cuando existe conocimiento a priori sobre los parámetros desconocidos de una Normal, seguramente queramos incorporarlos a nuestro análisis Bayesiano. Usualmente se asume independencia a priori entre los parámetros, y se toma como distribución a priori el producto de una Normal sobre la media y una Gamma inversa para la varianza. En este caso, ninguna de las distribuciones marginales a posteriori pertenece a la familia paramétrica, y la densidad conjunta a posteriori claramente no será producto de una Normal y una Gamma inversa, por tanto, no estamos afrontando un análisis conjugado propiamente dicho. Para designar esta situación se habla de análisis semiconjugado, ya que la distribución a priori sería conjugada si uno de los parámetros fuera conocido. Sin embargo, sí es posible estalecer una distribución a priori conjugada cuando los parámetros son desconocidos. La densidad conjunta a priori, que llamarenos densidad Normal-Gamma inversa, es el producto de la marginal Gamma inversa para ˜ σ2y la densidad (Normal) de ˜µdada ˜ σ2: f(˜µ, ˜ σ2) = f(˜ σ2)f(˜µ|˜ σ2) = IG(α, β)×N(µ0,1 k˜ σ2). Una vez obtenidos los datos y la media muestral ¯x, podemos aplicar el teorema de Bayes f(˜µ, ˜ σ2|¯x)∝βα Γ(α) 1 (˜σ2)α+1 exp(−β ˜σ2)×√k √2π˜ σ2exp(−k(˜µ−µ0)2 2˜ σ2) ×(˜σ)−ne−ns2/2˜σ2×e−n(˜µ−¯x)2/2˜σ2 ∝1 (˜σ2)α+n 2+1 exp(−1 ˜ σ2(β+ns2 2)) ×√k ˜σexp(−k(˜µ−µ0)2+n(˜µ−¯x)2 2˜ σ2). Haremos ahora algunas manipulaciones en el numerador dentro de la última exponencial: 52 5.4. Distribución predictiva k(˜µ−µ0)2+n(˜µ−¯x)2=k(˜µ2−2˜µµ0+µ2 0) + n(˜µ2−2˜µ¯x+ ¯x2) = (k+n)˜µ2−2˜µ(kµ0+n¯x) + kµ2 0+n¯x2 = (k+n)[˜µ2−2˜µ(kµ0+n¯x) 2+ (kµ0+n¯x k+n)2]−(kµ0+n¯x)2 k+n+kµ2 0+n¯x2 = (k+n)(˜µ−kµ0+n¯x k+n)2+kn(¯x−µ0)2 k+n. Usando esto en la expresión para la distribución conjunta a posteriori, finalmente tenemos: f(˜µ, ˜ σ2|¯x)∝1 (˜σ2)α+n 2+1 exp[−1 ˜ σ2(β+ns2 2+kn(¯x−µ0)2 2(k+n))] ×√k ˜σexp[−(k+n)((˜µ−kµ0+n¯x k+n)2 2˜σ2)]. Lo obtenido en nuestros cálculos corresponde a una distribución Normal-Gamma inversa: ˜µ, ˜ σ2|¯x∼IG(α+n 2, β +ns2 2+kn(¯x−µ0)2 2(k+n)) ×N(kµ0+n¯x k+n,˜σ2 k+n). 5.4. Distribución predictiva Estudiaremos la distribución predictiva a posteriori para una distribución normal con media desconocida y varianza conocida. Sean x1, ..., xnuna muestra aleatoria obtenida de una población normal N(˜µ, σ2). Nuestro propósito es predecir una observación futura independiente Y=Xn+1 de la misma distribución. Usando la distribución a priori de Jeffreys para la media, es decir, f(˜µ)∼1, teníamos la distribución a posteriori ˜µ∼N(¯x, σ2/n). La distribución a posteriori predictiva de Y es, usando la fórmula dada en el capítulo segundo, la siguiente: f(y|x) = Rf(y|˜µ)f(˜µ|x)dµ ∝Rexp{−1 2[(˜µ−y)2 σ2+n(˜µ−¯x)2 σ2]}dµ. Consideremos el siguiente resultado: A(x−a)2+B(x−b)2=C(x−c)2+AB C(a−b)2 donde C=A+B y c=(Aa+Bb)/C. Podemos usar este resultado para obtener (˜µ−y)2 σ2+n(˜µ−¯x)2 σ2=C(˜µ−c)2+(y−¯x)2 (1+ 1 n)σ2 Capítulo 5. Distribución normal 53 donde C= (n+ 1)/σ2yc= (y+n¯x)/(n+ 1). Observemos que, mientras el primer término de la suma es una forma cuadrática de ˜µ, el segundo término no depende del parámetro desconocido. Aplicándo lo anterior a nuestra expresión para la distribución a posteriori obtenemos f(y|x)∝Rexp{−C 2(˜µ−c)2−1 2 (y−¯x)2 (1+ 1 n)σ2}d˜µ =exp{−1 2 (y−¯x)2 (1+ 1 n)σ2}Rexp{−C 2(˜µ−c)2}d˜µ ∝exp{−1 2 (y−¯x)2 (1+ 1 n)σ2}. donde se ha usado que Rexp{−C 2(˜µ−c)2}d˜µ=√2πσ √n+1. Hemos obtenido la función de densidad de una normal de media ¯xy varianza σ2(1 + 1 n), y por tanto la distribución predictiva a posteriori es Y∼N(¯x, σ2(1 + 1 n)). 5.5. Test de hipótesis Estudiaremos algunos contrastes de hipótesis relativos a la distribución normal, y veremos la relación existente entre el enfoque Bayesiano y el p-valor clásico. 5.5.1. Media de una normal con varianza desconocida Consideremos el contraste siguiente: H0:µ≤µ0 H1:µ > µ0 y los estadísticos muestrales ¯x0ys2. El p-valor clásico es P(˜ ¯x−µ0 ˜s/√n≥¯x0−µ0 s/√n) = P(˜ t(τ)≥¯x0−µ0 s/√n). Un análisis bayesiano con distribución a priori no informativa da lugar a la distribución a posteriori t(¯x0, s2/n, τ)para ˜µ. Pt(˜µ≤µ0|¯x0, s2/n, τ) = P(˜µ−¯x0 s/√n≥µ0−¯x0 s/√n) =P(˜ t(τ)≥µ0−¯x0 s/√n). Ya que la distribución tde Student es simétrica con respecto al cero, el área a la derecha, en sentido clásico ¯x0−µ0 s/√n es la misma que el área a la izquierda bajo µ0−¯x0 s/√n. 54 5.5. Test de hipótesis 5.5.2. Varianza Para realizar el contraste siguiente H0:σ2≤σ2 0 H1:σ2> σ2 0 dados datos recogidos en los valores nys2 0, el p-valor clásico es P(˜s2≥s2 0|ν, σ2 0) = P(ν˜s2 σ2 0≥νs2 0 σ2 0|ν, σ2 0) = Pχ2(˜χ2(ν)≥νs2 0 σ2 0). Una aproximación Bayesiana al problema considerando una distribución a priori no informativa resulta en que la distribución a posteriori para (νs2 0 ˜σ2 0|ν, s2 0) es una distribución χ2de parámetro ν, y el interés estaría en calcular lo que sigue: P(˜σ2≤σ2 0|ν, s2 0) = P(νs2 0 ˜σ2>νs2 0 σ2 0) = Pχ2(˜χ2(ν)>νs2 0 σ2 0). Coincide con el p-valor en sentido clásico. Capítulo 6 Inferencia Bayesiana en R En este capítulo se definirán instrucciones para el programa R que permitirán realizar inferencia Bayesiana usando los procedimientos descritos en el presente trabajo. Se empleará el paquete “LearnBayes”, disponible en el CRAN de R. 6.1. Inferencia sobre una proporción Selección de una Beta a priori dados dos cuantiles Descripción: encuentra los parámetros de una densidad Beta que se corresponda con los valores de dos cuantiles de la distribución. Función: beta.select(quantile1, quantile2). Argumentos: quantile1: lista con las componentes p, el valor de la primera probabilidad, y x, el valor del primer cuantil. quantile2: lista con las componentes p, el valor de la segunda probabilidad, y x, el valor del segundo cuantil. Resultado: vector con los parámetros de la distribución Beta. Distribución a posteriori para muestra Binomial y Beta a priori Descripción: Calcula los parámetros y las probabilidades a posteriori para un problema de muestreo binomial donde la distribución a priori es una mezcla discreta de densidades Beta. Función: binomial.beta.mix(probs,betapar,data) 55 56 6.1. Inferencia sobre una proporción Argumentos: probs: vector de probabilidades de las componentes Beta de la distribución a priori. betapar: matriz en la que cada fila contiene los parámetros de una componente Beta de la distribución a priori. data: vector con el número de éxitos y el número de fracasos. Resultado: probs: vector de probabilidades de las componentes Beta a posteriori. betapar: matriz en la que cada fila contiene los parámetros de una componente Beta de la distribución a posteriori. Distribución a posteriori de una proporción con distribución a priori discreta Descripción: calcula la distribución a posteriori de una proporción a partir de una distribución a priori discreta. Función: pdisc(p, prior, data) Argumentos: p: vector con los valores de la proporción. prior: vector con probabilidades a priori. data :vector con el número de éxitos y el número de fracasos. Resultado: vector de probabilidades a posteriori. Distribución predictiva para una muestra binomial con distribución Beta a priori Descripción: calcula la distribución predictiva del número de éxitos para un experimento Binomial con dstribución a priori Beta para la proporción. Función: pbetap(ab, n, s) Argumentos: ab: vector de parámetros de la Beta a priori. n: tamaño de la muestra binomial futura. s: vector con el número de éxitos para el furturo experimento Binomial. Resultado: vector de probabilidades predecidas para los valores del vector s. Distribución predictiva para una muestra binomial con distribución discreta a priori Capítulo 6. Inferencia Bayesiana en R 57 Descripción: calcula la distribución predictiva del número de éxitos para un experimento Binomial con dstribución a priori discreta para la proporción. Función: pdiscp(p, probs, n, s) Argumentos: p: vector de valores de la proporción. probs: vector de probabilidades. n: tamaño de la muestra binomial futura. s: vector con el número de éxitos para el furturo experimento Binomial. Resultado: vector de probabilidades predecidas para los valores del vector s. Test Bayesiano sobre una proporción Descripción: contraste de hipótesis bayesiano de que una proporción sea igual a un valor especificado usando una beta a priori. Función: pbetat(p0,prob,ab,data) Argumentos: p0: valor de la proporción que se quiere contrastar. prob: probabilidad a priori de la hipótesis. ab: vector de parámetros de la Beta a priori bajo la hipótesis alternativa. data: vector con el número de éxitos y el número de fracasos. Resultado: bf: factor Bayes en favor de la hipótesis nula. post: probabilidad a posteriori de la hipótesis nula. 6.2. Distribución de Poisson Distribución a posteriori para una muestra Poisson con distribución Gamma a priori Descripción: calcula los parámetros y las probabilidades para un problema de muestreo de Poisson donde la distribución a priori es una mezcla discreta de densidades gamma.