Styjun

The code is in PyMC3, but this is a general problem. I want to find which matrix (combination of variables) gives me the highest probability. Taking the mean of the trace of each element is meaningless because they depend on each other.

Here is a simple case; the code uses a vector rather than a matrix for simplicity. The goal is to find a vector of length 2, where the each value is between 0 and 1, so that the sum is 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)

If you only want the highest likelihood parameter values, then you want the Maximum A Posteriori (MAP) estimate, which can be obtained using pymc3.find_MAP() (see starting.py for method details). If you expect a multimodal posterior, then you will likely need to run this repeatedly with different initializations and select the one that obtains the largest logp value, but that still only increases the chances of finding the global optimum, though cannot guarantee it.

It should be noted that at high parameter dimensions, the MAP estimate is usually not part of the typical set, i.e., it is not representative of typical parameter values that would lead to the observed data. Michael Betancourt discusses this in A Conceptual Introduction to Hamiltonian Monte Carlo. The fully Bayesian approach is to use posterior predictive distributions, which effectively averages over all the high-likelihood parameter configurations rather than using a single point estimate for parameters.