Modo de la distribución inversa de Wishart (muestreada frente a calculada)

Question

Recientemente he investigado el comportamiento muestral de las matrices de covarianza mediante simulación. Me di cuenta de que el modo de las matrices distribuidas inversas de Wishart simuladas difiere en cierto modo del modo "teórico" que puede calcularse directamente a partir de los parámetros de la distribución.

Supongamos que estoy interesado en una matriz $\Sigma$ de tamaño $p \times p$ con entradas conocidas. Denoto la distribución inversa de Wishart $W^{-1}(v,S)$ donde $v$ son los grados de libertad, y $S$ es la matriz de escala de tamaño $p \times p$ . Ahora, entiendo que el modo de la inversa de Wishart es

$$ \frac{S}{v+p+1} \; . $$

Por lo tanto, la simulación a partir de $W^{-1}(v,(v+p+1)\Sigma)$ debe dan como resultado una distribución de matrices con una moda que coincide con las entradas de $\Sigma$ . Para mi sorpresa, no fue así.

A continuación se muestra el código que he utilizado para generar las matrices en R.

sim.invWishart <- function(N, v, Scale){

  invSc <- solve(Scale)
  p <- ncol(Scale)

  invS <- rWishart(N,v,invSc)
  S <- apply(invS, 3, solve)
  S <- array(S,dim=c(p,p,N))

  return(S)
}

df <- 10
Sigma <- matrix(c(.8,.3,.1,.3,.8,.1,.1,.1,.4),3,3)
Scale <- Sigma*(df+ncol(Sigma)+1)

N <- 50000
S <- sim.invWishart(N, df, Scale)

# mode of entries
apply(S, 1:2, function(x){d<-density(x); d$x[which.max(d$y)]})

Esperaba que el modo de las entradas replicara Sigma . Sin embargo, la última llamada (modo a través de las denisidades del núcleo) da como resultado la siguiente matriz, que contiene el modo de cada entrada. También imprimí un histograma para mostrar mejor lo que quiero decir.

          [,1]      [,2]      [,3]
[1,] 1.1364036 0.3459420 0.1369788
[2,] 0.3459420 1.1941119 0.1168395
[3,] 0.1369788 0.1168395 0.5587053

En realidad, las entradas fuera de diagonal cumplieron bastante bien mis expectativas. Sin embargo, el modo de las entradas diagonales siempre está un poco desviado. Sorprendentemente, al simular a partir de $W^{-1}(v,v\Sigma)$ en cambio, los modos obtenidos de la simulación coinciden exactamente con $\Sigma$ ) aunque la fórmula del modo sugiera lo contrario.

¿Por qué? ¿Es mi procedimiento defectuoso (por ejemplo, la simulación, el cálculo de los modos), o es este comportamiento común de la inversa-Wishart? ¿Por qué la multiplicación por $v+p+1$ dan resultados erróneos, y por qué son (aparentemente) correctos con sólo $v$ en ese lugar?

Gracias de antemano.

Answer 1

Creo que es porque los histogramas muestran las distribuciones marginales, por tanto la moda marginal de cada componente de la matriz. La moda teórica es la moda conjunta, considerando todos los componentes juntos. En otras palabras, la mayor probabilidad no se consigue optmizando individualmente las probabilidades marginales. Lamentablemente, no conozco una forma de representar gráficamente la moda de la distribución. Pero, utilizando técnicas de optimización, deberías poder llegar al óptimo correcto.