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?
Inferencia bayesiana en el parámetro de correlación de un bivariante normal
- Preguntado el 2 de Octubre, 2017
- Cuando se hizo la pregunta
- 80 visitas
- Cuantas visitas ha tenido la pregunta
- 2 Respuestas
- Cuantas respuestas ha tenido la pregunta
- Solucionado
- Estado actual de la pregunta
Respuestas
¿Demasiados anuncios?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])}
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.
Preguntas relacionadas
- ¿Cuál es la diferencia entre un filtro de partículas (secuenciales de Monte Carlo) y un filtro de Kalman?
- ¿Por qué la Repisa de la chimenea de la prueba de preferencia sobre Moran yo?
- La corrección de los valores de p para pruebas múltiples, donde las pruebas son correlacionados (genética)
- Para asegurar la coherencia y la comparabilidad, se han preparado dos series cronológicas:
- La automatización de correlación estadística entre los "textos" y "datos"
- ¿Por qué es un Bayesiano no se permite mirar los residuos?
- Basada en la entropía de la refutación de Shalizi del Bayesiano hacia atrás de la flecha del tiempo de la paradoja?
- ¿Qué correlación hace una matriz singular y cuáles son las implicaciones de la singularidad o cerca de la singularidad?
- ¿Cuando (si) es un enfoque frecuentista sustancialmente mejor que un Bayesiano?
- ¿Causalidad implica correlación?
- De correlación intraclase y agregación
Preguntas Destacadas
- Una conjetura de paralelogramo dentro convexo central y simétrica de la curva de
- Encontrar un discontinuo subvariedad
- ¿Por qué es axi autobús populares para la fpga?
- PMF de estimación: la concentración de las desigualdades de la $l_1$ $l_\infty$ errores
- De referencia de la solicitud: Morita contextos
- Van dos cargas idénticas que se mueve a la misma velocidad que la experiencia de fuerza del campo magnético debido a cada uno de los otros?
En nuestra red
- En C#, para que una clase herede de otra clase y una interfaz?
- No puede realizar cambios a la carpeta específica
- Hay alguna forma de establecer automáticamente saliente dirección de correo electrónico basado en el correo electrónico del destinatario o de dominio de correo.aplicación?
- Estoy viajando desde Heathrow a Toledo en los Estados Unidos con 2 conexiones, necesito ver a mi equipaje antes de llegar a Toledo?
- ¿Cómo obtener ejecutable de dumpcap?