<p>If you only want the highest likelihood parameter values, then you want the Maximum A Posteriori (MAP) estimate, which can be obtained using <code>pymc3.find_MAP()</code> (see <a href="https://github.com/pymc-devs/pymc3/blob/master/pymc3/tuning/starting.py" rel="nofollow noreferrer"><code>starting.py</code></a> 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 <code>logp</code> value, but that still only increases the chances of finding the global optimum, though cannot guarantee it.</p>

<p>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 href="https://arxiv.org/abs/1701.02434" rel="nofollow noreferrer">A Conceptual Introduction to Hamiltonian Monte Carlo</a>. The fully Bayesian approach is to use <a href="https://docs.pymc.io/api/inference.html#pymc3.sampling.sample_posterior_predictive" rel="nofollow noreferrer">posterior predictive distributions</a>, which effectively averages over all the high-likelihood parameter configurations rather than using a single point estimate for parameters.</p>
How to Use MCMC with a Custom Log-Probability and Solve for a Matrix
<p>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.</p>

<p>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.</p>

<pre><code>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)
</code></pre>