1

我emcee.EnsembleSampler用来抽样一些可能性。我在数百万个不同的数据集上执行此操作，因此性能很重要。

一旦满足自相关时间的某些收敛标准，如果能够让一个链停止采样，那就太好了。

我在文档中或通过我的搜索找不到任何这样做的方法，所以希望有人有一个食谱。

4

1 回答 1

3

虽然我在当前的 emcee 文档中也找不到任何内容，但最新版本似乎添加了一个教程，恰好显示了您正在尝试做的确切事情：link

如果链接断开或事情发生变化，这里是示例：

import emcee
import numpy as np
np.random.seed(42)

# The definition of the log probability function
# We'll also use the "blobs" feature to track the "log prior" for each step
def log_prob(theta):
    log_prior = -0.5 *  np.sum((theta-1.0)**2 / 100.0)
    log_prob = -0.5 * np.sum(theta**2) + log_prior
    return log_prob, log_prior

# Initialize the walkers
coords = np.random.randn(32, 5)
nwalkers, ndim = coords.shape

# Set up the backend
# Don't forget to clear it in case the file already exists
filename = "tutorial.h5"
backend = emcee.backends.HDFBackend(filename)
backend.reset(nwalkers, ndim)

# Initialize the sampler
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob, backend=backend)

max_n = 100000

# We'll track how the average autocorrelation time estimate changes
index = 0
autocorr = np.empty(max_n)

# This will be useful to testing convergence
old_tau = np.inf

# Now we'll sample for up to max_n steps
for sample in sampler.sample(coords, iterations=max_n, progress=True):
    # Only check convergence every 100 steps
    if sampler.iteration % 100:
        continue

    # Compute the autocorrelation time so far
    # Using tol=0 means that we'll always get an estimate even
    # if it isn't trustworthy
    tau = sampler.get_autocorr_time(tol=0)
    autocorr[index] = np.mean(tau)
    index += 1

    # Check convergence
    converged = np.all(tau * 100 < sampler.iteration)
    converged &= np.all(np.abs(old_tau - tau) / tau < 0.01)
    if converged:
        break
    old_tau = tau
于 2018-05-11T15:05:49.987 回答