Семплер Гиббса не сходится

Я пытался понять выборку Гиббса в течение некоторого времени. Недавно я видел видео, которое имело большой смысл.

https://www.youtube.com/watch?v=a_08GKWHFWo

Автор использовал выборку Гиббса для сходимости средних значений (theta_1 и theta_2) двумерного нормального распределения, используя следующий процесс:

init: инициализировать theta_2 случайным значением.

Петля:

  1. образец theta_1, обусловленный theta_2 как N~(p(theta_2), [1-p**2])
  2. образец theta_2, обусловленный theta_1 как N~(p(theta_1), [1-p**2])

(повторять до сходимости).

Я попробовал это самостоятельно и столкнулся с проблемой:

import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

rv = multivariate_normal(mean=[0.5, -0.2], cov=[[1, 0.9], [0.9, 1]])

rv.mean
>>> 
array([ 0.5, -0.2])

rv.cov
>>>
array([[1. , 0.9],
       [0.9, 1. ]])

import numpy as np
samples = []

curr_t2 = np.random.rand()
def gibbs(iterations=5000):
    theta_1 = np.random.normal(curr_t2, (1-0.9**2), None)
    theta_2 = np.random.normal(theta_1, (1-0.9**2), None)
    samples.append((theta_1,theta_2))
    for i in range(iterations-1):
        theta_1 = np.random.normal(theta_2, (1-0.9**2), None)
        theta_2 = np.random.normal(theta_1, (1-0.9**2), None)
        samples.append((theta_1,theta_2))
gibbs()

sum([a for a,b in samples])/len(samples)
>>>
4.745736136676516

sum([b for a,b in samples])/len(samples)
>>>
4.746816908769834

Теперь я вижу, где я накосячил. Я обнаружил, что тета_1 зависит от фактического значения тета_2, а не от его вероятности. Точно так же я обнаружил, что тета_2 зависит от фактического значения тета_1, а не от его вероятности.

Где я застрял, так это в том, как мне оценить вероятность того, что тета примет любое заданное наблюдаемое значение?

Я вижу два варианта: плотность вероятности (на основе положения на нормальной кривой) И p-значение (интегрирование от бесконечности (и/или отрицательной бесконечности) до наблюдаемого значения). Ни одно из этих решений не кажется «правильным».

Как мне поступить?


python random mcmc sampling normal-distribution
person jbuddy_13    schedule 12.06.2020    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Возможно, мое видео было недостаточно четким. Алгоритм не сходится «по средним значениям», а сходится к выборкам из распределения. Тем не менее, средние значения выборок из распределений будут сходиться к их соответствующим средним значениям.

Проблема с вашими условными средствами. В видео я выбираю предельные средние, равные нулю, для сокращения записи. Если у вас ненулевые предельные средние значения, условное ожидание двумерной нормали включает предельные средние значения, корреляция и стандартные отклонения (которые равны 1 в вашей двумерной нормальности). Обновленный код

import numpy as np
from scipy.stats import multivariate_normal

mu1 = 0.5
mu2 = -0.2
rv = multivariate_normal(mean=[mu1, mu2], cov=[[1, 0.9], [0.9, 1]])

samples = []

curr_t2 = np.random.rand()
def gibbs(iterations=5000):
    theta_1 = np.random.normal(mu1 + 0.9 * (curr_t2-mu2), (1-0.9**2), None)
    theta_2 = np.random.normal(mu2 + 0.9 * (theta_1-mu1), (1-0.9**2), None)
    samples.append((theta_1,theta_2))
    for i in range(iterations-1):
        theta_1 = np.random.normal(mu1 + 0.9 * (theta_2-mu2), (1-0.9**2), None)
        theta_2 = np.random.normal(mu2 + 0.9 * (theta_1-mu1), (1-0.9**2), None)
        samples.append((theta_1,theta_2))

gibbs()

sum([a for a,b in samples])/len(samples)
sum([b for a,b in samples])/len(samples)
person jaradniemi    schedule 16.06.2020
comment
Спасибо, это супер полезно! - person jbuddy_13; 17.06.2020