Diagnóstico MCMC Geweke

14

Estou executando um amostrador Metropolis (C ++) e quero usar os exemplos anteriores para estimar a taxa de convergência.

Um diagnóstico fácil de implementar que encontrei é o diagnóstico Geweke , que calcula a diferença entre as duas médias da amostra divididas pelo erro padrão estimado. O erro padrão é estimado a partir da densidade espectral em zero.

Zn=θ¯Aθ¯B1nASθA^(0)+1nBSθB^(0),

onde , B são duas janelas dentro da cadeia de Markov. Eu fiz algumas pesquisas sobre o que são ^ S A θ ( 0 ) e ^ S B θ ( 0 ), mas entrei em uma confusão de literatura sobre densidade espectral de energia e densidade espectral de potência, mas não sou especialista nesses tópicos; Eu só preciso de uma resposta rápida: essas quantidades são iguais à variação da amostra? Se não, qual é a fórmula para calculá-los?UMABSθUMA^(0 0)SθB^(0 0)

Outra dúvida sobre esse diagnóstico de Geweke é como escolher ? A literatura acima disse que é algum θ ( X ) funcional e deveria implicar a existência de uma densidade espectral ^ S A θ ( 0 ) , mas por conveniência, acho que a maneira mais simples é usar a função de identidade (use as próprias amostras). Isso está correto?θθ(X)SθUMA^(0 0)

O pacote R coda tem uma descrição, mas também não especifica como calcular os valores S

Yang
fonte
você pode olhar nas entranhas da codafunção geweke.diagpara ver o que ele faz ...
Ben Bolker

Respostas:

4

Você pode procurar no código a geweke.diagfunção no codapacote para ver como a variação é calculada, através da chamada para a spectrum.ar0função.


Aqui está uma breve motivação do cálculo da densidade espectral de um processo AR ( ) em zero.p

A densidade espectral de um processo AR ( ) na frequência λpλ é dada pela expressão: ondeαjsão os parâmetros autorregressivos.

f(λ)=σ2(1-j=1pαjexp(-2πιjλ))2
αj

Essa expressão simplifica consideravelmente ao calcular a densidade espectral de um AR (p0 0

f(0 0)=σ2(1-j=1pαj)2

A computação ficaria assim (substituindo os estimadores usuais por parâmetros):

tsAR2 = arima.sim(list(ar = c(0.01, 0.03)), n = 1000)  # simulate an AR(2) process
ar2 = ar(tsAR2, aic = TRUE)  # estimate it with AIC complexity selection

# manual estimate of spectral density at zero
sdMan = ar2$var.pred/(1-sum(ar2$ar))^2

# coda computation of spectral density at zer0
sdCoda = coda::spectrum0.ar(tsAr2)$spec

# assert equality
all.equal(sdCoda, sdMan)
tchakravarty
fonte
0

Sxx(ω)Sxx(0 0)

xuhdev
fonte