Creating complex custom priors and likelihood

Questions v5
1

Hello again and thanks in advance for the help I am trying to set up my priors with a custom function in pymc v5. I figured out that I need to use pm.DensityDist, but I’m probably missing something because my example does not work.

I have 20 variables that each has 3 degrees of freedom. In pymc language, if I understand it correctly, this would be 20 coordinates with 3 model dimensions. For each coordinate, and for each of its dimensions, I want to define a custom prior distribution. For example, the prior for dimension 1 and 2 is a combination of Normal distributions, for dimension 3 it is a uniform distribution.

My priors are constructed from data. Each coordinate may have different input data, so its prior will be different. For example, coordinate 1 may have a prior for its 1’st dimension that consists of 5 Normal distributions added together, same for coordinate 2, and for coordinate 3 it would be the addition of uniform distributions.

Now if all of the above made sense to you, I would really appreciate some guidance on what is the best way to construct these priors and sample the posteriors. I am still finding my way with probabilistic programming so this is the point where I am a bit lost. I have a code below that constructs the prior but gives it as a list, so it is not really compatible with pyMC sampling algorithms. Just to clarify: I will later call these priors in a custom likelihood function and evaluate the log likelihood.

Here is an example code of what I have:

import numpy as np
import pymc as pm

coordinates = 20
dimensions = 3


def prior_means_and_std(dimensions: int = 3,
                        coordinates: int = 20):
    '''Here i just create a test dataset to then create my prior distributions.
    The test dataset just gives the mean and std for the first two dimensions, and the lower & upper bounds
    for the uniform distribution. For each coordinate, I randomly set a number of datapoints to simulate what
    happens with my real dataset. The prior for each coordinate will be a combination of Normal / Uniform
    distributions that is equal to the number of datapoints'''
    output_table = {}
    for coordinate in range(coordinates):
        number_of_datapoints = np.random.randint(10)  # assign randomly the number of data for this coordinate
        means = np.random.uniform(size=(number_of_datapoints, dimensions - 1)) * 10 - 5
        stds = np.random.uniform(size=(number_of_datapoints, dimensions - 1)) * 3
        lower_lim = np.random.uniform(size=(number_of_datapoints, 1))
        upper_lim = lower_lim + np.random.uniform(size=(number_of_datapoints, 1))
        output_table['Coordinate ' + str(coordinate)] = [means, stds, lower_lim, upper_lim]
    return output_table


def custom_prior(table: dict):
    priors = {}
    for key, item in table.items():
        means = item[0]
        stds = item[1]
        lower_lim = item[2]
        upper_lim = item[3]
        p1 = 0
        p2 = 0
        p3 = 0
        for index in range(means.shape[0]):
            p1 += pm.Normal.dist(mu=means[index, 0], sigma=stds[index, 0])
            p2 += pm.Normal.dist(mu=means[index, 1], sigma=stds[index, 1])
            p3 += pm.Uniform.dist(lower=lower_lim[index], upper=upper_lim[index])
        priors[key] = [p1, p2, p3]
    return priors


# Create sample values for the means, std and lower/upper bounds
prior_table = prior_means_and_std()

# Now create the priors by adding together normal and uniform distributions from pyMC
# this will now return a dictionary with 20 entries, where for each entry there is a list with the three priors.
# But of course, this is not really compatible with what I want to do next..
priors = custom_prior(prior_table)

# I can call these functions and sample them
pm.draw(priors['Coordinate 19'][2])



#  But what I actually want to do is to use pm.sample and a custom likelihood function

def custom_likelihood(parameter_realization, observations):
    '''My likelihood function involves a more complicated forward model that I skip here. But
    it will require a realization of the coordinates to run this forward model and return the
    log likelihood'''

    def forward_model(parameters):
        return parameters

    synthetic_response = forward_model(parameter_realization)
    # For simplicity, I assume that the likelihood is the Normal distribution and just git the log-likelihood without
    # the scaling in front
    return ((synthetic_response - observations) ** 2).sum()


