Почему закаленные mcmc подходят плохо?

Question

Я пытаюсь разместить простой прямой тип y=mx+c для некоторых синтетических данных, используя parallel-tempered mcmc.Моя цель - просто понять, как его использовать, чтобы позже я смог применить к более сложным моделям.Пример, который я пытаюсь использовать, является копией того, что уже было сделано в простом коде ведущего: http://dfm.io/emcee/current/user/line/, но вместо использования mcmc я хочу использовать mcmc с параллельным отпуском: http://dfm.io/emcee/current/user/pt/

Вот рабочий код:

import numpy as np
from emcee import PTSampler
import emcee

# Choose the "true" parameters.
m_true = -0.9594
b_true = 4.294
f_true = 0.534

# Generate some synthetic data from the model.
N = 50
x = np.sort(10*np.random.rand(N))
yerr = 0.1+0.5*np.random.rand(N)
y = m_true*x+b_true
y += np.abs(f_true*y) * np.random.randn(N)
y += yerr * np.random.randn(N)


def lnlike(theta, x, y, yerr):
    m, b, lnf = theta
    model = m * x + b
    inv_sigma2 = 1.0/(yerr**2 + model**2*np.exp(2*lnf))
    return -0.5*(np.sum((y-model)**2*inv_sigma2 - np.log(inv_sigma2)))
def lnprior(theta):
    m, b, lnf = theta
    if -5.0 < m < 0.5 and 0.0 < b < 10.0 and -10.0 < lnf < 1.0:
        return 0.0
    return -np.inf
def lnprob(theta, x, y, yerr):
    lp = lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(theta, x, y, yerr)

import scipy.optimize as op
nll = lambda *args: -lnlike(*args)
result = op.minimize(nll, [m_true, b_true, np.log(f_true)], args=(x, y, yerr))
m_ml, b_ml, lnf_ml = result["x"]
init = [0.5, m_ml, b_ml, lnf_ml]

ntemps = 10
nwalkers = 100
ndim = 3
from multiprocessing import Pool
pos = np.random.uniform(low=-1, high=1, size=(ntemps, nwalkers, ndim))
for i in range(ntemps):
    #initialize parameters near scipy optima
    pos[i:,] = np.array([result["x"] + 1e-4*np.random.randn(ndim) for i in range(nwalkers)])
pool = Pool(processes=4)
sampler=PTSampler(ntemps,nwalkers, ndim, lnlike, lnprior, loglargs=(x, y, yerr), pool=pool)# args=(x, y, yerr))
#burn-in 
sampler.run_mcmc(pos, 1000)
sampler.reset()
sampler.run_mcmc(pos, 10000, thin=10)
samples = sampler.chain.reshape((-1, ndim))
print('Number of posterior samples is {}'.format(samples.shape[0]))
#print best fit value together with errors
print(map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
                             zip(*np.percentile(samples, [16, 50, 84],
                                                axis=0))))

import corner
fig = corner.corner(samples, labels=["$m$", "$b$", "$\ln\,f$"],
                      truths=[m_true, b_true, np.log(f_true)])
fig.savefig("triangle.png")

Единственная проблема при запуске этого кода состоит в том, что я получаю оптимальные значения параметров, которые далеко от истинных значений.И увеличение количества ходунков или образцов не помогает ни в каком смысле.Может кто-нибудь, пожалуйста, посоветуйте, почему tempered-mcmc здесь не работает?

Обновление:

Я обнаружил полезный пакет под названием ptemcee (https://pypi.org/project/ptemcee/#description),хотя документация по этому пакету не существует. Похоже, этот пакет может быть полезен, любая помощь по реализации того же линейного фитинга с этим пакетом также будет высоко оценена.

Answer 1

Я изменил некоторые строки

import time
import numpy as np
from emcee import PTSampler
import corner
import matplotlib.pyplot as plt
import scipy.optimize as op

t1 = time.time()

np.random.seed(6) # To reproduce results
# Choose the "true" parameters.
m_true = -0.9594
b_true = 4.294
f_true = 0.534

# Generate some synthetic data from the model.
N = 50
x = np.sort(10 * np.random.rand(N))
yerr = 0.1 + 0.5 * np.random.rand(N)
y_1 = m_true * x + b_true
y = np.abs(f_true * y_1) * np.random.randn(N) + y_1
y += yerr * np.random.randn(N)

plt.plot(x, y, 'o')


# With emcee

def lnlike(theta, x, y, yerr):
    m, b, lnf = theta
    model = m * x + b
    inv_sigma2 = 1.0/(yerr**2 + model**2*np.exp(2*lnf))
    return -0.5*(np.sum((y-model)**2*inv_sigma2 - np.log(inv_sigma2)))
def lnprior(theta):
    m, b, lnf = theta
    if -5.0 < m < 0.5 and 0.0 < b < 10.0 and -10.0 < lnf < 1.0:
        return 0.0
    return -np.inf
def lnprob(theta, x, y, yerr):
    lp = lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(theta, x, y, yerr)


nll = lambda *args: -lnlike(*args)
result = op.minimize(nll, [m_true, b_true, np.log(f_true)], args=(x, y, yerr))
m_ml, b_ml, lnf_ml = result["x"]

init = [0.5, m_ml, b_ml, lnf_ml]

ntemps = 10
nwalkers = 100
ndim = 3

pos = np.random.uniform(low=-1, high=1, size=(ntemps, nwalkers, ndim))
for i in range(ntemps):
    pos[i:, :] = np.array([result["x"] + 1e-4*np.random.randn(ndim) for i in range(nwalkers)])

sampler = PTSampler(ntemps, nwalkers, ndim, lnlike, lnprior, loglargs=(x, y, yerr), threads=4) # args=(x, y, yerr))

#burn-in 
print(pos.shape)
sampler.run_mcmc(pos, 100)
sampler.reset()
sampler.run_mcmc(pos, 5000, thin=10)
samples = sampler.chain.reshape((-1, ndim))

print('Number of posterior samples is {}'.format(samples.shape[0]))

#print best fit value together with errors

p1, p2, p3 = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
                             zip(*np.percentile(samples, [16, 50, 84],
                                                axis=0)))

