Maximización de la probabilidad frente al muestreo MCMC: Comparación de parámetros y desviaciones

Question

Estoy trabajando en R. Utilizo lm() para maximizar la probabilidad en el primer análisis, y STAN para muestrear de la posterior en un segundo análisis.

require(rstan)

He fabricado algunos datos.

set.seed(123)

N      <- 1000
data   <- as.data.frame(sapply(1:3,function(x)rnorm(N)))
data$y <- 2 + data$V1 + 2*data$V2 + 4*data$V3 + rnorm(N,0,4)

He ajustado un modelo lineal a los datos. Más adelante presentaré los resultados.

lm1 <- lm(y ~ V1 + V2 + V3,data = data)

A continuación, utilizo STAN para muestrear de la distribución posterior de un modelo muy similar (el mismo modelo, pero STAN aplica automáticamente priores uniformes en los parámetros donde no se especifica un prior). A continuación se muestra el modelo de STAN, que se guarda en un archivo llamado LM.stan .

data{
    int N;
    real y[N];
    matrix[N,3] X;
}
parameters{
    real a;
    vector[3] b;
    real sigma;
}
model{
   y ~ normal( a + X*b, sigma );
}
generated quantities{
    real logLik;
    vector[N] lm;
    lm <- a + X*b;
    for( i in 1:N ) {
        logLik <- logLik + normal_log( y[i], lm[i], sigma );
    }
}

A continuación se muestra el código R que crea un objeto de datos con el que STAN puede lidiar, y compila y ejecuta el muestreador.

standata <- list(X = model.matrix(~-1+V1+V2+V3,data=data),
                 N = nrow(data),
                 y = data$y
            )
lm2 <- stan("LM.stan",verbose=F,iter=1000,data=standata,refresh=100)

Podemos examinar las medias y las desviaciones estándar de las distribuciones marginales de cada parámetro, y la probabilidad logarítmica media.

samples    <- do.call(cbind,extract(lm2,c("a","b","sigma","logLik")))
colnames(samples) <- c('a','b1','b2','b3','sig','logLik')
fitTab     <- data.frame(t(apply(samples,2,function(x) c(mean=mean(x),sd=sd(x)))))
fitTab
#                mean        sd
#a          1.9671178 0.1272709
#b1         0.9951333 0.1290366
#b2         1.9647771 0.1263366
#b3         4.2091996 0.1339271
#sig        3.9764835 0.0842643
#logLik -2798.2644784 1.6468066

Y para comparar, imprima las estimaciones de los parámetros de máxima verosimilitud y la máxima verosimilitud logarítmica encontrada usando lm()

summary(lm1)

#Call:
#lm(formula = y ~ V1 + V2 + V3, data = data)
#
#Residuals:
#     Min       1Q   Median       3Q      Max 
#-12.2247  -2.5338   0.0019   2.7010  11.5690 
#
#Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   1.9690     0.1257  15.664  < 2e-16 ***
#V1            0.9948     0.1272   7.823 1.31e-14 ***
#V2            1.9675     0.1249  15.750  < 2e-16 ***
#V3            4.2059     0.1285  32.740  < 2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 3.97 on 996 degrees of freedom
#Multiple R-squared:  0.5884,   Adjusted R-squared:  0.5872 
#F-statistic: 474.7 on 3 and 996 DF,  p-value: < 2.2e-16

logLik(lm1)
#'log Lik.' -2795.739 (df=5)

Observo que las estimaciones de los parámetros y los errores estándar son muy similares a las medias y las DE de las muestras posteriores. Sin embargo, las probabilidades logarítmicas son diferentes. Reconozco que la verosimilitud reportada para las muestras posteriores es una media, así que compruebo cuál es el mayor valor de verosimilitud de las muestras y cuáles son los valores de los parámetros asociados a él. Efectivamente, el mayor valor de verosimilitud es muy similar a la máxima verosimilitud de los datos dado el modelo OLS.

t(samples[samples[,'logLik'] == max(samples[,'logLik']), , drop=F])
#                [,1]
#a          1.9935327
#b1         0.9750932
#b2         1.9817911
#b3         4.2507953
#sig        3.9632800
#logLik -2795.8380670

Sin embargo, los valores de los parámetros que dieron lugar a la mayor probabilidad de las muestras difieren bastante de las estimaciones OLS.

Para resumir, los dos hechos extraños y probablemente relacionados que observo son que 1) las medias marginales de los parámetros se parecen mucho a las estimaciones de MCO, pero la media marginal de log-verosimilitud no se parece a la de máxima verosimilitud (de MCO), y 2) los parámetros asociados a la muestra de mayor verosimilitud (de STAN), cuyo valor es muy cercano al de la verosimilitud de MCO, no coinciden con los de las estimaciones de máxima verosimilitud (de MCO).

¿Por qué no debería sorprenderme esto (especialmente el hecho 2)?

Answer 1

Imagina que tu posterior es algo gaussiano. La expectativa de la norma euclidiana al cuadrado de la mejor de $N$ dibujar depende de la dimensión $d$ y es asintóticamente $O(1/N^{2/d})$ pero la media sigue siendo $O(1/N)$ . En cuanto tienes más de 2 dimensiones, la media converge más rápido. Se tiene $d=4$ por lo que la media de los puntos alrededor de la moda será una estimación mucho mejor de la moda que el punto más cercano a la misma.

Si eso no es inmediatamente intuitivo, imagine que está extrayendo de una normal multivariable estándar con covarianza de identidad y $d=1000$ . El vector más cercano a $0$ seguirá estando muy lejos en promedio; en dimensiones altas, la mayor parte de la masa de una gaussiana está lejos del centro. La media estará mucho más cerca, ya que los valores de las diferentes extracciones se anulan en cada dimensión.

Pruébalo en python

def draw(d):
   x = np.random.randn(1000*d).reshape((1000,d))
   return np.sum(np.mean(x,axis=0)**2) < np.min(np.sum(x**2,axis=1))

for d in range(1,4):
   print d, np.mean( [draw(d) for i in range(0,1000)] )

Edición: así que me puse a calcular el factor multiplicativo en la expectativa de la norma euclidiana mínima al cuadrado cuando se extrae de una normal multivariada estándar con dimensión $d$ . Es asintóticamente equivalente a

$$2\Gamma\left(1+\frac{d}{2}\right)^{\frac{2}{d}}\Gamma\left(1+\frac{2}{d}\right) n^{-2/d}$$

mientras que la media es, por supuesto, de $d n^{-1}$