Как использовать MCMC с пользовательской логарифмической вероятностью и найти матрицу - PullRequest
Новая
       120

Как использовать MCMC с пользовательской логарифмической вероятностью и найти матрицу

1 голос
/ 26 марта 2019

Код в PyMC3, но это общая проблема.Я хочу найти, какая матрица (комбинация переменных) дает мне наибольшую вероятность.Брать среднее значение следа каждого элемента бессмысленно, потому что они зависят друг от друга.

Вот простой случай;код использует вектор, а не матрицу для простоты.Цель состоит в том, чтобы найти вектор длины 2, где каждое значение находится между 0 и 1, так что сумма равна 1. 

import numpy as np
import theano
import theano.tensor as tt
import pymc3 as mc

# define a theano Op for our likelihood function
class LogLike_Matrix(tt.Op):
    itypes = [tt.dvector] # expects a vector of parameter values when called
    otypes = [tt.dscalar] # outputs a single scalar value (the log likelihood)

    def __init__(self, loglike):
        self.likelihood = loglike        # the log-p function

    def perform(self, node, inputs, outputs):
        # the method that is used when calling the Op
        theta, = inputs  # this will contain my variables

        # call the log-likelihood function
        logl = self.likelihood(theta)

        outputs[0][0] = np.array(logl) # output the log-likelihood

def logLikelihood_Matrix(data):
    """
        We want sum(data) = 1
    """
    p = 1-np.abs(np.sum(data)-1)
    return np.log(p)

logl_matrix = LogLike_Matrix(logLikelihood_Matrix)

# use PyMC3 to sampler from log-likelihood
with mc.Model():
    """
        Data will be sampled randomly with uniform distribution
        because the log-p doesn't work on it
    """
    data_matrix = mc.Uniform('data_matrix', shape=(2), lower=0.0, upper=1.0)

    # convert m and c to a tensor vector
    theta = tt.as_tensor_variable(data_matrix)

    # use a DensityDist (use a lamdba function to "call" the Op)
    mc.DensityDist('likelihood_matrix', lambda v: logl_matrix(v), observed={'v': theta})

    trace_matrix = mc.sample(5000, tune=100, discard_tuned_samples=True)

1 Ответ

1 голос
/ 01 апреля 2019

Если вам нужны только значения параметра с наибольшим правдоподобием, вам нужна оценка Maximum A Posteriori (MAP), которую можно получить с помощью pymc3.find_MAP() (подробности метода см. starting.py).Если вы ожидаете мультимодальный апостериор, то вам, вероятно, придется запускать его повторно с разными инициализациями и выбирать тот, который получает наибольшее значение logp, но это все еще только увеличивает шансы на нахождение глобального оптимума, хотя и не может этого гарантировать.

Следует отметить, что при больших измерениях параметров оценка MAP обычно не является частью типичного набора, т. Е. Она не представляет типичные значения параметров, которые могли бы привести к наблюдаемым данным.Майкл Бетанкур обсуждает это в Концептуальное введение в гамильтониан Монте-Карло .Полностью байесовский подход заключается в использовании апостериорных прогностических распределений , которые эффективно усредняются по всем конфигурациям параметров с высокой вероятностью, а не с использованием оценки по одной точке для параметров.

...