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

Questions
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