# Four different ways to deal with the log-likelihood of a probability density function (Python code included)

This is not really a question but more of a discussion. Please correct me where wrong and share your thoughts and past experience with regards to computing the likelihoods for continuous data models.

In In All Likelihood by Yudi Pawitan, there is a discussion on page 23 that briefly addresses the likelihood of continuous models.

Usually, the likelihood is given as $$L(\theta)= P(X=x)$$. Hence it is the probability of the observed data. It is easy to construct a likelihood when working with discrete data as a probability mass function gives you the probability of that observation.

The same does not hold for continuous data as a probability density function does not give you a probability of observing some given datum. I have seen four difference approaches to amending this.

# Method 1: Standard

Here we accept all the ambiguity involved in using a probability density function instead of a probability mass function. Hence,

$$L(\theta) = \prod_{t=1}^T f(x_t)$$ We will use the log-likelihood to prevent numerical undeflow: $$l(\theta) = \sum_{t=1}^T \ln (f(x_t))$$

Issues:

• Unbounded likelihood problems can occur. E.g. if the datum lies at the mean of a normal distribution with almost zero variance.
• The log-likelihood can be positive: this is not so much an issue as an annoyance. It just seems wrong.
• The log-likelihoods for the same data-set using probability mass functions can be compared because the log-likelihoods will lie in the same vicinity. Hence, the better model can be identified. This is not true for pdfs. If you compare a exponential pdf model and a normal pdf for the same data-set, then the log-likelihoods will lie in very different values of each other. Hence, this approach does not allow different models to be compared. This stems from the fact that likelihoods for pdfs are ambiguous.
• Due to the ambiguity in its definition, one cannot interpret the log-likelihood to have meaning beyond the use of a single model. By that, I mean we can optimize the log-likelihood of a model using this approach. It will give the best fit for the parameter. However, nothing more could reasonably said of it.

# Method 2: Interval Censoring

Hidden Markov Models for Time-series by Walter Zucchini and Ian MacDonald uses interval censoring for some of their HMM's that have pdf for the emission distributions. Here, we accept that data must lie in some fixed intervals, starting from zero e.g. $$0,0.01,0.02,0.03\dots$$

This is simple to implement, acknowledges a sense of uncertainty or limited precision in the measuring instrument, leads to bounded likelihood and has meaningful interpretation.

You get some datum such as $$3.33$$ but have censored your data into intervals of length $$0.1$$. Hence you would compute the lower and upper bound of the interval your datum lies in. Let $$F$$ denote the cumulative density function. The lower and upper bounds are $$\underline{x_i}$$ and $$\overline{x_i}$$. We then have:

$$l(\theta) = \sum_{t=1}^T \ln (F(\overline{x_t}) - F(\underline{x_t}))$$

We interpret this as the probability that a datum lies in a given interval. Our likelihood uses probabilities and is not ambiguous anymore. One can use this approach to compare different models as the log-likelihoods will lie in the same regions. Hence, the better model can be identified on the basis of this interval censored likelihood.

Note:

