Issue with setting up pm.Potential("likelihood", logl(param))

Questions
1

I am currently setting up a Hierarchical Bayesian Model for Calibration of a Material Model.

When I setup my model I define a customised likelihood.

When defining the line pm.Potential(“likelihood”, logl(param))

pymc3 appears to test my customised likelihood, yet it does it with parameters a long way from the prior I have defined. As a result my model crashes / spits out rubbish.

Any ideas?

2

Hi, can you perhaps share more details about your model and logl function?

What do you get when you run the following line?

model.check_test_point()

Where model is the name you gave to your pm.Model() context.

3

Hi @ricardoV94, Thanks for your reply.

So first a bit of detail about the problem. I am using a customised likelihood which takes the form.

class LogLike(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, data, stressVals, theta, bounds_lower, bounds_upper, J, N_MC):
        # add inputs as class attributes
        self.likelihood = loglike
        self.data = data
        self.stressVals = stressVals
        self.theta = theta
        self.bounds_lower = bounds_lower
        self.bounds_upper = bounds_upper
        self.J = J
        self.N_MC = N_MC

    def perform(self, node, inputs, outputs):
        # the method that is used when calling the Op
        phi, = inputs  # this will contain my variables
        # call the log-likelihood function
        logl = self.likelihood(phi, self.stressVals, self.data, self.theta, self.bounds_lower, self.bounds_upper, J, N_MC)
        outputs[0][0] = np.array(logl) # output the log-likelihood

If I check the values of phi when this function is called they are all 1’s except a single value of 0.25.

There are 18 parameters for this hierarchical Bayesian model, means and standard deviations of 9 material parameters for the model. Ranging from values of 210.0 to 0.1

The model is defined as follows

with basic_model:

    BoundedNormal1 = pm.Bound(pm.Normal, lower=0.0)
    BoundedNormal2 = pm.Bound(pm.Normal, lower=-1.0, upper=0.5)

  

    mu_E1 = BoundedNormal1('mu_E1', mu=156.0, sigma=3.0)
    mu_E2 = BoundedNormal1('mu_E2', mu=153.0, sigma=3.0)
    mu_G12 = BoundedNormal1('mu_G12', mu=216., sigma=3.0)
    mu_nu12 = BoundedNormal2('mu_nu12', mu=0.29, sigma=0.01)

    mu_sigy0 = BoundedNormal1('mu_sigy0', mu=0.36, sigma=0.03)
    mu_sigy90 = BoundedNormal1('mu_sigy90', mu=0.33, sigma=0.03)
    mu_sigy45 = BoundedNormal1('mu_sigy45', mu=0.41, sigma=0.03)

    mu_n = BoundedNormal1('mu_n', mu = 13.0, sigma=1.0)
    
    mu_sigmaf = BoundedNormal1('mu_sigmaf', mu = 1.0, sigma=0.1)
    
    sig_E1 = BoundedNormal1('sig_E1', mu=10.0, sigma=1.0)
    sig_E2 = BoundedNormal1('sig_E2', mu=10.0, sigma=1.0)
    sig_G12 = BoundedNormal1('sig_G12', mu=10.0, sigma=1.0)
    sig_nu12 = BoundedNormal2('sig_nu12', mu=0.02, sigma=0.002)

    sig_sigy0 = BoundedNormal1('sig_sigy0', mu=0.03, sigma=0.005)
    sig_sigy90 = BoundedNormal1('sig_sigy90', mu=0.03, sigma=0.005)
    sig_sigy45 = BoundedNormal1('sig_sigy45', mu=0.03, sigma=0.005)

    sig_n = BoundedNormal1('sig_n', mu = 1.0, sigma=0.1)

    sig_sigmaf = BoundedNormal1('sig_sigmaf', mu = 0.1, sigma=0.01)


    # convert stochastic parameters to a tensor vector
    param = tt.as_tensor_variable([mu_E1, mu_E2, mu_G12, mu_nu12, mu_sigy0, mu_sigy90, mu_sigy45, mu_n, mu_sigmaf, sig_E1, sig_E2, sig_G12, sig_nu12, sig_sigy0, sig_sigy90, sig_sigy45, sig_n, sig_sigmaf])

    # use a DensityDist (use a lamdba function to "call" the Op)
    pm.Potential("likelihood", logl(param))

The final line in the definition of pm.Potential("likelihood", logl(param)) causes the issue; apparently when it tests out my likelihood. With values of all 1 my loglikelihood give silly values.

As for your test here is the output from calling

model.check_test_point()

mu_E1_lowerbound__         -1336.74
mu_E2_lowerbound__         -1285.57
mu_G12_lowerbound__        -2570.07
mu_nu12_interval__         -1455.29
mu_sigy0_lowerbound__       -224.97
mu_sigy90_lowerbound__      -246.80
mu_sigy45_lowerbound__      -190.80
mu_n_lowerbound__            -72.92
mu_sigmaf_lowerbound__         1.38
sig_E1_lowerbound__          -41.42
sig_E2_lowerbound__          -41.42
sig_G12_lowerbound__         -41.42
sig_nu12_interval__        -9108.19
sig_sigy0_lowerbound__    -18813.62
sig_sigy90_lowerbound__   -18813.62
sig_sigy45_lowerbound__   -18813.62
sig_n_lowerbound__             1.38
sig_sigmaf_lowerbound__    -4046.31
Name: Log-probability of test_point, dtype: float64

Which support this!

Thanks for any suggestions! Probably something silly!

Tim

4

At a first glance it seems like you have a positive log likelihood for 2 of your parameters. That looks wrong, but maybe it’s just unnormalized PDFs?

5

so if I look at model.test_point these are all zero

Issue is for this value my model will be dividing by zero in some cases. But for my prior, there is zero probability of these parameters being zero.

Even if I change the lower bound on some parameters, test_point still gives zero values!

Am I being stupid . . :man_shrugging:

6

Edit: I misunderstood, yes the test_point gives you the initial point. You can change it by setting testval=x when you initialize your distribution. Try to set those that cannot be zero some other valid value.

7

Thanks for your quick reply! Ok, Let me have a dig, I might start with a simpler model to see where I am doing something silly :slight_smile: