4.3 Modelo Normal multivariante con media y varianza desconocida

Al igual que en la distribución normal univariada, cuando se desconoce tanto el vector de medias como la matriz de varianzas y covarianzas de la distribución, es necesario plantear diversos enfoques y situarse en el más conveniente. Nótese que en términos de parámetros, existen \(p\) parámetros correspondientes al vector de medias \(\boldsymbol \theta\) y \(\binom{p}{2}=\dfrac{p(p+1)}{2}\) parámetros correspondientes a la matriz de varianzas \(\boldsymbol \Sigma\). Pensando en la gran cantidad de parámetros que se deben modelar, es necesario tener en cuenta que el número de datos en la muestra aleatoria sea lo suficientemente grande. Suponiendo que el número de observaciones en la muestra aleatoria sea suficiente, existe otra situación que se debe surtir y es la asignación de las distribuciones previas para \(\boldsymbol \theta\) y \(\boldsymbol \Sigma\). En estos términos, es posible

  1. Suponer que la distribución previa \(p(\boldsymbol \theta)\) es independiente de la distribución previa \(p(\boldsymbol \Sigma)\) y que ambas distribuciones son informativas. Luego, utilizar un análisis de simulación condicional conjunta para extraer muestras provenientes de las respectivas distribuciones posterior.
  2. Suponer que la distribución previa para \(\boldsymbol \theta\) depende de \(\boldsymbol \Sigma\) y escribirla como \(p(\boldsymbol \theta\mid \boldsymbol \Sigma)\), mientras que la distribución previa de \(\boldsymbol \Sigma\) no depende de \(\boldsymbol \theta\) y se puede escribir como \(p(\boldsymbol \Sigma)\). El análisis posterior de este enfoque encuentra la distribución posterior de \(\boldsymbol \Sigma\mid \mathbf{Y}\) y con esta se encuentra la distribución posterior de \(\boldsymbol \theta\mid \boldsymbol \Sigma,\mathbf{Y}\).
  3. Suponer que la distribución conjunta previa para \(\boldsymbol \theta\) y \(\boldsymbol \Sigma\) es una distribución no informativa.

4.3.1 Parámetros independientes con distribuciones previas informativas

En este enfoque se supone que las distribuciones previas para los parámetros de interés son independientes e informativas. Hacemos siguiente observación para lograr que las resultantes distribuciones posterior sean conjugadas. \[\begin{align*} \sum_{i=1}^n(\mathbf{Y}_i-\theta)'\boldsymbol \Sigma^{-1}(\mathbf{Y}_i-\theta)&=traza \left(\sum_{i=1}^n(\mathbf{Y}_i-\theta)'\boldsymbol \Sigma^{-1}(\mathbf{Y}_i-\theta)\right)\\ &= \sum_{i=1}^n traza\left((\mathbf{Y}_i-\theta)'\boldsymbol \Sigma^{-1}(\mathbf{Y}_i-\theta)\right)\\ &= \sum_{i=1}^n traza\left(\boldsymbol \Sigma^{-1}(\mathbf{Y}_i-\theta)(\mathbf{Y}_i-\theta)'\right)\\ &= traza\left(\boldsymbol \Sigma^{-1}\sum_{i=1}^n(\mathbf{Y}_i-\theta)(\mathbf{Y}_i-\theta)'\right)\\ &= traza\left(\boldsymbol \Sigma^{-1}\mathbf{S}_{\boldsymbol \theta}\right) \end{align*}\]

Donde \(\mathbf{S}_{\boldsymbol \theta}=\sum_{i=1}^n(\mathbf{Y}_i-\boldsymbol \theta)(\mathbf{Y}_i-\boldsymbol \theta)'\). En cuanto a la asignación de las distribuciones previas, para el vector de medias \(\boldsymbol \theta\) es posible usar la distribución normal, esto es, \[\begin{equation*} \boldsymbol \theta\sim Normal_p(\boldsymbol \mu,\boldsymbol \Gamma) \end{equation*}\]

Por otro lado, la distribución para la matriz de varianzas \(\boldsymbol \Sigma\) es \[\begin{equation*} \boldsymbol \Sigma\sim Inversa-Wishart(\boldsymbol \Lambda,v) \end{equation*}\]

donde \(v\) denota los grados de libertad y \(\boldsymbol \Lambda\) la matriz de escala. Esto es, la función de densidad está dada por \[\begin{equation*} p(\boldsymbol \Sigma)\propto |\boldsymbol \Sigma|^{-\frac{v+p+1}{2}}\exp\left\{-\frac{1}{2}traza(\boldsymbol \Lambda\boldsymbol \Sigma^{-1})\right\} \end{equation*}\]

Asumiendo independencia previa, la distribución previa conjunta resulta estar dada por \[\begin{align} p(\boldsymbol \theta,\boldsymbol \Sigma)&=p(\boldsymbol \theta)p(\boldsymbol \Sigma)\notag\\ &\propto \mid \boldsymbol \Sigma\mid ^{-(v+p+1)/2}\notag\\ &\times \exp\left\{ -\frac{1}{2}\left[traza(\boldsymbol \Lambda\boldsymbol \Sigma^{-1})+ (\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Gamma^{-1}(\boldsymbol \theta-\boldsymbol \mu)\right]\right\} \end{align}\]

Una vez que se conoce la forma estructural de la distribución previa conjunta, es posible establecer la distribución posterior conjunta teniendo en cuenta la forma de la función de verosimilitud \(p(\mathbf{Y} \mid \boldsymbol \theta,\boldsymbol \Sigma)\) y la expresión equivalente para \(\sum_{i=1}^n(\mathbf{Y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{Y}_i-\boldsymbol \theta)\) mostrada al inicio de esta sección. Adicionalmente, acudiendo a la simetría de las matrices \(\boldsymbol \Lambda\), \(\boldsymbol \Sigma\) y \(\mathbf{S}_{\boldsymbol \theta}\), se tiene que \[\begin{align} p(\boldsymbol \theta,\boldsymbol \Sigma\mid \mathbf{Y})&\propto p(\boldsymbol \theta,\boldsymbol \Sigma)p(\mathbf{Y} \mid \boldsymbol \theta,\boldsymbol \Sigma)\notag\\ &\propto \mid \boldsymbol \Sigma\mid ^{-(v+n+p+1)/2}\notag\\ &\times \exp\left\{ -\frac{1}{2}\left[traza(\boldsymbol \Lambda\boldsymbol \Sigma^{-1}+\boldsymbol \Sigma^{-1}\mathbf{S}_{\boldsymbol \theta})+ (\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Gamma^{-1}(\boldsymbol \theta-\boldsymbol \mu)\right]\right\}\notag\\ &\propto \mid \boldsymbol \Sigma\mid ^{-(v+n+p+1)/2}\notag\\ &\times \exp\left\{ -\frac{1}{2}\left[traza(\boldsymbol \Sigma^{-1}(\boldsymbol \Lambda+\mathbf{S}_{\boldsymbol \theta}))+ (\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Gamma^{-1}(\boldsymbol \theta-\boldsymbol \mu)\right]\right\} \end{align}\]

Dado que la distribución posterior conjunta no tiene una forma estructural conocida, no es posible utilizar el método de integración analítica. Sin embargo, es posible obtener las distribuciones condicionales de cada uno de los parámetros suponiendo fijos los restantes y teniendo en cuenta que \[\begin{align*} p(\boldsymbol \theta\mid \boldsymbol \Sigma,\mathbf{Y})\propto p(\boldsymbol \theta,\underbrace{\boldsymbol \Sigma}_{fijo} \mid \mathbf{Y}) \ \ \ \ \ \ \ \ \ \text{y} \ \ \ \ \ \ \ \ \ \ p(\boldsymbol \Sigma\mid \boldsymbol \theta,\mathbf{Y})\propto p(\underbrace{\boldsymbol \theta}_{fijo},\boldsymbol \Sigma\mid \mathbf{Y}) \end{align*}\]

Resultado 4.10 La distribución posterior de la matriz de parámetros \(\boldsymbol \Sigma\) condicional a \(\boldsymbol \theta,\mathbf{Y}\) es \[\begin{equation*} \tag{4.8} \boldsymbol \Sigma\mid \boldsymbol \theta,\mathbf{Y} \sim Inversa-Wishart_{v+n}(\boldsymbol \Lambda+\mathbf{S}_{\boldsymbol \theta}) \end{equation*}\]


Prueba. La prueba es inmediata notando que \[\begin{align*} \boldsymbol \Sigma\mid \boldsymbol \theta,\mathbf{Y}&\propto\mid \boldsymbol \Sigma\mid ^{-(v+n+p+1)/2}\notag\\ &\times\exp\left\{ -\frac{1}{2}\left[traza(\boldsymbol \Sigma^{-1}(\boldsymbol \Lambda+\mathbf{S}_{\boldsymbol \theta}))+(\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Gamma^{-1}(\boldsymbol \theta-\boldsymbol \mu)\right]\right\} \end{align*}\]

Por lo tanto, factorizando convenientemente, se encuentra una expresión idéntica a la función de distribución de una variable aleatoria con distribución \(Inversa-Wishart_{v+n}(\boldsymbol \Lambda+\mathbf{S}_{\boldsymbol \theta})\).


Resultado 4.11 La distribución posterior del vector de parámetros \(\boldsymbol \theta\) condicional a \(\boldsymbol \Sigma,\mathbf{Y}\) es \[\begin{equation} \tag{4.9} \boldsymbol \theta\mid \boldsymbol \Sigma,\mathbf{Y} \sim Normal_p(\boldsymbol \mu_n,\boldsymbol \Gamma_n) \end{equation}\]

donde \(\boldsymbol \mu_n\) y \(\boldsymbol \Gamma_n\) están dadas por las expresiones (4.6) y (4.7), respectivamente.


Una vez encontradas las distribuciones posteriores condicionales de \(\boldsymbol \theta\) y \(\boldsymbol \Sigma\), se puede obtener la estimación de estos parámetros vía el muestreo de Gibbs, que en este caso se resume en los siguientes pasos:

  1. Fijar un valor inicial para \(\boldsymbol \theta\); lo denotamos por \(\boldsymbol \theta_{(1)}\).
  2. Simular un valor de la distribución de \(\boldsymbol \Sigma|\boldsymbol \theta,\mathbf{Y}\) en (4.8) donde el parámetro \(\mathbf{S_\mathbf{\boldsymbol \theta}}\) que depende de \(\boldsymbol \theta\), debe ser reemplazado por \(\boldsymbol \theta_{(1)}\) del paso anterior; este valor simulado se denotará por \(\boldsymbol \Sigma_{(1)}\)
  3. Simular un valor de la distribución de \(\boldsymbol \theta|\boldsymbol \Sigma,\mathbf{Y}\) en (4.9) donde en \(\mathbf{mu}_n\) y \(\boldsymbol \Gamma_n\) se debe reemplazar \(\boldsymbol \Sigma\) por \(\boldsymbol \Sigma\); este valor simulado se denota por \(\boldsymbol \theta\).
  4. Repetir los pasos (2) y (3) hasta completar un número de iteraciones suficientes para alcanzar la convergencia en ambos parámetros

Una vez tengamos los valores muestreados, se debe garantizar la convergencia y la correlación nula entre estos valores, con el fin de calcular las estimaciones. En el siguiente ejemplo ilustramos la implementación de este muestreo de Gibbs en R.

Ejemplo 4.3 Retomamos los datos del efecto de dos medicamentos soporíferos introducidos por Student (1908), los cuales fueron estudiados en el ejemplo 4.2 asumiendo que la matriz de varianzas y covarianzas era conocida. El vector de medias muestrales de estos datos están dados por \(\bar{y}=(0.75, 2.33)'\), y la matriz de varianzas y covarianzas muestrales está dada por \(\mathbf{S}=\begin{pmatrix}3.20&2.85\\2.85&4.01\end{pmatrix}\).

Ahora supongamos que tanto el vector de medias como la matriz de varianzas y covarianzas son desconocidos. Para el vector de medias, asumimos la distribución previa del ejemplo , es decir, \(\boldsymbol \mu=(0, 1)'\) y \(\boldsymbol \Gamma=\begin{pmatrix}2&0\\0&2\end{pmatrix}\). Para la matriz de varianzas y covarianzas asumimos la distribución inversa-Wishart con matriz de escala igual a \(\boldsymbol \Lambda=\begin{pmatrix}20&8\\8&20\end{pmatrix}\) y \(v=10\) grados de libertad. De esta forma, la estimación previa de \(\boldsymbol \Sigma\) viene dada por \(\frac{1}{v-2-1}\boldsymbol \Lambda=\begin{pmatrix}2.86&1.14\\ 1.14&2.86\end{pmatrix}\).

Ilustramos los códigos de STAN a continuación. 
NormalMultMediaCov <- '
data {
  int<lower=0> n;
  int<lower=0> P;
  vector[P] y[n];
  vector[P] mu;
  matrix[P, P] Gamma;
  matrix[P, P] Lambda;
  int<lower=0> v;
}
parameters {
  vector[P] theta;
  cov_matrix[P] Sigma;
}
transformed parameters {
  real diftheta;
  diftheta = theta[2] - theta[1];
}
model {
  theta ~ multi_normal(mu, Gamma);
  Sigma ~ inv_wishart(v, Lambda);
  for (i in 1:n)
    y[i] ~ multi_normal(theta, Sigma);
}
'

y <- structure(.Data = sleep[,1], .Dim=c(10,2))
n <- nrow(y)
P <- ncol(y)
Sigma  <- matrix(c(1, 0.6, 0.6, 2), 2, 2)
mu <- as.vector(c(0, 1))
Gamma <- matrix(c(2, 0, 0, 2), 2, 2)
v <- 10
Lambda <- matrix(c(20, 8, 8, 20), 2, 2)

sample_data <- list(y = y, n = n, P = P,
                    Sigma = Sigma, mu = mu,
                    Gamma = Gamma, v = v,
                    Lambda = Lambda)
set.seed(1234)
NormalMultMediaCovfit <- stan(model_code = NormalMultMediaCov,
                   data = sample_data, verbose = FALSE)

Con base en los resultados de las cadenas simuladas, se observa que la estimación bayesiana para el número de horas de sueño producidas por los dos medicamentos son 0.5772 y 2.1011, respectivamente. En cuanto a la estimación de la matriz de varianzas y covarianzas, ésta está dada por \(\hat{\boldsymbol \Sigma}=\begin{pmatrix}3.0132&2.1011\\2.1011&3.5222\end{pmatrix}\).

print(NormalMultMediaCovfit, digits = 4, 
      pars = c("theta", "Sigma"), probs = c(0.025, 0.975))
## Inference for Stan model: 8ef9836716ef08fd4cd1136cfef43904.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##              mean se_mean     sd    2.5%  97.5% n_eff   Rhat
## theta[1]   0.5482  0.0103 0.5001 -0.4113 1.5300  2337 1.0001
## theta[2]   2.0732  0.0122 0.5480  0.8840 3.0624  2033 1.0008
## Sigma[1,1] 3.0355  0.0241 1.1310  1.5812 5.8057  2206 0.9999
## Sigma[1,2] 2.1079  0.0232 1.0193  0.7354 4.7245  1928 1.0002
## Sigma[2,1] 2.1079  0.0232 1.0193  0.7354 4.7245  1928 1.0002
## Sigma[2,2] 3.5141  0.0280 1.3193  1.8198 6.8809  2214 0.9995
## 
## Samples were drawn using NUTS(diag_e) at Sun Jun 27 20:30:22 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Las figuras 4.7 y 4.8 muestran las distribuciones posteriores en este ejemplo.

Distribuciones posteriores para el vector de medias y la diferencia de medias.

Figura 4.7: Distribuciones posteriores para el vector de medias y la diferencia de medias.

Distribuciones posteriores para los elementos de la matriz de covarinzas.

Figura 4.8: Distribuciones posteriores para los elementos de la matriz de covarinzas.

A continuación se muestran los códigos necesarios para implementar el muestreo de Gibbs de forma manual en R.

library(MCMCpack)
library(mvtnorm)

y.bar <- colMeans(y)
n <- nrow(y)
nsim <- 10000

theta.pos <- matrix(NA, nsim, P)
Sigma.pos <- array(NA, c(nsim, P, P))

# Valor inicial de theta
theta.pos[1,] <- c(0, 1)

# Parámetros posteriores de Sigma
v.pos <- v + n
matrix.theta <- kronecker(matrix(rep(1, n)),
                          t(theta.pos[1, ])) 
S.theta <- t(y - matrix.theta) %*% (y-matrix.theta)
Lambda.pos <- Lambda + S.theta
# Simulación de la distribución posterior condicional de Sigma
Sigma.pos[1, , ] <- riwish(v.pos, Lambda.pos)

#####################
# muestreo de Gibbs #
#####################

for(i in 2:nsim){
  # Parámetros posteriores de theta 
  Gamma.n <- solve(solve(Gamma) + 
                     n * solve(Sigma.pos[i - 1, , ]))
  mu.n <- Gamma.n %*%
    (solve(Gamma) %*% mu + 
       n * solve(Sigma.pos[i - 1, , ]) %*% y.bar)
  # Simulación de la distribucion posterior condicional de theta
  theta.pos[i, ] <- rmvnorm(1, mu.n, Gamma.n)
  # Parámetros posteriores de Sigma
  matrix.theta <- kronecker(matrix(rep(1, n)),
                          t(theta.pos[i, ]))  
  S.theta <- t(y - matrix.theta) %*% (y - matrix.theta)
  Lambda.pos <- Lambda + S.theta
  # Simulación de la distribución posterior condicional de Sigma
  Sigma.pos[i, , ] <- riwish(v.pos, Lambda.pos)
}

Una vez finalizada la ejecución del muestreo de Gibbs, debemos examinar la calidad de los valores muestreados para asegurar que las estimaciones bayesianas sean obtenidas de una muestra de valores que hayan convergido, y en segundo lugar que no estén correlacionados. Para eso a continuación observamos la gráfica de los valores generados para algunos parámetros (en particular, consideramos los parámetros \(\theta_1\), \(\theta_2\), \(\sigma^2_1\) y \(\sigma_{12}\)), así como la gráfica de las autocorrelaciones muestrales.

Convergencia de las distribuciones posteriores y diagramas de la función de autocorrelación en las cadenas.

Figura 4.9: Convergencia de las distribuciones posteriores y diagramas de la función de autocorrelación en las cadenas.

Con estas gráficas, observamos que los valores muestreados han alcanzado la convergencia; además estos tienen correlaciones cercanas a cero. De esta forma, podemos usar los valores muestreados para calcular las estimaciones y los intervalos de credibilidad.

theta.Bayes <- colMeans(theta.pos)
Sigma.Bayes <- matrix(c(mean(Sigma.pos[,1,1]),
                        mean(Sigma.pos[,2,1]),
                        mean(Sigma.pos[,1,2]),
                        mean(Sigma.pos[,2,2])),
                      nrow = 2, ncol = 2)
theta.Bayes
## [1] 0.5638342 2.0993759
Sigma.Bayes
3.042341 2.093069
2.093069 3.490823

El procedimiento inferencial sobre la comparación entre los efectos de los dos medicamentos se puede realizar de la misma manera como ilustró el ejemplo 4.2.

4.3.2 Parámetros dependientes

Al igual que en el caso univariado, la inferencia posterior de los parámetros de interés debe ser llevada a cabo en dos etapas: En la primera, se debe establecer la distribución previa conjunta para ambos parámetros mediante \[\begin{equation*} p(\boldsymbol \theta,\boldsymbol \Sigma)=p(\boldsymbol \Sigma)p(\boldsymbol \theta\mid \boldsymbol \Sigma) \end{equation*}\]

Luego, en la segunda etapa es posible analizar posterior propiamente cada uno de los parámetros de interés puesto que \[\begin{equation*} p(\boldsymbol \theta,\boldsymbol \Sigma\mid \mathbf{Y})\propto p(\mathbf{Y} \mid \boldsymbol \theta,\boldsymbol \Sigma)p(\boldsymbol \theta,\boldsymbol \Sigma) \end{equation*}\]

Al igual que en el caso univariado, la anterior formulación conlleva a asignar una distribución previa para \(\boldsymbol \theta\) dependiente de la matriz \(\boldsymbol \Sigma\). Esto quiere decir que en la distribución \(p(\boldsymbol \theta\mid \boldsymbol \Sigma)\) el valor de \(\boldsymbol \Sigma\) se considera una constante fija y conocida. Siguiendo los lineamientos del capítulo anterior, una distribución previa para \(\boldsymbol \theta\) condicional a \(\boldsymbol \Sigma\) es \[\begin{equation*} p(\boldsymbol \theta\mid \boldsymbol \Sigma)\sim Normal_p(\boldsymbol \mu,\boldsymbol \Sigma/c_0) \end{equation*}\]

Donde \(c_0\) es una constante. Por otro lado, y siguiendo los argumentos de la sección anterior, una posible opción para la distribución previa de \(\boldsymbol \Sigma\), corresponde a \[\begin{equation*} p(\boldsymbol \Sigma)\sim Inversa-Wishart_{v_0}(\boldsymbol \Lambda) \end{equation*}\]

Resultado 4.12 La distribución previa conjunta de los parámetros \(\boldsymbol \theta\) y \(\boldsymbol \Sigma\) está dada por \[\begin{equation*} p(\boldsymbol \theta,\boldsymbol \Sigma) \propto \mid \boldsymbol \Sigma\mid ^{-(v_0+p)/2-1} \exp\left\{ -\frac{1}{2}\left[traza(\boldsymbol \Lambda_0\boldsymbol \Sigma^{-1})+ c_0(\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\boldsymbol \mu)\right]\right\} \end{equation*}\]


Prueba. La prueba es inmediata al multiplicar las densidades y asignar los términos que no dependen de los parámetros de interés a la constante de proporcionalidad.


Para encontrar las distribuciones posteriores de cada uno de los parámetros de interés se utilizan argumentos similares a los del capítulo anterior.

Resultado 4.13 La distribución posterior de \(\boldsymbol \theta\) condicional a \(\boldsymbol \Sigma,\mathbf{Y}\) está dada por \[\begin{equation*} \theta \mid \boldsymbol \Sigma,\mathbf{Y} \sim Normal_p(\boldsymbol \mu_n,\boldsymbol \Sigma/(n+c_0)) \end{equation*}\]

donde \[\begin{equation*} \mu_n=\frac{n\bar{\mathbf{Y}}+c_0\boldsymbol \mu}{n+c_0} \end{equation*}\]


Prueba. Utilizando propiedades de la distribución condicional, tenemos que \[\begin{align*} p(\boldsymbol \theta|\boldsymbol \Sigma,\mathbf{Y})&\propto p(\boldsymbol \theta, \boldsymbol \Sigma|\mathbf{Y})\\ &\propto \mid \boldsymbol \Sigma\mid ^{-(v_0+p)/2-1} \exp\left\{ -\frac{1}{2}\left[traza(\boldsymbol \Lambda_0\boldsymbol \Sigma^{-1})+ c_0(\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\boldsymbol \mu)\right]\right\}\\ &\ \ \ \ \ \ \ \ \mid \boldsymbol \Sigma\mid ^{-n/2}\exp\left\{-\frac{1}{2}\sum_{i=1}^n(\mathbf{y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{y}_i-\boldsymbol \theta)\right\}\\ &\propto \exp\left\{ -\frac{1}{2} c_0(\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\boldsymbol \mu)\right\}\exp\left\{-\frac{1}{2}\sum_{i=1}^n(\mathbf{y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{y}_i-\boldsymbol \theta)\right\} \end{align*}\]

La anterior expresión es la misma que apareción en el capítulo anterior para \(p(\boldsymbol \theta\mid\boldsymbol \Sigma,\mathbf{Y})\), en donde \(\boldsymbol \Sigma/c_0\) toma el valor de \(\boldsymbol \Gamma\). Así, teniendo en cuenta las ecuaciones (4.6) y (4.7), podemos afirmar que el vector de medias y la matriz de varianzas y covarianzas posterior están dadas por \[\begin{align} \boldsymbol \Gamma_n&=((\boldsymbol \Sigma/c_0)^{-1}+n\boldsymbol \Sigma^{-1})^{-1}=\frac{\boldsymbol \Sigma}{n+c_0}\\ \boldsymbol \mu_n&=\frac{\boldsymbol \Sigma}{n+c_0}((\boldsymbol \Sigma/c_0)^{-1}\boldsymbol \mu+n\boldsymbol \Sigma^{-1}\bar{\mathbf{y}})=\frac{n\bar{\mathbf{Y}}+c_0\boldsymbol \mu}{n+c_0} \end{align}\]


En cuanto a la distribución de \(\boldsymbol \Sigma\), se tiene el siguiente resultado:

Resultado 4.14 La distribución marginal posterior de la matriz de parámetros \(\boldsymbol \Sigma\) es \[\begin{equation*} \boldsymbol \Sigma\mid \mathbf{Y} \sim Inversa-Whishart_{n+v_0}(\boldsymbol \Lambda_n) \end{equation*}\]

Donde \[\begin{equation}\label{bLambda_n} \boldsymbol \Lambda_n=\boldsymbol \Lambda+(n-1)\mathbf{S}+\frac{c_0n}{c_0+n}(\boldsymbol \mu-\bar{\mathbf{y}})(\boldsymbol \mu-\bar{\mathbf{y}})' \end{equation}\]

con \(S\) la matriz de varianzas y covarianzas muestrales.


Prueba. \[\begin{align*} &\ \ \ \ \ p(\boldsymbol \Sigma\mid\mathbf{Y})\\ &=\int_{R^p} p(\boldsymbol \theta,\boldsymbol \Sigma\mid\mathbf{Y})d\boldsymbol \theta\\ &\propto \int_{R^p}\mid\boldsymbol \Sigma\mid^{-(v_0+p+n)/2-1}\exp\left\{-\frac{1}{2}\left[traza(\boldsymbol \Lambda\boldsymbol \Sigma^{-1})+c_0(\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\boldsymbol \mu)+\sum_{i=1}^n(\mathbf{y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{y}_i-\boldsymbol \theta)\right]\right\}d\boldsymbol \theta\\ &\propto \mid\boldsymbol \Sigma\mid^{-(v_0+p+n)/2-1} \exp\left\{-\frac{1}{2}\left[traza(\boldsymbol \Lambda\boldsymbol \Sigma^{-1})\right]\right\}\\ &\ \ \ \ \ \ \ \ \ \ \int_{R^p}\exp\left\{-\frac{1}{2}\left[c_0(\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\boldsymbol \mu)+\sum_{i=1}^n(\mathbf{y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{y}_i-\boldsymbol \theta)\right]\right\}d\boldsymbol \theta\\ &\propto \mid\boldsymbol \Sigma\mid^{-(v_0+p+n)/2-1} \exp\left\{-\frac{1}{2}\left[traza(\boldsymbol \Lambda\boldsymbol \Sigma^{-1})+c_0\boldsymbol \mu'\boldsymbol \Sigma^{-1}\boldsymbol \mu+\sum_{i=1}^n(\mathbf{y}_i-\bar{\mathbf{y}})'\boldsymbol \Sigma^{-1}(\mathbf{y}_i-\bar{\mathbf{y}})+n\bar{\mathbf{y}}'\boldsymbol \Sigma^{-1}\bar{\mathbf{y}}\right]\right\}\\ &\ \ \ \ \ \ \ \ \ \ \int_{R^p}\exp\left\{-\frac{1}{2}\left[c_0\boldsymbol \theta'\boldsymbol \Sigma^{-1}\boldsymbol \theta-2c_0\boldsymbol \mu'\boldsymbol \Sigma^{-1}\boldsymbol \theta-2n\bar{\mathbf{y}}'\boldsymbol \Sigma^{-1}\boldsymbol \theta+n\boldsymbol \theta'\boldsymbol \Sigma^{-1}\boldsymbol \theta\right]\right\}d\boldsymbol \theta\\ &\propto \mid\boldsymbol \Sigma\mid^{-(v_0+p+n)/2-1} \exp\left\{-\frac{1}{2}\left[traza\left((\boldsymbol \Lambda+c_0\boldsymbol \mu\boldsymbol \mu'+\sum_{i=1}^n(\mathbf{y}_i-\bar{\mathbf{y}})(\mathbf{y}_i-\bar{\mathbf{y}})'+n\bar{\mathbf{y}}\bar{\mathbf{y}}')\boldsymbol \Sigma^{-1}\right)\right]\right\}\\ &\ \ \ \mid\frac{\boldsymbol \Sigma}{c_0+n}\mid^{1/2}\exp\left\{\frac{1}{2}\frac{c_0\boldsymbol \mu'+n\bar{\mathbf{y}}'}{c_0+n}\left(\frac{\boldsymbol \Sigma}{c_0+n}\right)^{-1}\frac{c_0\boldsymbol \mu+n\bar{\mathbf{y}}}{c_0+n}\right\}\\ &\ \ \ \ \ \ \ \underbrace{\int_{R^p}\mid\frac{\boldsymbol \Sigma}{c_0+n}\mid^{-1/2}\exp\left\{-\frac{1}{2}\left(\boldsymbol \theta-\frac{c_0\boldsymbol \mu+n\bar{\mathbf{y}}}{c_0+n}\right)'\left(\frac{\boldsymbol \Sigma}{c_0+n}\right)^{-1}\left(\boldsymbol \theta-\frac{c_0\boldsymbol \mu+n\bar{\mathbf{y}}}{c_0+n}\right)\right\}d\boldsymbol \theta}_{\text{Igual a 1}} \end{align*}\]

Por otro lado, \[\begin{align*} &\ \ \ \ \frac{c_0\boldsymbol \mu'+n\bar{\mathbf{y}}'}{c_0+n}\left(\frac{\boldsymbol \Sigma}{c_0+n}\right)^{-1}\frac{c_0\boldsymbol \mu+n\bar{\mathbf{y}}}{c_0+n}\\ &=\frac{1}{c_0+n}(c_0\boldsymbol \mu'+n\bar{\mathbf{y}}')\boldsymbol \Sigma^{-1}(c_0\boldsymbol \mu+n\bar{\mathbf{y}})\\ &=traza\left(\frac{1}{c_0+n}(c_0\boldsymbol \mu+n\bar{\mathbf{y}})(c_0\boldsymbol \mu'+n\bar{\mathbf{y}}')\boldsymbol \Sigma^{-1}\right)\\ &=traza\left(\left(\frac{c_0^2\boldsymbol \mu\boldsymbol \mu'}{c_0+n}+\frac{2c_0n\bar{\mathbf{y}}\boldsymbol \mu'}{c_0+n}+\frac{n^2\bar{\mathbf{y}}\bar{\mathbf{y}}'}{c_0+n}\right)\boldsymbol \Sigma^{-1}\right) \end{align*}\]

Reemplazando la anterior expresión en \(p(\boldsymbol \Sigma\mid\mathbf{Y})\), se tiene que \[\begin{align*} &\ \ \ \ \ p(\boldsymbol \Sigma\mid\mathbf{Y})\\ &\propto \mid\boldsymbol \Sigma\mid^{-(v_0+p+n+1)/2}\exp\left\{-\frac{1}{2}traza\left[\left(\boldsymbol \Lambda+(n-1)\mathbf{S}+\frac{c_0n}{c_0+n}(\boldsymbol \mu-\bar{\mathbf{y}})(\boldsymbol \mu-\bar{\mathbf{y}})'\right)\boldsymbol \Sigma^{-1}\right]\right\} \end{align*}\]

la cual corresponde a la distribución deseada.


En términos de simulación de densidades, para obtener las estimaciones bayesians de \(\boldsymbol \theta\) y \(\boldsymbol \Sigma\) se debe primero simular valores de \(\boldsymbol \Sigma\) de la distribución \(p(\boldsymbol \Sigma\mid \mathbf{Y})\) y luego, se debe utilizar estos valores para simular valores de \(\boldsymbol \theta\) de la distribución \(p(\boldsymbol \theta\mid \boldsymbol \Sigma,\mathbf{Y})\).

Una forma equivalente de obtener las estimaciones es calcular directamente la esperanza teórica de las distribuciones posteriores marginales de \(\boldsymbol \theta\) y de \(\boldsymbol \Sigma\).

Del resultado 4.14, podemos concluir que la estimación bayesiana de la matriz de varianza y covarianzas \(\boldsymbol \Sigma\) está dada por \[\begin{equation*} \hat{\boldsymbol \Sigma}=\dfrac{\boldsymbol \Lambda+(n-1)\mathbf{S}+\frac{c_0n}{c_0+n}(\boldsymbol \mu-\bar{\mathbf{y}})(\boldsymbol \mu-\bar{\mathbf{y}})'}{n+v_0-p-1} \end{equation*}\]

Teniendo en cuenta que la estimación previa de \(\boldsymbol \Sigma\) viene dada por \(\hat{\boldsymbol \Sigma}_{pre}=\frac{\boldsymbol \Lambda}{v_0-p-1}\), podemos ver que la estimación bayesiana de \(\boldsymbol \Sigma\) está conformada por tres componentes: la estimación previa \(\hat{\boldsymbol \Sigma}_{pre}\), la estimación clásica \(\mathbf{S}\) y una medida de discrepancia entre la estimación previa y la clásica de \(\boldsymbol \theta\). Para encontrar correctas formas de escoger los parámetros previas de \(\boldsymbol \Sigma\), por ahora ignoramos el último componente, y vemos que la estimación previa \(\hat{\boldsymbol \Sigma}_{pre}\) y la estimación clásica \(\mathbf{S}\) entran al cómputo de la estimación bayesiana con los pesos de \(v_0-p-1\) y \(n-1\), de esta forma, podemos escoger \(v_0\) tal que \(v_0-p\) represente el número de la información previa, y el valor de \(\boldsymbol \Lambda\) se puede calcular a partir de \(v_0\) y \(\hat{\boldsymbol \Sigma}_{pre}\).

El siguiente resultado muestra la distribución posterior marginal de \(\boldsymbol \theta\).

Resultado 4.15 La distribución marginal posterior del parámetro \(\boldsymbol \theta\) es la distribución \(t\) de Student multivariante tal que \[\begin{equation*} \boldsymbol \theta\mid \mathbf{Y} \sim t_{n+v_0-p+1}\left(\boldsymbol \mu_n, \frac{\boldsymbol \Lambda_n}{(c_0+n)(n+v_0-p+1)}\right) \end{equation*}\]

con \(\boldsymbol \mu_n=\frac{c_0\boldsymbol \mu+n\bar{\mathbf{y}}}{c_0+n}\) y \(\boldsymbol \Lambda_n\) dado en la ecuación ().


Prueba. \[\begin{align*} &\ \ \ \ p(\boldsymbol \theta\mid\mathbf{Y})\\ &=\int_{R^p\times R^p}p(\boldsymbol \theta,\boldsymbol \Sigma\mid\mathbf{Y})d\boldsymbol \Sigma\\ &=\int_{R^p\times R^p} \mid\boldsymbol \Sigma\mid^{-(v_0+p+n)/2-1}\exp\left\{-\frac{1}{2}\left[traza(\boldsymbol \Lambda\boldsymbol \Sigma^{-1})+c_0(\boldsymbol \theta-\boldsymbol \mu)'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\boldsymbol \mu)+\sum_{i=1}^n(\mathbf{y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{y}_i-\boldsymbol \theta)\right]\right\}d\boldsymbol \Sigma\\ &=\int_{R^p\times R^p} \mid\boldsymbol \Sigma\mid^{-(v_0+p+n+2)/2}\exp\left\{-\frac{1}{2}traza\left[\boldsymbol \Lambda+c_0(\boldsymbol \theta-\boldsymbol \mu)(\boldsymbol \theta-\boldsymbol \mu)'+\sum_{i=1}^n(\mathbf{y}_i-\boldsymbol \theta)(\mathbf{y}_i-\boldsymbol \theta)'\right]\boldsymbol \Sigma^{-1}\right\}d\boldsymbol \Sigma\\ &\propto \big|\boldsymbol \Lambda+c_0(\boldsymbol \theta-\boldsymbol \mu)(\boldsymbol \theta-\boldsymbol \mu)'+\sum_{i=1}^n(\mathbf{y}_i-\boldsymbol \theta)(\mathbf{y}_i-\boldsymbol \theta)'\big|^{-\frac{v_0+n+1}{2}}\\ &=\big|\boldsymbol \Lambda+c_0(\boldsymbol \theta-\boldsymbol \mu)(\boldsymbol \theta-\boldsymbol \mu)'+\sum_{i=1}^n(\mathbf{y}_i-\bar{\mathbf{y}})(\mathbf{y}_i-\bar{\mathbf{y}})'+n(\bar{\mathbf{y}}-\boldsymbol \theta)(\bar{\mathbf{y}}-\boldsymbol \theta)'\big|^{-\frac{v_0+n+1}{2}}\\ &=\big|\boldsymbol \Lambda+(n-1)\mathbf{S}+\frac{c_0n}{c_0+n}(\boldsymbol \mu-\bar{\mathbf{y}})(\boldsymbol \mu-\bar{\mathbf{y}})'+(c_0+n)(\boldsymbol \theta-\boldsymbol \mu_n)(\boldsymbol \theta-\boldsymbol \mu_n)'\big|^{-\frac{v_0+n+1}{2}}\\ &=\big|\boldsymbol \Lambda_n+(c_0+n)(\boldsymbol \theta-\boldsymbol \mu_n)(\boldsymbol \theta-\boldsymbol \mu_n)'\big|^{-\frac{v_0+n+1}{2}}\\ &\propto\big|\mathbf{I}_p+(c_0+n)\boldsymbol \Lambda_n^{-1}(\boldsymbol \theta-\boldsymbol \mu_n)(\boldsymbol \theta-\boldsymbol \mu_n)'\big|^{-\frac{v_0+n+1}{2}}\\ &=\big|1+(c_0+n)(\boldsymbol \theta-\boldsymbol \mu_n)'\boldsymbol \Lambda_n^{-1}(\boldsymbol \theta-\boldsymbol \mu_n)\big|^{-\frac{v_0+n+1}{2}}\\ &=\big|1+\frac{1}{n+v_0-p+1}(\boldsymbol \theta-\boldsymbol \mu_n)'\left(\frac{\boldsymbol \Lambda_n}{(c_0+n)(n+v_0-p+1)}\right)^{-1}(\boldsymbol \theta-\boldsymbol \mu_n)\big|^{-\frac{v_0+n+1}{2}} \end{align*}\]

Esta expresión obtenida corresponde a la forma de la distribución \(t\) de Student multivariado. En el desarrollo se utilizó la propiedad \(|\mathbf{I}+\mathbf{A}\mathbf{B}|=|\mathbf{I}+\mathbf{B}\mathbf{A}|\) para matrices \(\mathbf{A}\) y \(\mathbf{B}\) de tamaños compatibles para las multiplicaciones.


El anterior resultado indica que la estimación bayesiana del parámetro \(\boldsymbol \theta\) está dada por \[\begin{equation*} \hat{\boldsymbol \theta}=\mathbf{\mu}_n=\dfrac{n\hat{\mathbf{Y}}+c_0\mathbf{\mu}}{n+c_0}=\dfrac{n}{n+c_0}\hat{\mathbf{Y}}+\dfrac{c_0}{n+c_0}\mathbf{\mu} \end{equation*}\]

donde se puede observar que \(\hat{\boldsymbol \theta}\) se acercará a la estimación clásica \(\hat{\mathbf{y}}\) cuando \(n\) es grande comparado a \(c_0\), de lo contrario se acercará a la estimación previa \(\mathbf{\mu}\). La varianza posterior para el \(i\)-ésimo componente de \(\boldsymbol \theta\) está dada por \[\begin{equation*} var(\theta_i|\mathbf{Y})=\dfrac{\lambda_{ii}}{(c_0+n)(n+v_0-p+1)}\dfrac{n+c_0-p+1}{n+c_0-p-1}\approx\dfrac{\lambda_{ii}}{(c_0+n)(n+v_0-p+1)} \end{equation*}\]

donde \(\lambda_{ii}\) denota el \(i\)-ésimo elemento en la diagonal de la matriz \(\boldsymbol \Lambda_n\).

Ejemplo 4.4 Peña (2002) reporta las mediciones de 6 variables indicadoras de desarrollo en 91 paises en los años noventa. Para este ejemplo, utilizamos tres variables: la tasa de natalidad, la tasa de mortalidad y la mortalidad infantil en algunos paises de Suramérica y Asia mostrados en la tabla 4.1. Específicamente, usaremos los datos de los paises de Suramérica como datos muestrales y los de Asia para extraer la información previa.
Tabla 4.1: Tasa de natalidad, tasa de mortalidad, mortalidad infantil en algunos países.
Natalidad Mortalidad MortalidadInfantil
20.7 8.4 25.7
46.6 18.0 111.0
28.6 7.9 63.0
23.4 5.8 17.1
27.4 6.1 40.0
32.9 7.4 63.0
29.0 23.2 43.0
34.8 6.6 42.0
32.9 8.3 109.9
18.0 9.6 21.9
27.5 4.4 23.3

A continuación se muestra el proceso necesario para obtener la estimación bayesiana del vector de medias y de la matriz de covarianzas con los datos del ejemplo.

# Datos muestrales
y.sam <- data.frame(
  Nata = c(20.7, 46.6, 28.6, 23.4, 27.4, 
           32.9, 29, 34.8, 32.9,18,27.5), 
  Mort = c(8.4, 18, 7.9, 5.8, 6.1, 7.4,
           23.2, 6.6, 8.3, 9.6, 4.4),
  Infa = c(25.7, 111, 63, 17.1, 40, 63,
           43, 42, 109.9, 21.9, 23.3))
# Datos de la información previa
y.pre <- data.frame(
  Nata = c(21.2, 30.5, 28.6, 31.6, 
           36.1, 39.6, 17.8),
  Mort = c(6.7, 10.2, 9.4, 5.6,
           8.8, 14.8, 5.2),
  Infa=c(32, 91, 75, 24, 68, 128, 7.5))

n <- nrow(y.sam)
P <- ncol(y.sam)
# Estimación clásica de los parámetros
y.bar <- colMeans(y.sam) 
S <- var(y.sam)

# Estimación previa de los parámetros
mu <- colMeans(y.pre)
c0 <- nrow(y.pre)
v0 <- P + nrow(y.pre)
Lambda <- var(y.pre) * (v0 - P - 1)

# Parámetros de las distribuciones posteriores marginales
mu.n <- (n * y.bar + c0 * mu)/(n + c0)
Lambda.n <- Lambda + (n - 1) * S + 
  matrix(mu - y.bar) %*% 
  t(matrix(mu - y.bar)) * c0 * n
var.theta <- Lambda.n/((c0 + n) * (n + v0 - P + 1))
mu.n
##      Nata      Mort      Infa 
## 29.288889  9.244444 54.744444
var.theta
Nata Mort Infa
Nata 2.7991957 0.7890556 10.675752
Mort 0.7890556 1.3526977 2.195668
Infa 10.6757519 2.1956683 85.548053 
Lambda.n
Nata Mort Infa
Nata 957.3249 269.8570 3651.1071
Mort 269.8570 462.6226 750.9186
Infa 3651.1071 750.9186 29257.4343

De los anteriores cálculos, se puede ver que la distribución posterior de \(\boldsymbol \theta\) está dada por \[\begin{equation*} \boldsymbol \theta\mid\mathbf{Y}\sim t_{19}\left(\begin{pmatrix}29.3\\9.2\\54.7\end{pmatrix}, \begin{pmatrix}2.80&0.79&10.68\\0.79&1.35&2.20\\10.68&2.20&85.55\end{pmatrix}\right) \end{equation*}\]

Usando propiedades de la distribución multivariante \(t\) de Student, tenemos que \(\theta_1\sim t_{19}(29.29, 2.80)\), \(\theta_2\sim t_{19}(9.24, 1.35)\) y \(\theta_3\sim t_{19}(0.62, 85.55)\), de allí se puede encontrar fácilmente los intervalos de credibilidad para cada uno de estos tres parámetros.

En cuanto a la distribución posterior de \(\boldsymbol \Sigma\), ésta está dada por \[\begin{equation*} \boldsymbol \Sigma\mid\mathbf{Y}\sim Inversa-Wishart_{21}\left(\begin{pmatrix}957&270&3651\\ 270&463&751\\ 3651&751&29257\end{pmatrix}\right) \end{equation*}\]

La estimación bayesiana de \(\boldsymbol \Sigma\) viene dada por \(\hat{\boldsymbol \Sigma}=\begin{pmatrix}56.3&15.9&214.8\\15.9&27.2&44.2\\214.8&44.2&1721.0\end{pmatrix}\). Por propiedades de la distribución inversa-Wishart, se puede concluir que los elementos diagonales de \(\boldsymbol \Sigma\) tienen distribución inversa-Gamma. Por ejemplo, se tiene que \(\sigma^2_{1}\sim Inversa-Gamma(21/2, 56.3/2)\) y cualquier inferencia que se desear realizar sobre \(\sigma^2_{1}\) es posible a partir de esta distribución.

Aparte de los análisis anteriores, también podemos realizar ejercicios de comparación y verificar la posible independencia entre parejas de variables. Por ejemplo, si queremos verificar la hipótesis de que la tasa de natalidad es dos veces la tasa de mortalidad, esto es \(\theta_1=2 \times \theta_2\). Una forma de confirmar o refutar esta hipótesis es hallar el intervalo de credibilidad de \(\theta_1 - 2\theta_2\) que se puede expresar como \((1,-2,0)\begin{pmatrix}\theta_1\\ \theta_2\\ \theta_3\end{pmatrix}\). Por propiedades de la distribución \(t\) de Student multivariante, tenemos que \(\theta_1 - 2 \theta_2\) tiene distribución \(t\) de Student univariada con los mismos grados de libertad que \(\boldsymbol \theta\). La esperanza de esta distribución está dada por \((1,-2,0)\begin{pmatrix}29.3\\9.2\\54.7\end{pmatrix}=10.8\) y la escala está dada por \((1,-2,0)\begin{pmatrix}2.80&0.79&10.68\\0.79&1.35&2.20\\10.68&2.20&85.55\end{pmatrix}\begin{pmatrix}1\\-2\\0\end{pmatrix}=5.054\), esto es, \[\begin{equation*} \theta_1-2\theta_2 \mid \mathbf{Y}\sim t_{19}(10.8, 5.054) \end{equation*}\]

De esta forma, un intervalo de credibilidad para \(\theta_1-2\theta_2\) viene dado por los percentiles 2.5% y 97.5% de la anterior distribución, que a la vez son iguales a los percentiles 2.5% y 97.5% de la distribución \(t\) de Student estandarizada multiplicado por \(\sqrt{5.054}\) y sumando 10.8. Este intervalo es igual a \((6.095, 15.505)\). Al observar que este intervalo no contiene el valor 0, podemos concluir que no es válido afirmar que la tasa de natalidad sea dos veces la tasa de mortalidad.

En el anterior análisis, vemos que el intervalo de credibilidad para \(\theta_1-2\theta_2\) contiene solo valores positivos, lo cual es un indicio de que la variable \(\theta_1-2\theta_2\) tenga la mayor parte de la función de densidad ubicada en el eje positivo. De hecho podemos indagar por \(Pr(\theta_1-2\theta_2>0)\), la cual se puede calcular de la distribución \(t_{19}(10.8, 5.054)\) encontrada anteriormente. Esta probabilidad es 1 - pt((0 - 10.8)/sqrt(5.054), 19) dando como resultado 0.9999383; de donde se muestra una fuerte evidencia de que la tasa de natalidad es superior a dos veces la tasa de mortalidad.

Los anteriores resultados fueron obtenidos directamente de las distribuciones posteriores marginales \(p(\boldsymbol \theta|\mathbf{Y})\) y \(p(\boldsymbol \Sigma|\mathbf{Y})\). De forma equivalente también se puede usar las técnicas de simulación con base en las distribuciones \(p(\boldsymbol \theta,\boldsymbol \Sigma|\mathbf{Y})\) y \(p(\boldsymbol \Sigma|\mathbf{Y})\). A continuación se muestran los códigos necesario en R:

nsim <- 1000
theta.pos <- matrix(NA, nsim, P)
Sigma.pos <- array(NA, c(nsim, P, P))

for(i in 1:nsim){
  Sigma.pos[i,,] <- riwish(n + v0, Lambda.n)
  theta.pos[i,] <- rmvnorm(1, mu.n, 
                           Sigma.pos[i, , ]/(n + c0))
}

# Estimaciones finales
theta.final <- colMeans(theta.pos)
Sigma.final <- matrix(c(mean(Sigma.pos[, 1, 1]),
                        mean(Sigma.pos[, 1, 2]),
                        mean(Sigma.pos[, 1, 3]),
                        mean(Sigma.pos[, 2, 1]),
                        mean(Sigma.pos[, 2, 2]),
                        mean(Sigma.pos[, 2, 3]),
                        mean(Sigma.pos[, 3, 1]),
                        mean(Sigma.pos[, 3, 2]), 
                        mean(Sigma.pos[, 3, 3])), 
                      nrow = P, ncol = P)
theta.final
## [1] 29.245516  9.288723 54.650509
Sigma.final
55.88544 15.59074 211.03550
15.59074 26.47268 42.52096
211.03550 42.52096 1687.49974

Podemos ver que los resultados obtenidos con los dos métodos son totalmente coincidentes. En cuanto al intervalo de credibilidad para \(\theta_1-2\theta_2\), se puede calcular con

quantile(
  theta.pos[,1] - 2 * theta.pos[, 2], 
  c(0.025, 0.975))
##      2.5%     97.5% 
##  6.120592 15.045185

También podemos calcular \(Pr(\theta_1-2\theta_2>0)\) como

sum(theta.pos[, 1] > 2 * theta.pos[, 2])/nsim
## [1] 1

Podemos ver que estos resultados son muy similares a los obtenidos usando \(p(\boldsymbol \theta|\mathbf{Y})\).

4.3.3 Parámetros no informativos

A. Gelman et al. (2003) afirma que la distribución previa no informativa de Jeffreys conjunta para \(\boldsymbol \theta,\boldsymbol \Sigma\), en este caso está dada por la siguiente expresión \[\begin{equation*} p(\boldsymbol \theta,\boldsymbol \Sigma)\propto \mid \boldsymbol \Sigma\mid ^{-(p+1)/2} \end{equation*}\]

La distribución posterior conjunta para \(\boldsymbol \theta,\boldsymbol \Sigma\) está dada por

\[\begin{equation*} p(\boldsymbol \theta,\boldsymbol \Sigma\mid \mathbf{Y})\propto \mid \boldsymbol \Sigma\mid ^{-(p+n+1)/2} \exp\left\{ -\frac{1}{2}\sum_{i=1}^n (\mathbf{Y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{Y}_i-\boldsymbol \theta)\right\} \end{equation*}\]

De la anterior distribución, podemos encontrar la distribución condicional posterior de \(\boldsymbol \theta\) dada en el siguiente resultado.

Resultado 4.16 La distribución posterior del vector de parámetros \(\boldsymbol \theta\) condicional a \(\boldsymbol \Sigma,\mathbf{Y}\) es \[\begin{equation*} \boldsymbol \theta\mid \boldsymbol \Sigma,\mathbf{Y}\sim Normal_p(\bar{\mathbf{y}},\boldsymbol \Sigma/n) \end{equation*}\]


Prueba. Algunas simples operaciones algebráicas muestran que: \[\begin{align*} p(\boldsymbol \theta\mid \boldsymbol \Sigma,\mathbf{Y}) &\propto \exp\left\{ -\frac{1}{2}\sum_{i=1}^n (\mathbf{Y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{Y}_i-\boldsymbol \theta)\right\}\\ &\propto \exp\left\{ -\frac{n}{2}(\boldsymbol \theta-\bar{\mathbf{Y}})'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\bar{\mathbf{Y}})\right\} \end{align*}\] Por lo tanto, factorizando convenientemente, se encuentra una expresión idéntica a la función de distribución de una variable aleatoria con distribución \(Normal_p(\bar{y},\boldsymbol \Sigma/n)\).


En cuanto a la estimación de \(\boldsymbol \Sigma\), en el siguiente resultado encontramos su distribución posterior.

Resultado 4.17 La distribución marginal posterior de la matriz de parámetros \(\boldsymbol \Sigma\) es \[\begin{equation*} \boldsymbol \Sigma\mid \mathbf{Y}\sim Inversa-Whishart_{n-1}(\mathbf{S}) \end{equation*}\] donde \(\mathbf{S}=\sum_{i=1}^n(\mathbf{y}_i-\bar{\mathbf{y}})(\mathbf{y}_i-\bar{\mathbf{y}})'\)


Prueba. En primer lugar recordamos la expresión \[\begin{equation*} \textbf{S}_{\boldsymbol \theta}=\sum_{i=1}^n(\mathbf{y}_i-\boldsymbol \theta)(\mathbf{y}_i-\boldsymbol \theta)'=\mathbf{S}+n(\boldsymbol \theta-\bar{\mathbf{y}})(\boldsymbol \theta-\bar{\mathbf{y}})' \end{equation*}\]

Por otro lado, recurriendo a las propiedades del operador \(traza\), e integrando la distribución posterior conjunta con respecto a \(\boldsymbol \theta\), se tiene que \[\begin{align*} p(\boldsymbol \Sigma\mid \mathbf{Y})&=\int p(\boldsymbol \theta,\boldsymbol \Sigma\mid \mathbf{Y}) \ d\boldsymbol \theta\\ &= \mid \boldsymbol \Sigma\mid ^{-(p+n+1)/2}\int\exp\left\{ -\frac{1}{2}\sum_{i=1}^n (\mathbf{Y}_i-\boldsymbol \theta)'\boldsymbol \Sigma^{-1}(\mathbf{Y}_i-\boldsymbol \theta)\right\} \ d\boldsymbol \theta\\ &= \mid \boldsymbol \Sigma\mid ^{-(p+n+1)/2}\int \exp\left\{ -\frac{1}{2}traza(\boldsymbol \Sigma^{-1}\mathbf{S}_{\boldsymbol \theta})\right\} \ d\boldsymbol \theta\\ &= \mid \boldsymbol \Sigma\mid ^{-(p+n+1)/2}\int \exp\left\{ -\frac{1}{2}traza(\boldsymbol \Sigma^{-1} (\mathbf{S}+n(\boldsymbol \theta-\bar{\mathbf{y}})(\boldsymbol \theta-\bar{\mathbf{y}})'))\right\} \ d\boldsymbol \theta\\ &= \mid \boldsymbol \Sigma\mid ^{-(p+n)/2}\exp\left\{ -\frac{1}{2}traza(\boldsymbol \Sigma^{-1}\mathbf{S})\right\}\\ &\hspace{2cm}\times \int \mid \boldsymbol \Sigma\mid ^{-1/2}\exp\left\{ -\frac{n}{2}traza(\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\bar{\mathbf{y}})(\boldsymbol \theta-\bar{\mathbf{y}})')\right\} \ d\boldsymbol \theta\\ &= \mid \boldsymbol \Sigma\mid ^{-(p+n)/2}\exp\left\{ -\frac{1}{2}traza(\boldsymbol \Sigma^{-1}\mathbf{S})\right\}\\ &\hspace{2cm}\times\int \mid \boldsymbol \Sigma\mid ^{-1/2}\exp\left\{ -\frac{n}{2}traza((\boldsymbol \theta-\bar{\mathbf{y}})'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\bar{\mathbf{y}}))\right\} \ d\boldsymbol \theta\\ &= \mid \boldsymbol \Sigma\mid ^{-(p+n)/2}\exp\left\{ -\frac{1}{2}traza(\boldsymbol \Sigma^{-1}\mathbf{S})\right\}\\ &\hspace{2cm}\times \int\underbrace{ \mid \boldsymbol \Sigma\mid ^{-1/2}\exp\left\{ -\frac{n}{2}(\boldsymbol \theta-\bar{\mathbf{y}})'\boldsymbol \Sigma^{-1}(\boldsymbol \theta-\bar{\mathbf{y}})\right\}}_{Normal_p(\bar{\mathbf{y}},\boldsymbol \Sigma/n)} \ d\boldsymbol \theta\\ &= \mid \boldsymbol \Sigma\mid ^{-(p+n)/2}\exp\left\{ -\frac{1}{2}traza(\boldsymbol \Sigma^{-1}\mathbf{S})\right\} \end{align*}\]

Por lo tanto, factorizando convenientemente, se encuentra una expresión idéntica a la función de distribución de una variable aleatoria con distribución \(Inversa-Whishart_{n-1}(\mathbf{S})\).


El anterior resultado indica que la estimación bayesiana de \(\boldsymbol \Sigma\) cuando se utiliza una previa no informativa está dada por \[\begin{equation*} \hat{\boldsymbol \Sigma}=\frac{\sum_{i=1}^n(\mathbf{y}_i-\bar{\mathbf{y}})(\mathbf{y}_i-\bar{\mathbf{y}})'}{n-p-2} \end{equation*}\]

Esta expresión es muy similar a la estimación clásica de la matriz de varianzas y covarianzas dada por \(\frac{\sum_{i=1}^n(\mathbf{y}_i-\bar{\mathbf{y}})(\mathbf{y}_i-\bar{\mathbf{y}})'}{n-1}\). Se puede observar que a medida que \(n\) se aumente, las dos expresiones darán resultados muy similares, pero siempre la estimación bayesiana será mayor a la estimación clásica, especialmente en situaciones donde el tamaño muestral es pequeño.

Para obtener la estimación de \(\boldsymbol \Sigma\) junto con la estimación de \(\boldsymbol \theta\), podemos proceder de la siguiente forma para obetener valores simulados de \(\boldsymbol \theta\) y \(\boldsymbol \Sigma\) y así, obtener las estimaciones respectivas. Si el número de iteraciones se fija como \(G\), entonces se procede a:

  1. Simular \(G\) valores de la distribución de \(\boldsymbol \Sigma|\mathbf{Y}\); estos valores se denotan por \(\boldsymbol \Sigma_{(1)},\boldsymbol \Sigma_{(2)},\cdots,\boldsymbol \Sigma_{(G)}\).
  2. Para cada valor de \(\boldsymbol \Sigma_{(g)}\), con \(g=1,\cdots,G\), simular un valor de la distribución de \(\boldsymbol \theta|\boldsymbol \Sigma,\mathbf{Y}\); es decir, de la distribución \(N_p(\bar{\mathbf{y}}, \boldsymbol \Sigma/n)\), donde \(\boldsymbol \Sigma\) se reemplaza por \(\boldsymbol \Sigma_{(g)}\). De esta forma, se obtienen los valores \(\boldsymbol \theta_{(1)},\boldsymbol \theta_{(2)},\cdots,\boldsymbol \theta_{(G)}\).

El siguiente ejemplo ilustra la forma de obtener las estimaciones siguiendo el anterior procedimiento.

Ejemplo 4.5 Retomamos los datos del efecto de aumento en horas de sueño de dos medicamentos soporíferos utilizados en los ejemplos 4.2 y 4.3. Los siguientes códigos en R ilustran el procedimiento computacional para obtener valores de la distribución posterior conjunta de \(\boldsymbol \theta\) y \(\boldsymbol \Sigma\). 
library(MCMCpack)
library(mvtnorm)

y <- as.matrix(
  data.frame(M1 = sleep[1:10, 1], 
             M2 = sleep[-(1:10), 1]))
n <- nrow(y)

y.bar <- colMeans(y)
S <- var(y) * (n - 1)

nsim <- 1000
theta.pos <- matrix(NA, nsim, 2)
Sigma.pos <- array(NA, c(nsim, 2, 2))

for(i in 1:nsim){
  #simulacion de la distribucion posterior condicional de Sigma
  Sigma.pos[i, , ] <- riwish(n - 1, S)
  #simulacion de la distribucion posterior condicional de theta
  theta.pos[i, ] <- rmvnorm(1, y.bar, Sigma.pos[i, , ]/n)
}

Dado que en el cálculo no se hizo uso de valores iniciales y por la forma de las distribuciones posteriores de \(p(\boldsymbol \theta|\boldsymbol \Sigma,\mathbf{Y})\) y \(\boldsymbol \Sigma|\mathbf{Y}\), los valores muestrados en la diferentes iteraciones no guardan relación entre sí; por ende, podemos usar directamente todos los valores simulados para realizar el cálculo de las estimaciones bayesianas.

theta.Bayes <- colMeans(theta.pos)
Sigma.Bayes <- matrix(c(mean(Sigma.pos[, 1, 1]),
                        mean(Sigma.pos[, 2, 1]),
                        mean(Sigma.pos[, 1, 2]),
                        mean(Sigma.pos[, 2, 2])), 
                      ncol = 2, nrow = 2)
theta.Bayes
## [1] 0.7890492 2.3443686
Sigma.Bayes
4.985367 4.337360
4.337360 6.002711

Por otro lado, la estimación clásica de los parámetros está dada por

y.bar
##   M1   M2 
## 0.75 2.33
var(y)
M1 M2
M1 3.200556 2.848333
M2 2.848333 4.009000

Por consiguiente, podemos observar que, en cuanto al parámetro \(\boldsymbol \theta\), la estimación bayesiana es igual a la estimación clásica; mientras que el determinante de la estimación bayesiana de \(\boldsymbol \Sigma\) es mucho mayor que el de la estimación clásica, esto ocurre en situaciones cuando el tamaño muestral es pequeño. En cuanto a la estimación por intervalo de los efectos promedios de los dos medicamentos, tenemos que:

quantile(theta.pos[, 1], c(0.025, 0.975))
##       2.5%      97.5% 
## -0.6388661  2.1506794
t.test(y[, 1])$conf.int
## [1] -0.5297804  2.0297804
## attr(,"conf.level")
## [1] 0.95
quantile(theta.pos[, 2], c(0.025, 0.975))
##      2.5%     97.5% 
## 0.7742381 3.9578339
t.test(y[,2])$conf.int
## [1] 0.8976775 3.7623225
## attr(,"conf.level")
## [1] 0.95

De las anteriores salidas observamos que los resultados obtenidos con el enfoque bayesiano, aunque no son exactamente iguales a los obtenidos con el enfoque clásico, sí son muy similares. En cuanto a la estimación por intervalo de las varianzas y covarianzas. Primero consideramos la varianza del primer medicamento denotada por \(\sigma^2_1\). La distribución posterior de la matriz de varianzas y covarianzas está dada por \[\begin{equation*} \boldsymbol \Sigma=\begin{pmatrix}\sigma^2_1&\sigma_{12}\\\sigma_{21}&\sigma^2_{2} \end{pmatrix}\sim Inversa-Wishart_9(\mathbf{S}) \end{equation*}\]

con \(\mathbf{S}=\sum_{i=1}^{10}(\mathbf{y}_i-\bar{\mathbf{y}})(\mathbf{y}_i-\bar{\mathbf{y}})'=\begin{pmatrix}28.81&25.64\\25.64&36.08\end{pmatrix}\). Usando las propiedades de la distribución Inversa-Wishart, se puede concluir que la distribución marginal posterior de \(\sigma^2_1\) está dada por \(Inversa-Gamma(\alpha=\frac{9-1}{2}, \beta=\frac{28.81}{2})\), y su intervalo de credibilidad se puede calcular directamente de dicha distribución, o equivalentemente usando los percentiles muestrales de los valores de \(\sigma^2_1\) simulados El intervalo obtenido por estos dos medios son muy similares como se puede ver a continuación.

library(pscl)
qigamma(0.025, alpha = 8/2, beta = 28.81/2)
## [1] 1.643042
qigamma(0.975, alpha = 8/2, beta = 28.81/2)
## [1] 13.21723
quantile(Sigma.pos[, 1, 1], c(0.025, 0.975))
##     2.5%    97.5% 
##  1.60053 14.00017

El intervalo de confianza del 95% se puede obtener con el siguiente código10.

c(9 * var(y[, 1]) / qchisq(0.975, 9), 
  9 * var(y[, 2]) / qchisq(0.025, 9))
## [1]  1.514238 13.361406

En comparación con el intervalo de credibilidad, el intervalo de confianza está ubicado levemente hacia la izquierda del eje real, esto se debe a que la estimación clásica de la varianza siempre será menor a la estimación bayesiana con una previa no informativa.

Referencias

———. 2003. Bayesian Data Analysis. 2.ª ed. Chapman; Hall/CRC.
Peña, D. 2002. Análisis de datos multivariantes. McGraw-Hill.
Student. 1908. «The Probable Error of a Mean». Biometrika 6 (1): 1-25.
Zhang, H., y H. A. Gutiérrez. 2010. Teoría estadística. Aplicación y métodos. Universidad Santo Tomás.

  1. Consultar Zhang y Gutiérrez (2010, sec.3.2.1) para mayor información.↩︎