lower_bound = interval_length*(datum//interval_length)
upper_bound = lower_bound + interval_length


# Method 3: Finite Precision

This is very similar to interval censoring and is used in the Yudi Pawitan book. Finite precision says that a data meaurement has error $$\epsilon$$. Hence that datum should lie in the interval $$[x_t - 0.5\epsilon,x_t+0.5\epsilon]$$. In contrast to interval censoring, the intervals are dynamic and not fixed. This approach might be better for a measuring instrument such as a timer when compared to interval censoring. What do you think?

Hence, we have:

$$l(\theta) = \sum_{t=1}^T \ln (F(x_t+0.5\epsilon) - F(x_t-0.5\epsilon))$$

I like this method as it seems to solve all the problems.

# Method 4: Correct likelihood

On page for of this paper, one can find the "correct likelihood" remedy for our problem (although purely intended to fix the unbounded likelihood problem).

Our log-likelihood can be written as $$l(\theta) = \sum_{t=1}^T \ln ( \frac{F(x_t+0.5\epsilon) - F(x_t-0.5\epsilon)}{\epsilon})$$

I did not understand this as first (I thought it looked like Method 3 but with division) but then realised that it will simply taking the derivative of the cumulative density function using the central finite difference formula in hope that it would be bounded if the epsilon is large enough.

It seems to give similar values to the original pdf. This is to be expected as as it is a finite difference approximation of it with the purpose of being robust against unbounded values if $$\epsilon$$ is large enough.

I see this approach as only fixing the unbounded likelihood issue. It does not fix the other problems. Using this approach, you will get very different log-likelihood values when comparing different distributions on the same data-set.

# Experiment:

I use all four methods on a data-set generated by an exponential distribution with rate 4. All four methods give the same log-likelihood profile shape. A global maximum is found at 4. This means all four methods can be used to fit an exponential distribution via numerical optimization. Note, however, the positive log-likelihoods.

The Python Code to play with is given below. Note, you can alter it to use another pdf instead of say the Exponential and use it on the same data-set. This will show what I mean by Method 3 and 2 generate consistent log-likelihood values for different models near their optimal parameter setting.


import numpy as np
from scipy.stats import expon
import matplotlib.pyplot as plt

#-------------------------------------------------------------------------------

class expon_ll(object):

def __init__(self,data,rate,interval_size):
self.data = data
self.rate = rate
self.interval_size = interval_size

@property
def scale(self):
return 1/self.rate  # scale gets updated if rate does

def standard(self):
'''
Standard log-likelihood computation using a probability density.
However, a density is NOT a probability hence the likelihood is
ambiguous.
'''
ll = 0
for datum in self.data:
ll += np.log(expon.pdf(datum,scale=self.scale))
return ll

def interval_censored(self):
'''
Fixed intervals, starting from zero, are defined. The probability if
being in the relevant fixed interval is computed.
'''
ll = 0
for datum in self.data:
lower_bound = self.interval_size*(datum//self.interval_size)
upper_bound = lower_bound + self.interval_size
# Compute the probability of being in that fixed interval.
prob = expon.cdf(upper_bound,scale=self.scale) -\
expon.cdf(lower_bound,scale=self.scale)
ll += np.log(prob)
return ll

def finite_precision(self):

ll = 0
for datum in self.data:
lower_bound = datum - 0.5*self.interval_size
# We know the data to be positive
# Hence the lower bound cannot be negative.
# If it is lower than zero, we have gained certainty
if lower_bound < 0:
lower_bound = 0
upper_bound = datum + 0.5*self.interval_size
# Probabiity of being within some margin of error/precision
prob = expon.cdf(upper_bound,scale=self.scale) -\
expon.cdf(lower_bound,scale=self.scale)
ll += np.log(prob)
return ll

def correct(self):
'''
The so called correct likelihood proposed to fix the unbounded
likelihood problem. It is interpreted to be a finite difference
formula to approximate the unbounded pdf by the derivative of the
cdf at the point of the datum. Central finite differences are preferred,
howver, we know the data to be positive. Hence resort to forward
differences if the lower bound of the precision interval is negative.
'''

ll = 0

# Helper functions
def central_diff(f_plus,f_minus,delta):
return (f_plus-f_minus)/(2*delta)

def forward_diff(f_plus,f,delta):
return (f_plus-f)/delta

# For the finite differences
delta = 0.5*self.interval_size

for datum in self.data:
lower_bound = datum - delta
# We know the data to be positive
# Hence the lower bound cannot be negative.
# If it is, truncate at zero
derivative = None
# Use forward differences
if lower_bound < 0:
lower_bound = 0
upper_bound = datum + delta
f_plus = expon.cdf(upper_bound,scale=self.scale)
f = expon.cdf(datum,scale=self.scale)
derivative = forward_diff(f_plus,f,delta)
# Use preferred central differnces
else:
upper_bound = datum + delta
f_plus = expon.cdf(upper_bound,scale=self.scale)
f_minus = expon.cdf(lower_bound,scale=self.scale)
derivative = central_diff(f_plus,f_minus,delta)

ll += np.log(derivative)
return ll

#------------------------------------------------------------------------------

if __name__ == '__main__':

# Create some data
rate = 3
epsilon = 1e-1
np.random.seed()
X = np.random.exponential(scale=1/rate,size=100)

# Record likelihood profile over different rate parameters for
# each of the different methods
N = 100
RATES = np.linspace(0.01,5,N)
log_likelihoods = np.zeros((4,N))  # Record here
ll = expon_ll(X,rate,epsilon)    # Class instance
for n,rate in enumerate(RATES):
ll.rate = rate    # only need to update the rate of the class instance
log_likelihoods[0,n] = ll.standard()
log_likelihoods[1,n] = ll.interval_censored()
log_likelihoods[2,n] = ll.correct()
log_likelihoods[3,n] = ll.finite_precision()

# PLOT:

fig,ax = plt.subplots(2,2)

ax[0,0].plot(RATES,log_likelihoods[0])
ax[0,0].grid(True)
ax[0,0].set_title('Standard (use density)')

ax[1,0].plot(RATES,log_likelihoods[1])
ax[1,0].grid(True)
ax[1,0].set_title('Interval censored')

ax[0,1].plot(RATES,log_likelihoods[2])
ax[0,1].grid(True)
ax[0,1].set_title('\"Correct\" likelihood')

ax[1,1].plot(RATES,log_likelihoods[3])
ax[1,1].grid(True)
ax[1,1].set_title('Finite precision')
plt.tight_layout()
plt.show()



# Conclusion:

Interval censoring and finite precision would seem to be the best way of evaluating the log-likelihood for continuous data models as one gets an almost equivalent probabiity mass function. This results in unambiguous log-likelihoods and likelihoods. Models can have their likelihoods interpreted and models can be compared on the basis of their log-likelihoods as they should lie in the same region. Furthermore, you will never get positive log-likelihoods via these two approaches. Unbounded likelihood is also a thing of the past with these two approaches.

I would use the finite-precision method over interval-censoring because it enables one to include the measurement error quoted by your measuring instrument.

Would you agree with this?