时间戳

首页 > 解决方案 > 如何使用具有自定义对数概率的 MCMC 并求解矩阵

pymc3 - 如何使用具有自定义对数概率的 MCMC 并求解矩阵

问题描述

代码在 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)

标签： pymc3pymcmcmc

解决方案

如果您只需要最高似然参数值，那么您需要最大 A 后验 (MAP) 估计，该估计可以使用pymc3.find_MAP()（参见starting.py方法详细信息）获得。如果您期望多峰后验，那么您可能需要使用不同的初始化重复运行它并选择获得最大值的logp那个，但这仍然只会增加找到全局最优值的机会，尽管不能保证它。

应该注意的是，在高参数维度上，MAP 估计通常不是典型集合的一部分，即，它不代表会导致观察到的数据的典型参数值。Michael Betancourt 在A Conceptual Introduction to Hamiltonian Monte Carlo中讨论了这一点。完全贝叶斯方法是使用后验预测分布，它有效地对所有高似然参数配置进行平均，而不是对参数使用单点估计。

推荐阅读