Stan

Eu estava revisando a documentação do Stan, que pode ser baixada aqui . Eu estava particularmente interessado na implementação do diagnóstico Gelman-Rubin. O artigo original Gelman & Rubin (1992) define o potencial fator de redução de escala (PSRF) da seguinte maneira:

Deixe que $X_{i,1}, \dots , X_{i,N}$ ser o $i$ th cadeia de Markov amostrado, e que não haja geral $M$ cadeias independentes amostrados. Deixe $\bar{X}_{i\cdot}$ ser a média do $i$ th da cadeia, e $\bar{X}_{\cdot \cdot}$ ser o global significativo. Definir,

W = \frac{1}{M} \sum_{m = 1}^{M} s_{m}^{2},

$W = \dfrac{1}{M} \sum_{m=1}^{M} {s^2_m},$ onde

E defina

s_{m}^{2} = \frac{1}{N - 1} \sum_{t = 1}^{N} ({\bar{X}}_{m t} - {\bar{X}}_{m \cdot})^{2} .

$s^2_m = \dfrac{1}{N-1} \sum_{t=1}^{N} (\bar{X}_{m t} - \bar{X}_{m \cdot})^2\,.$

B

$B$

B = \frac{N}{M - 1} \sum_{m = 1}^{M} ({\bar{X}}_{m \cdot} - {\bar{X}}_{\cdot \cdot})^{2} .

$B = \dfrac{N}{M-1} \sum_{m=1}^{M} (\bar{X}_{m \cdot} - \bar{X}_{\cdot \cdot})^2 \,.$

Definir O PSRF é estimado com

\hat{V} = (\frac{N - 1}{N}) W + (\frac{M + 1}{M N}) B .

$\hat{V} = \left(\dfrac{N-1}{N} \right)W + \left( \dfrac{M+1}{MN} \right)B\,.$

em que

\sqrt{\hat{R}}

$\sqrt{\hat{R}}$

Onde

\hat{R} = \frac{\hat{V}}{W} \cdot \frac{d f + 3}{d f + 1},

$\hat{R} = \dfrac{\hat{V}}{W} \cdot \dfrac{df+3}{df+1}\,,$

d f = 2 \hat{V} / V a r (\hat{V})

$df = 2\hat{V}/Var(\hat{V})$

A documentação Stan na página 349 ignora o prazo com e também remove o multiplicativo prazo. Essa é a fórmula deles, $df$ $(M+1)/M$

O estimador de variância é Finalmente, a estatística de redução de escala potencial é definido por
${\hat{var}}^{+} (θ | y) = \frac{N - 1}{N} W + \frac{1}{N} B .$ $\widehat{\text{var}}^{+}(\theta \, | \, y) = \frac{N-1}{N} W + \frac{1}{N} B\,.$ $\hat{R} = \sqrt{\frac{{\hat{var}}^{+} (θ | y)}{W}} .$ $\hat{R} = \sqrt{\frac{\widehat{\text{var}}^{+}(\theta \, | \, y) }{W}}\,.$

Pelo que pude ver, eles não fornecem uma referência para essa mudança de fórmula e nem discutem isso. Geralmente não é muito grande e geralmente pode ser tão baixo quanto , portanto não deve ser ignorado, mesmo que o termo possa ser aproximado com 1. $M$ $2$ $(M+1)/M$ $df$

Então, de onde vem essa fórmula?

Edição: Encontrei uma resposta parcial para a pergunta "de onde vem essa fórmula? ", Em que o livro Bayesian Data Analysis de Gelman, Carlin, Stern e Rubin (Segunda edição) tem exatamente a mesma fórmula. No entanto, o livro não explica como / por que é justificável ignorar esses termos?

mcmc convergence gibbs metropolis-hastings stan Greenparker
fonte

Ainda não há um artigo publicado, e a fórmula provavelmente mudará nos próximos meses.

Ben Goodrich

@ BenGoodrich Obrigado pelo comentário. Você pode dizer mais alguma coisa sobre a motivação de usar essa fórmula? E por que exatamente a fórmula mudará?

Greenparker

A fórmula atual do R-hat dividido é a maneira mais comum de aplicá-la ao caso em que existe apenas uma cadeia. As próximas mudanças são principalmente para lidar com o fato de que a distribuição posterior marginal subjacente pode não ser normal ou ter uma média e / ou variação.