some_random_data = np.random.random(size=(coordinates, dimensions))
# Define the model
with pm.Model() as model:
    sampled_parameters = pm.DensityDist('params',
                                        logp=lambda value: custom_prior(), # here is one thing that is definitely done wrong
                                        shape=(coordinates,))
    observations = pm.Data('observed_data', some_random_data)
    likelihood = pm.DensityDist('likelihood',
                                logp=custom_likelihood(sampled_parameters, observations),
                                shape=(coordinates,))

# Sample the posterior
with model:
    trace = pm.sample(draws=1000, chains=4)
2

I would really appreciate some guidance on what is the best way to construct these priors and sample the posteriors.

Without knowing too much yet about the specifics of the problem, my recommendation would be to try to express the model in terms of a single prior distribution (remember - this prior can be multivariate) for each type of distribution you use, i.e. one for the normal, one for the uniform, one for the mixture of normals.

5 Normal distributions added together,

Here, I assume that you mean a mixture of normals - does that sound right?

Finally, this will all work most easily within PyMC if you can express the forward model within Aesara/PyTensor and/or the pymc math functions. These typically map 1-1 to functions available in Numpy. Do you think that is possible?

1 Like
3

Thanks for you quick reply. So the problem i’m solving is that I’m trying to create a geological model of faults, given observations in boreholes. I have several boreholes and in these I observe the depth and orientation of a fault. The orientation is defined by two angles, the dip and azimuth.

So each fault has a dip, azimuth and depth. The dip and azimuth are probably off by ± 10 degrees, and I define this prior uncertainty using a normal distribution. The depth can be off by a meter or so, probably better described by a uniform prior.

Now each of these observations can make a plane in 3D, and the plane is extended to all other boreholes. If it intersects at a place where another fault is observed, with a similar orientation and depth, then it is a probable match. Hence, the likelihood is defined by forward modeling all these observations and seeing where they cross other boreholes, and if they match existing observations. It’s almost like the game of connecting the dots.

To anwser your comments:

Without knowing too much yet about the specifics of the problem, my recommendation would be to try to express the model in terms of a single prior distribution (remember - this prior can be multivariate) for each type of distribution you use, i.e. one for the normal, one for the uniform, one for the mixture of normals.

I do not think I can have a single prior, since all these fault observations may be independent in the end of the day. That’s why I define 3 dimensions (dip, azimuth, depth) and several coordinates, corresponding to the faults.

Here, I assume that you mean a mixture of normals - does that sound right?

Sometimes a fault observation is defined by a bundle of faults that are nearby, probably informative of the overall structure. So several obervations may exist in a nearby region. I then bundle all these observations together by creating a mixture of normals (addind several normals together in pyMC) to characterize this uncertatiny. I guess you are pointing out that I can define directly a mixture of normals instead of adding together normals?

Finally, this will all work most easily within PyMC if you can express the forward model within Aesara/PyTensor and/or the pymc math functions. These typically map 1-1 to functions available in Numpy. Do you think that is possible?

This is possible. My problem is defined with very simple numpy computations and some for loops. Is this offering me anything more that pure speedup of the code? Perhaps gradient computations? This would be great, but i’d like to figure out how to properly set the priors first.

Thank you for the feedback!

4

To recap on this topic.

I managed to get what I need (as suggested by @ckrapu) with pm.Mixture. Here is a function that takes in a list of lists, where the inner lists can have any length. The output is a set of normal distributions equal to the length of the input list, and where if the inner list has more than 1 values, several normals are mixed. I use a fixed standard deviation but this is easily changed.

def mixture_of_normals(std: float,
                       list_of_means: list):
    logp = []
    for index, values in enumerate(list_of_means):
        inner_logp = []
        weights = np.ones(len(values)) / len(values)
        var_name = 'variable_' + str(index + 1)
        if len(values) > 1:
            for j, mean in enumerate(values):
                inner_logp.append(pm.Normal.dist(mu=mean, sigma=std))
            logp.append(pm.Mixture(var_name,
                                   w=weights,
                                   comp_dists=inner_logp))
        elif len(values) == 1:
            logp.append(pm.Normal(var_name, mu=values[0], sigma=std))
    return at.as_tensor_variable(logp)

Now I just need to translate my custom likelihood function to PyTensor so that I can run all this in parallel :slight_smile: Thanks for the input!