Custom likelihood function (again) with DensityDist and Theano function

Questions
1

Hi everyone,
I need to fit some experimental data using a model that can be expressed as: y(t)=N(S_0.f(t, \theta),\sigma^2 where S_0 is a scaling parameter, f(t, \theta) is a non linear function of the time t and of the model parameters \theta, and \sigma is the standard deviation of the noise.
Using non informative priors on S_0 and \sigma^2, the integration of the likelihood over S_0 and \sigma^2 yields the marginalized likelihood, which is expressed as:
p(\mathbf{y} \vert \theta) \propto [\mathbf{y}^T.\mathbf{y} - (\mathbf{y}^T.\mathbf{f})^2/(\mathbf{f}^T.\mathbf{f})]^{-N/2}
\mathbf{y} is the vector containing the N observations of y(t), sampled at time t, and \mathbf{f} is the vector containing the model f(t,\theta) sampled at the same time t for a given set of parameters \theta.

To keep it simple, lets assume that f is a mono exponential decay model, where \theta is the characteristic time:
f(t,\theta)=(1-\exp(-t/\theta)).
Furthermore, lets assume that the prior on \theta is normal prior with large standard deviation: p(\theta)\ \sim N(0,100).
How can I build such a custom likelihood function? I have been looking on this forum and on the documentation of PyMC3, so I guess that the solution will use the function DensityDist and an implementation of the function f using Theano. Here is my code:

import pymc3 as pm
import numpy as np
from matplotlib import pyplot as plt
import theano.tensor as T

def logp(theta, time, y):
    f = (1 - T.exp(-time/ theta))
    yty = T.dot(y, y)
    ytf = T.dot(y, f)
    ftf = T.dot(f, f)
    N = T.shape(time)[0]
    return -N / 2 * (yty - ytf ** 2 / ftf)

if __name__ == '__main__':

    # Generate noisy data
    S0=10
    sigma=1
    theta=3
    time=np.linspace(0,9,10)
    mu=S0*(1-np.exp(-time/theta))
    y=np.random.normal(mu,sigma)

    with pm.Model() as model:
        theta=pm.Normal('theta',mu=0, sd=100)
        likelihood=pm.DensityDist('Likelihood', logp, observed={'theta':theta, 'time':time, 'y':y})

        # trace = pm.sample(1000, step=pm.Metropolis())
        trace = pm.sample(1000)

    pm.traceplot(trace)
    plt.show()

The code throw an error if I use the default NUTS algorithm
ValueError: Bad initial energy: nan. The model might be misspecified.
“”"

The above exception was the direct cause of the following exception:

ValueError: Bad initial energy: nan. The model might be misspecified.

The above exception was the direct cause of the following exception:

If I change the sampling method to Metropolis, there is no error but the analysis of the trace shows that \theta is constant on the full length of the chain.
So now I am stuck and I don’t know how to go further.
Could someone explain me how I can build my likelihood and use it in PyMC3?

Thanks in advance,

Best

2

The error would be improve in the next release, but currently you need to do some manual check by hand when you see

For example, here when you do

In [12]: model.check_test_point()
Out[12]: 
theta        -5.52
Likelihood     NaN
Name: Log-probability of test_point, dtype: float64

Which means there is some problem with your logp function, now you should check the implementation:

In [13]: f = (1 - T.exp(-time/ theta))

In [14]: f.tag.test_value
Out[14]: array([nan,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])

So there is a problem of f, the problem is that if theta is zero you cannot use it as a denominator. You should restrict theta to a non-zero real number.

3

Thanks, it works now thanks to your advice!
Additionally, I choose a HalfNormal prior on theta and everything works fine!
I modified the code like this to add a condition on the theta parameter.

import pymc3 as pm
import numpy as np
from matplotlib import pyplot as plt
import theano.tensor as T

def logp(theta, time, y):
    N = y.size
    yty = T.dot(y,y)
    if T.gt(theta, 0):
        f = (1 - T.exp(-time / theta))
        ytf = T.dot(y, f)
        ftf = T.dot(f, f)
        logp = -N / 2 * (yty - ytf ** 2 / ftf)
    else:
        logp = -N / 2 * yty

    return logp

if __name__ == '__main__':

    # Generate here noisy data
    S0=10
    sigma=0.1
    theta=3
    time=np.linspace(0,9,10)
    mu=S0*(1-np.exp(-time/theta))
    y=np.random.normal(mu, sigma)

    with pm.Model() as model:
        theta=pm.HalfNormal('theta', sd=100)
        likelihood=pm.DensityDist('Likelihood', logp, observed=dict(theta=theta,
                                                                    time=time,
                                                                    y=y))
        trace = pm.sample(10000, cores=2)

    pm.traceplot(trace)
    plt.show()

I have an additional question regarding the if statement, I’ve seen that there is a switch and ifelse function in theano (here). Would you recommend using it instead of what I implemented?

Thanks!

4

I always use switch, but I havent test the differences personally so I dont know as well.

5

@tboutelier
Hello,
I am discovering bayesian analysis and I read your post with interest. I plan to write a custom likelihood function. I have not exactly the same likelihood as yours but it shares some similarities. And, I wonder why, in your code, the logp return is -N / 2 * (yty - ytf ** 2 / ftf) and not -N / 2 * T.log(yty - ytf ** 2 / ftf)
Best regards
Stephane

1 Like