Ben Goodrich

@BenGoodrich Sim, entendo por que a STAN divide Rhat. Mas mesmo neste caso

, e assim a constante

, que não é ignorada.

M = 2

$M = 2$

(M + 1) / M = 3 / 2

$(M+1)/M = 3/2$

Greenparker

Respostas:

\hat{σ} = \frac{n - 1}{n} W + \frac{1}{n} B

$\hat{\sigma} = \frac{n-1}{n}W+ \frac{1}{n}B$

\hat{σ}

$\hat{\sigma}$

{\hat{σ}}_{+}

$\hat{\sigma}_+$

{\hat{v a r}}^{+}

$\widehat{\rm var}^+$

\hat{R} = \frac{m + 1}{m} \frac{{\hat{σ}}_{+}}{W} - \frac{n - 1}{m n},

$\hat{R} = \frac{m+1}{m}\frac{\hat{\sigma}_+}{W} - \frac{n-1}{mn},$ which can be rearranged as

\hat{R} = \frac{{\hat{σ}}_{+}}{W} + \frac{{\hat{σ}}_{+}}{W m} - \frac{n - 1}{m n} .

$\hat{R} = \frac{\hat{\sigma}_+}{W} + \frac{\hat{\sigma}_+}{Wm}- \frac{n-1}{mn}.$ We can see that the effect of second and third term are negligible for decision making when

n

$n$ is large. See also the discussion in the paragraph before Section 3.1 in Brooks & Gelman (1998).

Gelman & Rubin (1992) also had the term with df as df/(df-2). Brooks & Gelman (1998) have a section describing why this df corretion is incorrect and define (df+3)/(df+1). The paragraph before Section 3.1 in Brooks & Gelman (1998) explains why (d+3)/(d+1) can be dropped.

It seems your source for the equations was something post Brooks & Gelman (1998) as you had (d+3)/(d+1) there and Gelman & Rubin (1992) had df/df(-2). Otherwise Gelman & Rubin (1992) and Brooks & Gelman (1998) have equivalent equations (with slightly different notations and some terms are arranged differently). BDA2 (Gelman, et al., 2003) doesn't have anymore terms $\frac{\hat{\sigma}_+}{Wm}- \frac{n-1}{mn}$ . BDA3 (Gelman et al., 2003) and Stan introduced split chains version.

My interpretation of the papers and experiences using different versions of $\hat{R}$ is that the terms which have been eventually dropped can be ignored when $n$ is large, even when $m$ is not. I also vaguely remember discussing this with Andrew Gelman years ago, but if you want to be certain of the history, you should ask him.

Usually M is not too large, and can often be as low so as 2

I really do hope that this is not often the case. In cases where you want to use split- $\hat{R}$ convergence diagnostic, you should use at least 4 chains split and thus have M=8. You may use less chains, if you already know that in your specific cases the convergence and mixing is fast.

Additional reference:

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

Aki Vehtari
fonte

Yes it has the same

{\hat{σ}}^{2}

$\hat{\sigma}^2$ as you mention, but their

\hat{R}

$\hat{R}$ statistic is

({\hat{σ}}^{2} + B / m n) / W * d f_{t e r m}

$(\hat{\sigma}^2 + B/mn)/W * df_{term}$ (look at the equation on top of page 495 in the Stat Science official version), which introduces the

(m + 1) / m

$(m+1)/m$ term I was talking about. In addition, look at the code and description in the R package coda, which has had the GR diagnostic since 1999.

Greenparker

I'm confused. The article via the link you provided and the article from Stat Science web pages has only pages 457-472.I didn't check now, but years ago and last year when I checked coda, it didn't have the current recommended version.

Aki Vehtari

Note that I edited my answer. Gelman & Brooks (1998) has that (m+1)/m term more clearly, and it seems you missed the last term which mostly cancels the effect of (m+1)/m term for decision making. See that paragraph before section 3.1.

Aki Vehtari

Sorry about that, that was a typo. It's page 465, and Gelman and Rubin have the same exact definition as Brooks and Gelman (which you state above). Equation 1.1 in Brooks and Gelman is exactly what I wrote down as well (when you rearrange some terms).

Greenparker

"We can see that the effect of second and third term are negligible for decision making when n is large", so what you are saying is that the expression in BDA and hence STAN comes from essentially ignoring these terms for large n?

Greenparker