0

我想使用 pymc 使用 MH 链从自定义日志可能性中进行采样。但我似乎无法让它工作，也无法在网上找到一个像样的例子。

def getPymcVariable(data):

  def logp(value):
    ...
    return ljps # returns a float

  def random():
    return np.random.rand(numDims);

  dtype = type(random());
  initPt = [0.45, 0.24, 0.68];

  ret = pymc.Stochastic(logp = logp,
                        doc = 'SNLS RV',
                        name = 'SNLS',
                        parents = {},
                        random = random,
                        trace = True,
                        value = initPt,
                        dtype = dtype,
                        observed = False,
                        cache_depth = 2,
                        plot = True,
                        verbose = 0 );
  return ret


if __name__ == '__main__':

    data = np.loadtxt('../davisdata.txt');
    numExperiments = 30;
    numSamples = 10000;

    SNLS = getPymcVariable(data)
    model = pymc.Model([SNLS]);
    mcmcModel = pymc.MCMC(model);
    mcmcModel.use_step_method(pymc.Metropolis, SNLS, proposal_sd=1);
    mcmcModel.sample(numSamples, burn=0, thin=1);
    x = mcmcModel.trace('SNLS')[:]
    np.savetxt(fileName, x);

它是一个 3 维变量，具有统一的先验和在 logp() 中计算的对数似然。我想运行一个具有自适应提案分布的 MH 链。但是每次我运行采样器时它只返回一个均匀分布（事实上，它只是从上面的随机函数返回样本 - 当我将它修改为 0.5*np.random.rand(numDims) 它返回一个 Unif( (0, 1)^3) 分布。)

但是，我知道正在调用 logp 函数。

我还有几个问题： - 上面随机函数的目的是什么？我最初认为这是一个先验，但看起来不像。

4

1 回答 1

1

在 PyMC2 中，我发现使用内置发行版和@potential装饰器来完成此类任务要简单得多。这是一个最小的例子：

X = pm.Uniform('X', 0, 1, value=[0.45, 0.24, 0.68])

@pm.potential
def SNLS(X=X):
    logp = -X[0]**2 / X[1]
    logp += -X[1]**2 / X[2]  # or whatever...
    return logp

您可以选择自适应 Metropolis 步进方法，如下所示：

m = pm.MCMC([X, SNLS])
m.use_step_method(pm.AdaptiveMetropolis, X)

这是一个将这些放在一起并绘制结果的笔记本。

于 2014-06-03T13:26:46.510 回答