Блог

Код находится в 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)
  

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

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