Stan $\hat{R}$ contra Gelman-Rubin $\hat{R}$ definición

Question

Estaba revisando la documentación de Stan que puede descargarse de aquí . Me interesó especialmente su aplicación del diagnóstico Gelman-Rubin. El documento original Gelman y Rubin (1992) definen el factor potencial de reducción de escala (PSRF) de la siguiente manera:

Sea $X_{i,1}, \dots , X_{i,N}$ sea el $i$ ª cadena de Markov muestreada, y que haya globalmente $M$ cadenas independientes muestreadas. Sea $\bar{X}_{i\cdot}$ sea la media de los $i$ ª cadena, y $\bar{X}_{\cdot \cdot}$ sea la media global. Definir, $$W = \dfrac{1}{M} \sum_{m=1}^{M} {s^2_m}, $$ donde $$s^2_m = \dfrac{1}{N-1} \sum_{t=1}^{N} (X_{m t} - \bar{X}_{m \cdot})^2\,. $$ Y definir $B$ $$B = \dfrac{N}{M-1} \sum_{m=1}^{M} (\bar{X}_{m \cdot} - \bar{X}_{\cdot \cdot})^2 \,.$$

Defina $$\hat{V} = \left(\dfrac{N-1}{N} \right)W + \left( \dfrac{M+1}{MN} \right)B\,.$$ El PSRF se estima con $\sqrt{\hat{R}}$ donde $$ \hat{R} = \dfrac{\hat{V}}{W} \cdot \dfrac{df+3}{df+1}\,,$$ donde $df = 2\hat{V}/Var(\hat{V})$ .

La documentación de Stan en la página 349 ignora el término con $df$ y también elimina el $(M+1)/M$ término multiplicativo. Esta es su fórmula,

El estimador de la varianza es $$\widehat{\text{var}}^{+}(\theta \, | \, y) = \frac{N-1}{N} W + \frac{1}{N} B\,.$$ Por último, la estadística de reducción de escala potencial se define por $$ \hat{R} = \sqrt{\frac{\widehat{\text{var}}^{+}(\theta \, | \, y) }{W}}\,. $$

Por lo que he podido ver, no proporcionan una referencia para este cambio de fórmula, y tampoco lo discuten. Normalmente $M$ no es demasiado grande, y a menudo puede ser tan bajo como $2$ Así que $(M+1)/M$ no debe ignorarse, aunque el $df$ puede aproximarse con 1.

¿De dónde procede esta fórmula?

---

EDITAR: He encontrado una respuesta parcial a la pregunta " ¿de dónde procede esta fórmula? ", ya que el Análisis bayesiano de datos, de Gelman, Carlin, Stern y Rubin (Segunda edición) tiene exactamente la misma fórmula. Sin embargo, el libro no explica cómo/por qué está justificado ignorar esos términos?

Answer 1

He seguido el enlace específico dado para Gelman & Rubin (1992) y tiene $$ \hat{\sigma} = \frac{n-1}{n}W+ \frac{1}{n}B $$ como en las versiones posteriores, aunque $\hat{\sigma}$ sustituido por $\hat{\sigma}_+$ en Brooks & Gelman (1998) y con $\widehat{\rm var}^+$ en BDA2 (Gelman et al, 2003) y BDA3 (Gelman et al, 2013).

BDA2 y BDA3 (no pude comprobar ahora BDA1) tienen un ejercicio con pistas para demostrar que $\widehat{\rm var}^+$ es una estimación insesgada de la cantidad deseada.

Gelman & Brooks (1998) tienen la ecuación 1.1 $$ \hat{R} = \frac{m+1}{m}\frac{\hat{\sigma}_+}{W} - \frac{n-1}{mn}, $$ que puede reordenarse como $$ \hat{R} = \frac{\hat{\sigma}_+}{W} + \frac{\hat{\sigma}_+}{Wm}- \frac{n-1}{mn}. $$ Podemos ver que el efecto del segundo y tercer término son insignificantes para la toma de decisiones cuando $n$ es grande. Véase también la discusión en el párrafo anterior a la Sección 3.1 en Brooks & Gelman (1998).

Gelman & Rubin (1992) también tenían el término con df como df/(df-2). Brooks & Gelman (1998) tienen una sección que describe por qué esta corrección de df es incorrecta y definen (df+3)/(df+1). En el párrafo anterior a la sección 3.1 de Brooks & Gelman (1998) se explica por qué se puede prescindir de (d+3)/(d+1).

Parece que tu fuente para las ecuaciones era algo posterior a Brooks & Gelman (1998) ya que allí tenías (d+3)/(d+1) y Gelman & Rubin (1992) tenían df/df(-2). Por lo demás, Gelman y Rubin (1992) y Brooks y Gelman (1998) tienen ecuaciones equivalentes (con notaciones ligeramente diferentes y algunos términos dispuestos de forma distinta). BDA2 (Gelman, et al., 2003) ya no tiene términos $\frac{\hat{\sigma}_+}{Wm}- \frac{n-1}{mn}$ . BDA3 (Gelman et al., 2003) y la versión de cadenas divididas introducida por Stan.

Mi interpretación de los documentos y experiencias utilizando diferentes versiones de $\hat{R}$ es que los términos que finalmente se han suprimido pueden ignorarse cuando $n$ es grande, incluso cuando $m$ no lo es. También recuerdo vagamente haber discutido esto con Andrew Gelman hace años, pero si quieres estar seguro de la historia, deberías preguntarle a él.

Normalmente M no es demasiado grande, y a menudo puede ser tan bajo como 2

Espero de verdad que esto no ocurra a menudo. En los casos en que desee utilizar split- $\hat{R}$ diagnóstico de convergencia, debe utilizar al menos 4 cadenas divididas y así tener M=8. Puede utilizar menos cadenas, si ya sabe que en sus casos específicos la convergencia y la mezcla son rápidas.

Referencia adicional:

• Brooks y Gelman (1998). Journal of Computational and Graphical Statistics, 7(4)434-455.