¿Cómo puedo prácticamente "marginar un parámetro molesto"?

Question

Genero mis datos con este modelo:

$y=ax+b+\nu$

donde $\nu$ es un valor aleatorio de una variable aleatoria que sigue una distribución normal con media igual a cero y desviación típica igual a $\sigma$ . El siguiente script R genera mis datos.

set.seed(42)
a_t=1 #true value for a
b_t=1 #true value for b
sigma_t=1 #true value for sigma
x=1:10
noise=rnorm(length(x),0,sigma_t)
y=a_t*x+b_t+noise
plot(x,y)
abline(a_t,b_t)

Ahora me gustaría "recuperarme" de $x$ y $y$ los valores de $a$ , $b$ y $\sigma$ utilizando un enfoque bayesiano. Seguiré el aproximación de malla 1,2

1. Define la cuadrícula. Esto significa que se decide cuántos puntos se van a utilizar en la estimación de la posterior, y luego se hace una lista del parámetro en la cuadrícula.
2. Calcular el valor de la prioridad en cada valor del parámetro en la cuadrícula.
3. Calcule la probabilidad para cada valor de los parámetros.
4. Calcular la posterior no estandarizada en cada valor del parámetro, multiplicando la anterior por la verosimilitud. En este paso, estamos la densidad posterior por una distribución de probabilidad discreta discreta en la cuadrícula.
5. Estandarizar la posterior, dividiendo cada valor por la suma de todos los valores. Por último, utilizando el comando de R sample tomar una muestra aleatoria con reemplazamiento de la distribución discreta.

# define the grid
a=seq(0,2,by=.05)
b=seq(0,2,by=.05)
sigma=seq(.1,2,by=.01)

pars <- expand.grid(a=a,b=b,sigma=sigma)

# define the priors
pars$a_prior <- dunif(pars$a,min=0,max=2)
pars$b_prior <- dunif(pars$b,min=0,max=2)
pars$sigma_prior <- dunif(pars$sigma,min=0,max=2)
pars$prior <- pars$a_prior * pars$b_prior * pars$sigma_prior

# compute the likelihood
sz <- nrow(pars)
for(i in 1:sz){
    likelihoods <- dnorm(y,pars$a[i]*x+pars$b[i],pars$sigma[i])
    pars$likelihood[i] <- prod(likelihoods)
    print(i/sz)
}

# compute the unnormalized posterior
pars$unnormalized_posterior <- pars$likelihood*pars$prior
# normalize the posterior
pars$posterior <- pars$unnormalized_posterior/sum(pars$unnormalized_posterior)

# get the sample from the posterior
sample_indices <- sample(1:nrow(pars), size=1e4, replace = TRUE, prob=pars$posterior)
pars_sample <- pars[sample_indices,c('a','b','sigma')]

# draw a 2d histogram for a and b
library(ggplot2)
p <- ggplot(pars_sample,aes(a,b))
p <- p+stat_bin2d(bins=10)
plot(p)

Al final del procedimiento tengo una matriz de tres columnas (o el pars_sample dataframe en el script R) que es una muestra de la distribución posterior $p(a,b,\sigma|x,y)$ ; ahora me gustaría pensar en $\sigma$ como un parámetro molesto (¿tiene sentido?) y me gustaría "marginarlo", es decir, me gustaría obtener

$p(a,b|x,y)=\int p(a,b,\sigma|x,y) d \sigma$

Hogg et al. 3 escribió:

Esta marginación puede realizarse por integración numérica directa (piensa en cuadricular los parámetros molestos, evaluar la probabilidad posterior en cada punto de la cuadrícula y sumar). en cada punto de la cuadrícula, y sumando)

pero no tengo ni idea de cómo proceder en R.

1 MCELREATH, Richard. Replanteamiento estadístico: Un curso bayesiano con ejemplos en R y Stan. CRC press , 2020. Página 40.

2 ALBERT, Jim. Cálculo bayesiano con R. Nueva York, NY : Springer Nueva York, 2009. Página 27.

3 HOGG, David W., BOVY, Jo y LANG, Dustin. Recetas de análisis de datos: Ajuste de un modelo a los datos. arXiv:1008.4686

Answer 1

Es tan sencillo como soltar el $\sigma$ columna de su pars_sample dataframe.

El marco de datos de tres columnas es una muestra de la distribución conjunta de los tres parámetros. Dejando $\sigma$ deja una muestra de la distribución conjunta de los dos parámetros restantes. El procedimiento consiste en generar a, b y $\sigma$ teniendo en cuenta la influencia de los tres valores, pero luego ignorar $\sigma$ para calcular la distribución marginal de a y b.

Esto es diferente del procedimiento que usted cita de Hogg et al. Su método es una modificación de su paso anterior. En lugar de muestrear a partir de la posterior conjunta completa, usted reduciría la posterior conjunta de tres vías a una posterior marginal de dos vías (a y b) y muestrearía a partir de ella.

Los dos métodos son válidos. Deberías elegir el que sea más conveniente desde el punto de vista computacional en términos de precisión deseada, tiempo de cálculo, complejidad, etc. Desde tu punto de vista, parece que lo más sencillo es generar muestras de los tres parámetros, soltar el parámetro $\sigma$ y luego calcular la distribución marginal a y b.