Complex model on PYMC3 - Nonlinear Hierarchical model?

Questions
1

Hi!

I’d been studying the concepts of probabilistic programming using Pymc3 in the past couple of week and I have a question on how to implement a particular model, I will describe it next:

I have a set of N data elements, and a matrix NxN defining the interactions between each other of those elements (the matrix is symmetric). I also have a property E for each element, I want to model this property as it follows:

eqq

What the model says is that, the energy of a element M is a weighted sum over all the other elements in the data set. The function d(M, M_i) is known ( is the NxN matrix I described before)

I want to model the weights alpha_i and sigma using Bayesian inference.

My try so far goes like this:

    indx = np.arange(0, N   )
    # The index for the N elements 
    def estimator( alpha ,x , sigma , ):
        # My idea is that this estimator returns the estimated value for the energy for a given set alpha, x , and sigma
        return T.sum(  alpha*T.exp( -1/(2.0*sigma**2 )*x**2  )  )


    with pm.Model() as hi_model:

        # Positive prior on sigma, N positive priors for the alphas 
        sigma = pm.HalfCauchy( 'sigma' , beta = 1 , shape = 1 )
        alpha1 = pm.HalfCauchy( 'alphas' , beta = 1 , shape = N  )
        # An error parameter 
        sigma_e = pm.HalfCauchy('sigma_e' ,  beta = 1  )
        tau_e = sigma_e**-2
        # The energy estimator build using the function above. 
        # D2s is an NxN matrix, whose column i is a vector of length N
        #representing the interaction of the i element with all the other N elements 
        energy_est =  estimator( alpha1 , D2s[indx] , sigma  )  
        ee = pm.Deterministic( 'muu' , energy_est )
        y_like = pm.Normal('y_like', mu= ee , tau=tau_e, observed= energy_smp[:N]  )

Pymc3 allow me to build and sample from this model without errors, but I’m not being able to use the sampled alphas an sigma to “predict” the N energies I’m giving as observed, I’m off by many orders of magnitude. So I’m concerned if the model is being understood as I want or if I’m missing something and thats not the way to build this kind of model.

Thanks in advance for the help :slight_smile:

2

I suggest the first thing to try is setting tau_e below to a large value (i.e., setting the sigma_e to a deterministic small constant)

Otherwise, it is possible that most of the uncertainty is explained by tau_e, and the information is not propagated to alpha and sigma.

Another more direct approach is try to expressed E_est as a likelihood, and wrap it with pm.DensityDist

3

Hi , thanks for your answer. I will try the suggested idea with sigma and see what happens.

Can you refer me to some example on how to use pm.DensityDist , and maybe a little guide on how to use it in this case.

Thanks.

4

You can do a search in this discourse, there are a few related discussions: https://discourse.pymc.io/search?q=DensityDist
Also there are a few examples in the docs.pymc.io