print(p1, '\n', p2, '\n', p3)

fig = corner.corner(samples, labels=["$m$", "$b$", "$\ln\,f$"],
                      truths=[m_true, b_true, np.log(f_true)])

t2 = time.time()

print('It took {:.3f} s'.format(t2 - t1))

plt.show()

Число, которое я получаю с corner:

Важная строка

sampler = PTSampler(ntemps, nwalkers, ndim, lnlike, lnprior, loglargs=(x, y, yerr), threads=4)

Я использовал threads=4 вместо Pool.

Посмотрите внимательно на эту строку print(p1, '\n', p2, '\n', p3), она печатает значения m_true, b_true и f_true, которые вы получите:

(-1.277782877669762, 0.5745273177144817, 2.0813620981463297) 
(4.800481378230051, 3.1747356851201163, 2.245189235990341) 
(-0.9391847529845194, 1.1196053087321716, 3.6017609114364273)

Для f вам нужно np.exp(-0.93918), что составляет 0.3909, что близко к 0.534. Полученные значения близки (-1.277 по сравнению с -0.9594 и 4.8 по сравнению с 4.294), хотя ошибки неплохие (кроме f). Я имею в виду, вы ожидаете получить точные цифры? При использовании этого метода на моем компьютере требуется 111 с, это нормально?

Давайте попробуем что-то другое. Давайте проясним: проблема не проста, когда добавлено f_true. Я буду использовать pymc3 (вам не нужно знать, как использовать pymc3, я хочу проверить результаты, найденные emcee).

import time
import numpy as np
import corner
import matplotlib.pyplot as plt
import pymc3 as pm

t1 = time.time()

np.random.seed(6)
# Choose the "true" parameters.
m_true = -0.9594
b_true = 4.294
f_true = 0.534

# Generate some synthetic data from the model.
N = 50
x = np.sort(10 * np.random.rand(N))
yerr = 0.1 + 0.5 * np.random.rand(N)
y_1 = m_true * x + b_true
y = np.abs(f_true * y_1) * np.random.randn(N) + y_1
y += yerr * np.random.randn(N)

plt.plot(x, y, 'o')


with pm.Model() as model: # model specifications in PyMC3 are wrapped in a with-statement
    # Define priors
    f = pm.HalfCauchy('f', beta=5)
    m = pm.Normal('m', 0, sd=20)
    b = pm.Normal('b', 0, sd=20)

    mu2 = b + m * x
    sigma2 = yerr**2 + f**2 * (y_1)**2

    post = pm.Normal('y', mu=mu2, sd=pm.math.sqrt(sigma2), observed=y)

with model:
    trace = pm.sample(2000, tune=2000)

print(pm.summary(trace))
pm.traceplot(trace)

all_values = np.stack([trace.get_values('b'), trace.get_values('m'), trace.get_values('f')], axis=1)

fig2 = corner.corner(all_values, labels=["$b$", "$m$", "$f$"],
                      truths=[b_true, m_true, f_true])

t2 = time.time()

print('It took {:.3f} s'.format(t2 - t1))

plt.show()

Резюме

       mean        sd  mc_error   hpd_2.5  hpd_97.5        n_eff      Rhat
m -0.995545  0.067818  0.001174 -1.123187 -0.857653  2685.610018  1.000121
b  4.398158  0.332526  0.005585  3.767336  5.057909  2746.736563  1.000201
f  0.425442  0.063884  0.000904  0.311037  0.554446  4195.591204  1.000309

Важной частью является столбец mean, вы видите, что значения, найденные pymc3, близки к истинным значениям. Столбцы hpd_2.5 и hpd_97.5 являются ошибками для f, b и m. И это заняло 14 с.

Число, которое я получаю с corner, равно

Вы скажете, что результаты emcee не совсем хороши, но если вы действительно хотите большей точности, вам нужно изменить эту функцию:

def lnprior(theta):
    m, b, lnf = theta
    if -5.0 < m < 0.5 and 0.0 < b < 10.0 and -10.0 < lnf < 1.0:
        return 0.0
    return -np.inf

Знаменитый до . В этом случае он плоский, а так как много приоры ...