Inferencia bayesiana en el parámetro de correlación de un bivariante normal

Question

Supongamos que los datos de $y_i$'s de los siguientes normal bivariante

$$y_i\sim \mathcal{N}\bigg(\mu\left[ {\begin{array}{cc} \sigma_{11} & \sqrt{\sigma_{11}\sigma_{22}}\rho \\ \sqrt{\sigma_{11}\sigma_{22}}\rho & \sigma_{22} \\ \end{array} } \right]\bigg).$$

Supongamos que $\mu$, $\sigma_{11}$ y $\sigma_{22}$ son de todos conocidos y uno que quiere aprender de la distribución posterior de los $\rho$ bajo algunos antes de la distribución, por ejemplo,

$$\dfrac{\rho+1}{2}\sim beta(2,2).$$

Mi pregunta es, ¿puede la parte posterior directamente muestreado? Es allí cualquier conjugado antes de que puede resultar en algunos manejable posterior?

He trabajado a través de la tediosa tarea de matemáticas y tienen las siguientes

$$L(y_1,\ldots,y_n|\rho)\propto(1-\rho^2)^{-\frac{n}{2}}\exp\bigg\{-\dfrac{\sum_{i=1}^{n}\tilde{y}_{i1}^2 - 2\rho\tilde{y}_{i1}\tilde{y}_{i2}+\tilde{y}_{i2}^2}{2(1-\rho^2)}\bigg \},$$ donde$\tilde{y}_{i1} = (y_{i1}-\mu_1)/\sqrt{\sigma_{11}}$$\tilde{y}_{i2} = (y_{i2}-\mu_2)/\sqrt{\sigma_{22}}$. Sin embargo, esto no me recuerdan a las de cualquier posible conjugar antes.

O, si no hay conjugado previa disponible, cualquiera sugieren un buen rechazo sampler estrategia? Lo que podría eficiente de rechazo de la propuesta de distribución?

Alguna sugerencia?

Answer 1

Desde $$L(y_1,\ldots,y_n|\rho)\propto(1-\rho^2)^{-\frac{n}{2}}\exp\bigg\{-\dfrac{\sum_{i=1}^{n}\tilde{y}_{i1}^2 - 2\rho\tilde{y}_{i1}\tilde{y}_{i2}+\tilde{y}_{i2}^2}{2(1-\rho^2)}\bigg \}$$ es una función de $\rho$ de la forma $$(1-\rho^2)^{-\alpha}\exp\bigg\{-\dfrac{\beta}{1-\rho^2}-\dfrac{\gamma\rho}{1-\rho^2}\bigg \}\qquad (1)$$this leads to an exponential family choice of a conjugate prior (with $\alpha,\beta>0$ and $|\gamma|<\beta$). While this is not a standard distribution, as far as I know, an accept-reject solution may be available, using the bound$$\dfrac{\sum_{i=1}^{n}\tilde{y}_{i1}^2 - 2\rho\tilde{y}_{i1}\tilde{y}_{i2}+\tilde{y}_{i2}^2}{2(1-\rho^2)}\ge\dfrac{\sum_{i=1}^{n}\min[(\tilde{y}_{i1}-\tilde{y}_{i2})^2,(\tilde{y}_{i1}+\tilde{y}_{i2})^2]}{2(1-\rho^2)}$$ puede ayudar, aunque no he sido capaz de encontrar una manera sencilla de simular de $$f(\rho)\propto(1-\rho^2)^{-\alpha}\exp\bigg\{-\dfrac{\beta}{1-\rho^2}\bigg\}$$ Obviamente, puesto que el anterior densidad (1) es acotado, un ataque de fuerza bruta aceptar-rechazar método basado en un Uniforme de las obras, si bien lentamente [en la imagen debajo de la tasa de aceptación es de 2.35%!]:

targ=function(x,a,b,c){
  (1-x*x)^{-a}*exp(-b/(1-x*x)-c*x/(1-x*x))}

upb=function(a,b,c){
  return(optimise(targ,maximum=TRUE,a=a,b=b,c=c,inte=c(-1,1))$obj)}

simz=function(n,a=1,b=1,c=0){
  bon=upb(a,b,c)
  rejcz=integrate(targ,low=-1,upp=1,a=a,b=b,c=c)$val/2/bon
  uniz=runif(ceiling(2*n/rejcz),min=-1,max=1)
  vuniz=runif(ceiling(2*n/rejcz))
  samplz=uniz[vuniz<targ(uniz,a,b,c)/bon]
  return(samplz[1:n])}

Answer 2

Parece que una de Laplace aproximación funciona bastante bien. A continuación se define el logaritmo de la probabilidad y su gradiente. Tenga en cuenta que tengo que cambiar la variable, por lo que el apoyo está en la recta real para una mejor aproximación de Laplace de rendimiento. Yo uso una transformación logit, es decir, $\rho = \dfrac{2}{e^{-x}+1}-1$.

likfcn <- function(x, a, b, c, log = FALSE) {
  rho = 2 / (exp(-x) + 1) - 1
  aux = 1 - rho^2
  llik = -a * log(aux) - b / aux - c * rho / aux
  if (log) {
    return(llik)
  } else {
    return(exp(llik))
  }
}

llikGradient <- function(x, a, b, c) {
  rho = 2 / (exp(-x) + 1) - 1
  aux = 1 - rho^2
  llik_gradient_rho = (2 * a * rho * aux - 2 * b * rho - c - c * rho^2) / aux^2
  return(llik_gradient_rho * 2 * exp(x) / (1 + exp(x))^2)
}

Entonces me simular algunos datos de algunos de los sencillos bivariante de gauss y calcular el correspondiente a, b y c.

# simulation data
set.seed(2017)
suppressMessages(require(mvtnorm))
Sigma = matrix(c(1, .5, .5, 1), 2, 2)
n = 20
y = rmvnorm(n, c(0, 0), Sigma)
a = n / 2
b = sum(y[, 1]^2 + y[, 2]^2) / 2
c = -sum(y[, 1] * y[, 2])

Entonces me maximizar el logaritmo de la probabilidad y obtener el estado de hesse en el máximo. Con la información que se puede crear una aproximación de Laplace.

fn <- function(x) {
  -likfcn(x, a, b, c, log = TRUE)
}
gr <- function(x) {
  -llikGradient(x, a, b, c)
}

optim_res = optim(par = 0, fn = fn, gr = gr, method = "BFGS", lower = -Inf, upper = Inf,hessian = TRUE)

# laplace approximation
laplace_mean = optim_res$par
laplace_var = 1 / optim_res$hessian

Entonces puedo comparar la aproximación de Laplace a la verdadera posterior.

# compare the laplace approximation to the true posterior
x_rho = seq(-.99, .99, .01)
targ=function(x,a,b,c){
  (1-x*x)^{-a}*exp(-b/(1-x*x)-c*x/(1-x*x))}
normalizing_const = integrate(targ,a,b,c,low=-1,upp=1)$val
y_dens_true = targ(x_rho,a,b,c) / normalizing_const
plot(x_rho, y_dens_true, 'l')
x = log((x_rho+1)/2/(1-(x_rho+1)/2))
y_dens_laplace = dnorm(x, laplace_mean, sqrt(laplace_var)) * (1 / (1+x_rho) + 1 / (1-x_rho))
lines(x_rho, y_dens_laplace, col = "red")

donde la línea negra es la densidad de la trama de la verdad y la roja es la aproximación de Laplace. Como puede verse, Laplace aproximación funciona bien, especialmente tomando nota de que el tamaño de los datos de $n=20$ no es grande.