Defining non-standard likelihood

Questions
1

Hi, I am new to Pymc3, I am trying to sample from the posterior of a 1d parameter(phi) using Metropolis. The prior of the parameter is normal and the likelihood is a non-standard distribution. The likelihood is calculated on a vector of 1d residuals. I was following this discussion How to set up a custom likelihood function for two variables - #4 by junpenglao.

np.random.seed(123)

def Likelihood(phi,d,value):
	logp = 0.0;
	for val in value:
		logp = logp + (d/2-1)*np.log(val) - 0.5*func1(val, phi)
	return logp
		
def func1(value, phi):
	if value <= phi:
		return value
	else :
		return 2*np.sqrt(phi*value) - phi


data = np.loadtxt("residualsnew.txt", dtype=np.float64)


with Model() as func_model:
# Priors for unknown model parameters
	phi = Normal('phi', mu=3.0, sd=2.0)
	like = DensityDist('like',Likelihood, observed=dict(phi = phi, d = 1, value = data))
        #step1 = find_MAP(model = Huber_model)                                               
	step2 = Metropolis([phi])
	trace = sample(10000, step= [step2], tune = true, tune_interval = 500, chains = 5)

The error I am getting is

    if value <= phi:
  File "/home/navlab-shounak/.local/lib/python3.8/site-packages/theano/tensor/var.py", line 67, in __bool__
    raise TypeError("Variables do not support boolean operations.")
TypeError: Variables do not support boolean operations.
2

You’ll need to use the conditionals made available in theano.tensor . You will also likely need to use scan rather than a python for loop.

3

Do I replace value <= phi by tt.le(value, phi) ?

4

For your func1(), you will probably end up using tt.switch() instead of the if/then. For your likelihood(), you can probably just calculate the logp of every datapoint in value and then tt.sum() the logps. No need to loop (I don’t think).

5

I changed that, but still getting this error

File "/home/pymc3 codes/pyenv/lib/python3.8/site-packages/theano/graph/fg.py", line 342, in import_var
    raise MissingInputError("Undeclared input", variable=var)
theano.graph.fg.MissingInputError: Undeclared input

I saw this was a version problem with pymc3 and theano. Currently I have
python 3.8.12
pymc3 3.11.2
theano-pymc 1.1.2

6

Can you post the code that is raising this error?

7 
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import arviz as az

# Intialize random number generator
np.random.seed(123)

def Likelihood(phi,d,value):
	return (d/2-1)*pm.math.log(value) - 0.5*func1(value, phi)
		
def func1(value, phi):
	b = pm.math.switch(value <= phi, value, 2*pm.math.sqrt(phi*value) - phi)
	return b

data = np.loadtxt("residualsnew.txt", dtype=np.float64)

with pm.Model() as model:
	phi = pm.Normal('phi', mu=4.0, sd=1.0)
	like = pm.DensityDist('like',Likelihood, observed=dict(phi = phi, d = 1, value = data))
	trace = pm.sample(10000, return_inferencedata=False)
	
with model:
    az.plot_trace(trace, var_names=['phi']);
plt.show()
8

It looks like you may be mixing the Potential and DensityDist approaches outlined in Junpeng’s original post. Changing this line:

like = pm.DensityDist('like',
                      Likelihood,
                      observed=dict(phi = phi,
                                    d = 1,
                                    value = data)
                     )

to this:

like = pm.Potential('like', Likelihood(phi, 1, data))

at least gets it sampling. Is that what you intended?

1 Like
9

Yes,thank you ! It is working now.

1 Like