Вычисление отношения Хастингса, g (x | x ') / g (x' | x) для асимметричных распределений предложений в алгоритме MH?
Я понимаю алгоритм Метрополиса. Меня смущает алгоритм MH, в котором можно использовать асимметричные распределения предложений.

Я понимаю, что P (x) и P (x ') представляют собой вероятность / плотность вероятности x и x' согласно целевому распределению. Точно так же я понимаю, что g (x | x ') / g (x' | x) - это термин, используемый для исправления асимметричного распределения предложений. Меня не смущает его цель; Я не понимаю его исполнения.
В качестве игрушечной задачи я разработал семплер экспоненциального распределения. Есть два варианта: в одном используется симметричное распределение предложений, равномерное распределение. И один , что это не так : именно, Beta(a=3,b=2) - 0.5
. Я выбрал это распределение, потому что (A) оно асимметрично и в основном положительно (но иногда отрицательно из-за члена -0,5).
Понятия не имею, как найти g(x|x')/g(x'|x)
.
Код:
def target(x,lam):
return int(x>0) * lam * np.exp(-x * lam)
def exponential_MH(hops,lam=3):
states = []
burn_in = int(hops*0.2)
current = lam
for i in range(hops):
states.append(current)
# movement = current + random.uniform(-1,1) # does not require asymmetric correction
movement = current + np.random.beta(a=3,b=2)-0.5 # requires asymmetric correction
acceptance = target(x=movement,lam=lam)/target(x=current,lam=lam)
event = random.uniform(0,1)
if acceptance > event:
current = movement
return states[burn_in:]
lam = 1
exp_samples = exponential_MH(hops=10_000,lam=lam)
lines = np.linspace(0,5,10_000)
exp_curve = [lam*np.exp(-l*lam) for l in lines]
plt.hist(exp_samples,normed=1,bins=20)
plt.plot(lines,exp_curve)
График распределения ассиметричного предложения:

Симметричный график распределения предложений (с перевернутыми комментариями):

Чтобы ответить на этот вопрос, отредактируйте код, чтобы отразить действительное g(x|x')/g(x'|x)
для распределения предложения, из Beta(a=3,b=2) - 0.5
которого взяты возмущения.
Ответы
Взгляните на обновленный код и график ниже. Обратите внимание, что g (x | x ') / g (x' | x) по сути является мерой того, насколько вероятно, что эти возмущения будут видны при распределении предложений, которое было определено как Beta(a=3,b=2) -0.5
.
Сначала найдите разницу между текущими и предлагаемыми событиями. Во-вторых, установите -0,5; мы назовем эти несмещенные возмущения (где -0,5 - смещение.), в-третьих, найдите вероятность каждого возмущения (curr-> prop & prop-> curr). Наконец, верните соотношение как correction
.
Мы воспользуемся этим исправлением и умножим его на другие термины в acceptance
определении переменной. Вот и все!
def target(x,lam):
return int(x>0) * lam * np.exp(-x * lam)
def correct(prop,curr,a=3,b=2):
x0 = curr - prop + 0.5
x1 = prop - curr + 0.5
b0 = beta.pdf(x=x0, a=a, b=b)
b1 = beta.pdf(x=x1, a=a, b=b)
return b0/b1
def exponential_MH(hops,lam=3):
states = []
burn_in = int(hops*0.2)
current = 1
for i in range(hops):
states.append(current)
movement = current + np.random.beta(a=3,b=2)-0.5 # requires assymetric correction
correction = correct(curr=current,prop=movement)
acceptance = target(x=movement,lam=lam)/target(x=current,lam=lam)*correction
event = random.uniform(0,1)
if acceptance > event:
current = movement
return states[burn_in:]
