Um método de amostragem de "importância Gibbs" funcionaria?

8

Suspeito que esta seja uma pergunta bastante incomum e exploratória, por isso, tenha paciência comigo.

Gostaria de saber se é possível aplicar a idéia de amostragem importante à amostragem de Gibbs. Aqui está o que quero dizer: na amostragem de Gibbs, alteramos o valor de uma variável (ou bloco de variáveis) de cada vez, amostrando a partir da probabilidade condicional, dadas as demais variáveis.

No entanto, pode não ser possível ou fácil amostrar a partir da probabilidade condicional exata. Então, em vez disso, coletamos amostras de uma distribuição de proposta e usamos, por exemplo, Metropolis-Hastings (MH).q

Por enquanto, tudo bem. Mas eis um caminho divergente: o que acontece se, em vez de usar o MH, usarmos a mesma idéia usada na amostragem de importância, ou seja, coletamos amostras de e mantemos um peso de importância da amostra atual?qp/q

Mais detalhadamente: suponha que temos as variáveis e uma distribuição fatorada modo que . Mantemos a probabilidade da proposta usada para amostrar o valor atual de cada variável . Em cada etapa, alteramos um subconjunto das variáveis ​​e atualizamos (apenas os fatores de e que são afetados). Tomamos as amostras e seu peso de importância para calcular qualquer estatística em que estamos interessados.x1,,xnϕ1,,ϕmpi=1mϕiqixip(x)/q(x)pq

Esse algoritmo estaria correto? Caso contrário, existem razões claras por que não? Intuitivamente, faz sentido para mim, pois parece estar fazendo a mesma coisa que a amostragem de importância, mas com amostras dependentes.

Eu implementei isso para um modelo de caminhada aleatória gaussiana e observei que os pesos se tornam cada vez menores (mas não monotonicamente); portanto, as amostras iniciais acabam tendo muita importância e dominam a estatística. Tenho certeza de que a implementação não é de buggy, porque a cada passo eu comparo o peso atualizado com um cálculo explícito de força bruta. Observe que os pesos não diminuem indefinidamente para zero, porque são onde e são produtos de um número finito de densidades, e cada amostra é obtida de uma distribuição Normal que raramente será zero.p/qpq

Então, estou tentando entender por que os pesos caem dessa maneira e se isso é uma conseqüência do fato de esse método não estar realmente correto.


Aqui está uma definição mais precisa do algoritmo, aplicada a uma caminhada aleatória gaussiana nas variáveis . O código segue abaixo.X1,,Xn

O modelo é simplesmente , com fixo em .XiN(Xi1,σ2),i=1,,nX00

O peso da amostra atual é , onde são as densidades gaussianas e são as distribuições das quais os valores atuais foram amostrados. Inicialmente, simplesmente amostramos os valores de maneira direta, então e o peso inicial é .ip(xi)iq(xi)pqq=p1

Em cada etapa, eu escolho para alterar. Eu um novo valor para de , portanto essa densidade se torna a nova distribuição de proposta usada para .j{1,,n}xjXjN(Xj1,σ2)Xj

Para atualizar o peso, divido pelo densidades e de valor antigo de acordo com a e , e multiplicar pelo densidades e de novo valor de acordo com a e . Isso atualiza o numerador do peso.p(xj|xj1)p(xj+1|xj)xjxj1xj+1p(xj|xj1)p(xj+1|xj)xjxj1xj+1p

Para atualizar o denominador , multiplico o peso pela proposta antiga (removendo-o do denominador) e divido-o por .qq(xj)q(xj)

(Como eu do normal centrado em , é sempre igual a então eles são cancelados e a implementação não usá-los).xjxj1q(xj)p(xj|xj1)

Como mencionei antes, no código eu comparo esse cálculo de peso incremental com o cálculo explícito real apenas para ter certeza.


Aqui está o código para referência.

println("Original sample: " + currentSample);
int flippedVariablesIndex = 1 + getRandom().nextInt(getVariables().size() - 1);
println("Flipping: " + flippedVariablesIndex);
double oldValue = getValue(currentSample, flippedVariablesIndex);
NormalDistribution normalFromBack = getNormalDistribution(getValue(currentSample, flippedVariablesIndex - 1));
double previousP = normalFromBack.density(oldValue);
double newValue = normalFromBack.sample();
currentSample.set(getVariable(flippedVariablesIndex), newValue);
double previousQ = fromVariableToQ.get(getVariable(flippedVariablesIndex));
fromVariableToQ.put(getVariable(flippedVariablesIndex), normalFromBack.density(newValue));
if (flippedVariablesIndex < length - 1) {
    NormalDistribution normal = getNormalDistribution(getValue(currentSample, flippedVariablesIndex + 1));
    double oldForwardPotential = normal.density(oldValue);
    double newForwardPotential = normal.density(newValue);
    // println("Removing old forward potential " + oldForwardPotential);
    currentSample.removePotential(new DoublePotential(oldForwardPotential));
    // println("Multiplying new forward potential " + newForwardPotential);
    currentSample.updatePotential(new DoublePotential(newForwardPotential));
}

