Metropolis Hastings с вероятностью Custom Log в Pymc

Я хочу использовать 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);

Это трехмерная переменная, имеет однородный априор и логарифмическую вероятность, вычисленную в logp(). Я хочу запустить сеть MH с адаптивным распределением предложений. Но каждый раз, когда я запускаю сэмплер, он просто возвращает равномерное распределение (на самом деле, он просто возвращает выборки из приведенной выше случайной функции — когда я изменил ее на 0,5*np.random.rand(numDims), она вернула Unif((0, 1)^3) распределение. )

Однако я знаю, что вызывается функция logp.

У меня есть еще несколько вопросов: - Какова цель случайной функции выше? Я сначала подумал, что это приора, но не похоже.


mcmc pymc
person tetradeca7tope    schedule 02.06.2014    source источник

Ответы (1)


arrow_upward
1
arrow_downward

В 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

Вы можете выбрать метод адаптивного шага мегаполиса следующим образом:

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

Вот записная книжка, в которой все это собрано вместе и представлены результаты.

person Abraham D Flaxman    schedule 03.06.2014