Confiabilidade do modo de uma amostra do MCMC

12

Em seu livro Doing Bayesian Data Analysis, John Kruschke afirma que, ao usar JAGS de R

... a estimativa do modo de uma amostra do MCMC pode ser bastante instável porque a estimativa é baseada em um algoritmo de suavização que pode ser sensível a solavancos e ondulações aleatórios na amostra do MCMC. ( Fazendo a análise de dados bayesiana , página 205, seção 8.2.5.1)

Embora eu tenha uma noção do algoritmo Metropolis e formas exatas como a amostragem de Gibbs, não estou familiarizado com o algoritmo de suavização que também é mencionado e por que isso significaria que a estimativa do modo da amostra do MCMC é instável. Alguém é capaz de fornecer uma visão intuitiva do que o algoritmo de suavização está fazendo e por que torna instável a estimativa do modo?

Morgan Ball
fonte
2
John Kruschke falando sobre algoritmo para estimativa de modo que é baseado na estimativa de densidade do kernel.
Andrey Kolyadin
2
Este link pode ser útil.
Andrey Kolyadin
A menos que eu esteja errado em ser novo nesta área de estatística, o JAGS gera um conjunto de amostras da distribuição posterior, em vez de uma função de densidade de probabilidade, portanto, não tenho certeza se a estimativa da densidade kerna está presente. Obrigado pelo link embora.
Morgan Ball
Acho que isso provavelmente tem mais a ver com a forma como você obtém um modo de uma grande amostra de uma variável contínua, onde pode não haver mais de um valor específico e, portanto, é necessário agrupar (ou suavizar) a amostra.
Morgan Ball
1
Você pode obter o modo como valor com densidade máxima na estimativa de densidade do kernel. (Pelo menos é o que eu faço, e se não me engano J. Kruschke usar mesma abordagem em seus exemplos)
Andrey Kolyadin

Respostas:

8

Eu não tenho o livro em mãos, então não tenho certeza do método de suavização que Kruschke usa, mas, para intuição, considere esse gráfico de 100 amostras de um normal padrão, juntamente com as estimativas de densidade de kernel gaussiana usando várias larguras de banda de 0,1 a 1,0. (Resumidamente, os KDEs gaussianos são uma espécie de histograma suavizado: eles estimam a densidade adicionando um gaussiano para cada ponto de dados, com média no valor observado.)

Você pode ver que, mesmo quando a suavização cria uma distribuição unimodal, o modo geralmente fica abaixo do valor conhecido de 0.

insira a descrição da imagem aqui

Além disso, aqui está um gráfico do modo estimado (eixo y) por largura de banda do kernel usado para estimar a densidade, usando a mesma amostra. Esperemos que isso dê alguma intuição sobre como a estimativa varia com os parâmetros de suavização.

insira a descrição da imagem aqui

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Feb  1 09:35:51 2017

@author: seaneaster
"""

import numpy as np
from matplotlib import pylab as plt
from sklearn.neighbors import KernelDensity

REAL_MODE = 0
np.random.seed(123)

def estimate_mode(X, bandwidth = 0.75):
    kde = KernelDensity(kernel = 'gaussian', bandwidth = bandwidth).fit(X)
    u = np.linspace(-3,3,num=1000)[:, np.newaxis]
    log_density = kde.score_samples(u)
    return u[np.argmax(log_density)]

X = np.random.normal(REAL_MODE, size = 100)[:, np.newaxis] # keeping to standard normal

bandwidths = np.linspace(0.1, 1., num = 8)

plt.figure(0)
plt.hist(X, bins = 100, normed = True, alpha = 0.25)

for bandwidth in bandwidths:
    kde = KernelDensity(kernel = 'gaussian', bandwidth = bandwidth).fit(X)
    u = np.linspace(-3,3,num=1000)[:, np.newaxis]
    log_density = kde.score_samples(u)
    plt.plot(u, np.exp(log_density))

bandwidths = np.linspace(0.1, 3., num = 100)
modes = [estimate_mode(X, bandwidth) for bandwidth in bandwidths]
plt.figure(1)
plt.plot(bandwidths, np.array(modes))
Sean Easter
fonte
5

Sean Easter deu uma boa resposta; aqui está como é realmente feito pelos scripts R que acompanham o livro de Kruschke. A plotPost()função é definida no script R nomeado DBDA2E-utilities.R. Ele exibe o modo estimado. Dentro da definição da função, existem estas duas linhas:

mcmcDensity = density(paramSampleVec)
mo = mcmcDensity$x[which.max(mcmcDensity$y)]

A density()função vem com o pacote de estatísticas básicas de R e implementa um filtro de densidade de kernel do tipo descrito por Sean Easter. Possui argumentos opcionais para a largura de banda do kernel de suavização e para o tipo de kernel a ser usado. O padrão é um kernel gaussiano e possui alguma mágica interna para encontrar uma boa largura de banda. A density()função retorna um objeto com um componente nomeado yque possui as densidades suavizadas em vários valores x. A segunda linha de código, acima, apenas encontra o xvalor onde yé máximo.

John K. Kruschke
fonte