// println("Removing old backward potential " + previousP);
currentSample.removePotential(new DoublePotential(previousP));
// println("Multiplying (removing from divisor) old q " + previousQ);
currentSample.updatePotential(new DoublePotential(previousQ));

println("Final sample: " + currentSample);
println();

// check by comparison to brute force calculation of weight:
double productOfPs = 1.0;
for (int i = 1; i != length; i++) {
    productOfPs *= getNormalDistribution(getValue(currentSample, i - 1)).density(getValue(currentSample, i));
}
double productOfQs = Util.fold(fromVariableToQ.values(), (p1, p2) -> p1*p2, 1.0);
double weight = productOfPs/productOfQs;
if (Math.abs(weight - currentSample.getPotential().doubleValue()) > 0.0000001) {
    println("Error in weight calculation");
    System.exit(0);
}
user118967
fonte
A amostragem de importância não fornece amostras da distribuição de destino (neste caso, os condicionais completos de ). Portanto, a dinâmica do kernel de Markov que gera convergência do MCMC não se sustenta. Sem olhar para o seu código, eu não posso ver porque os pesos vão 0.ϕi
Greenparker
Obrigado. Acho que vou ter que me aprofundar nos teoremas da convergência do MCMC. Incluí o código por precaução, é bastante simples. Obrigado.
user118967
1
Em vez de incluir o código bruto (ou além disso), você pode explicar como está implementando o algoritmo? Qual é a distribuição alvo, quais são as condicionais completas, o que é a distribuição proposta, como você está combinando os pesos, etc etc
Greenparker
Obrigado. Eu fiz isso, por favor, deixe-me saber se isso é confuso em algum lugar.
user118967
@ Xi'an: aqui, a amostragem de importância está sendo aplicada ao movimento de uma única variável. Em vez de aceitar a proposta ou não, como em Metropolis Hastings, sempre a aceitamos, mas mantemos uma medida de importância desse flip, dividindo a probabilidade p pela proposta q para a variável que está sendo invertida.
user118967

Respostas:

4

Essa é uma ideia interessante, mas vejo várias dificuldades com ela:

  1. ao contrário da amostragem de importância padrão, ou mesmo da amostragem de importância metropolitana, a proposta não está atuando no mesmo espaço que a distribuição de destino, mas em um espaço de menor dimensão, portanto a validação não é clara [e pode ser imposta a manter pesos nas iterações e, portanto, enfrentar a degeneração]
  2. as constantes de normalização ausentes nos condicionais completos são alteradas a cada iteração, mas não são contabilizadas [veja abaixo]
  3. os pesos não são limitados, pois ao longo das iterações, eventualmente haverá simulações com um peso muito grande, a menos que se mantenha o controle da última ocorrência de uma atualização para o mesmo índice , que pode colidir com a validação markoviana do amostrador Gibbs . A execução de um experimento modesto com e iterações mostra uma variedade de pesos de até .jn=2T=1037.656397e-073.699364e+04

Para entrar em mais detalhes, considere um alvo bidimensional , incluindo a constante de normalização adequada, e implemente a importância do amostrador de Gibbs com as propostas e . Pesos de importância corretos [no sentido de produzir a expectativa correta, isto é, um estimador imparcial, para uma função arbitrária de ] para simulações sucessivas são onde e são os marginais de . Ou equivalente p(,)qX(|y)qY(|x)( X , Y ) p ( x t , y(X,Y)

p(xt,yt1)qX(xt|yt1)mY(yt1)orp(xt1,yt)qY(yt|xt1)mX(xt1)
mX()mY()p(,)
pX(xt|yt1)qX(xt|yt1)orpY(yt|xt1)qY(yt|xt1)
Nos dois casos, isso requer as densidades marginais [intratáveis] de e abaixo do destino .XYp(,)

Vale a pena comparar o que acontece aqui com o algoritmo Metropolis de importância paralela . (Veja, por exemplo, Schuster und Klebanov, 2018. ) Se o destino for novamente e a proposta for , a importância ponderada está correto [para produzir uma estimativa imparcial] e não atualiza o peso anterior, mas começa do zero a cada iteração.p(,)q(,|x,y)

p(x,y)q(x,y|x,y)

(C.) Uma correção para a importância original da proposta de Gibbs é propor um novo valor para todo o vetor, por exemplo, , da proposta de Gibbs , porque então o peso da importância está correto [faltando uma possível normalização constante que agora é verdadeiramente constante e não carrega das iterações anteriores de Gibbs] .(x,y)qX(xt|yt1)qY(yt|xt)

p(xt,yt)qX(xt|yt1)qY(yt|xt)

Uma observação final: para o destino de caminhada aleatória considerado no código, a simulação direta é possível em cascata: simule , depois com , etc.X1X2X1

Xi'an
fonte