# Maximum Likelihood Estimation of a dataset

I am coding a Maximum Likelihood Estimation of a given dataset (Data.csv). The goal is to estimate the mean and sigma. Note that the log of the dataset is well approximated by a normal distribution.

While working on the code, I have faced some issues that drive me crazy. The questions are listed below:

• Maximizing the criterion function provides a postive log-likelihood. I expected the log-likelihood to be negative (see also 3D plot). Does a positive log-likelihood make any sense at all?
• Why is the iteration method 'L-BFGS-B' (highly) sensitive to initalizations for mu_init and sig_init, while the method 'SLSQP' is not? Is there a logic behind this? Moreover, I noticed that the method 'L-BFGS-B' even provides inaccurate results when the initialization is far from the 'true' value.
• Is there a way to calculate the Hessian when the method 'SLSQP' is employed?

• I get small values (almost zero) when computing the log-likelihood, is there a way to deal with this issue of log(0)?

The csv-file containg the data can be found here: Data.csv

# ----------------------------------------------------------------------
# 0. Import the necessary libraries
# ----------------------------------------------------------------------
import numpy as np
import pandas as pd
import scipy.stats as sts
import matplotlib.pyplot as plt
import scipy.optimize as opt # minimizing procedure
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# Display settings
pd.set_option('display.width', 400) # display width in size
pd.set_option('display.max_columns', 50) # display columns

df.set_index('ID', inplace=True) # set ID as index

# Data data
Data = df['Data']                                     # Load data
Data = Data.dropna()
Data = np.log(Data)                                   # Take log() transformation

# Check for right answers estimate
print('mean Data (check)= ', Data.mean(), 'std dev. Data (check) = ', Data.std() )
# ----------------------------------------------------------------------
# 2. Define function that generates probabilities from a normal pdf
# ----------------------------------------------------------------------
def norm_pdf(xvals, mu, sigma):
pdf_vals = ((1 / (sigma * np.sqrt(2 * np.pi)) *
np.exp(- (xvals - mu) ** 2 / (2 * sigma ** 2)))) # normal pdf distribution
return pdf_vals
# ----------------------------------------------------------------------
# 3. Define function that provides the log likelihood
# ----------------------------------------------------------------------
def log_lik_norm(xvals, mu, sigma):
# Generates values from a normal distributed pdf
pdf_vals = norm_pdf(xvals, mu, sigma)
# Take log of normal distributed pdf values
ln_pdf_vals = np.log(pdf_vals)
# Summation of normal distributed pdf values (= loglikelihood)
log_lik_val = ln_pdf_vals.sum()
return log_lik_val
# ----------------------------------------------------------------------
# 4. Define criterion function
# ----------------------------------------------------------------------
def crit(params, *args): # Provides the negative log-likelihood
# Parameters
mu, sigma = params
# Arguments
xvals = args
# Log-likelihood
log_lik_val = log_lik_norm(xvals, mu, sigma)
# Maximizing -> NOTE: Minimizing is computationally more stable
neg_log_lik_val = -log_lik_val
return neg_log_lik_val
# ----------------------------------------------------------------------
# 6. Maximization procedure of criterion function
# ----------------------------------------------------------------------
# Initialization of parameters
mu_init = 3
sig_init = 1 # initalization is sensitive for 'L-BFGS-B' -> first guess: method='SLSQP'
params_init = np.array([mu_init, sig_init])
# Arguments
mle_args = (Data)
# Minimizing procedure (constrained)
results_cstr = opt.minimize(crit, params_init, args=(mle_args), method='SLSQP',
bounds=((None, None), (1e-10, None)))
print(results_cstr)
# Maximum Likelihood Estimators
mu_MLE, sig_MLE = results_cstr.x
print('mu_MLE=', mu_MLE, ' sig_MLE=', sig_MLE)
# ----------------------------------------------------------------------
# 7. Plot the results of MLE
# ----------------------------------------------------------------------
# Plot the MLE estimated distribution
dist_Data = np.linspace(0, 10, 500) # generates chunk of data for pdf
plt.plot(dist_Data, norm_pdf(dist_Data, mu_MLE, sig_MLE),
linewidth=2, color='k', label='3: $$\mu$$={},$$\sigma$$={}'.format(mu_MLE,sig_MLE))
plt.legend(loc='upper left')
#
plt.show()

# Print log-likelihoods for comparison
print('MLE log-likelihood 3: ', log_lik_norm(Data, mu_MLE, sig_MLE)) #MLE (clearly maxmized)

print('MLE log-likelihood check: ', log_lik_norm(Data, Data.mean(), Data.std())) #MLE (clearly maxmized)
# ----------------------------------------------------------------------
# 8. Plot log likelihood around M.L. Estimates for mu and sigma
# ----------------------------------------------------------------------
cmap1 = matplotlib.cm.get_cmap('summer')
#
mu_vals = np.linspace(0.75*mu_MLE, 1.25*mu_MLE, 50) # values around mu_MLE
sig_vals = np.linspace(0.75*sig_MLE, 1.25*sig_MLE, 50) # values around sig_MLE
lnlik_vals = np.zeros((50, 50))
for mu_ind in range(50):
for sig_ind in range(50):
lnlik_vals[mu_ind, sig_ind] = \
log_lik_norm(Data, mu_vals[mu_ind],
sig_vals[sig_ind])

mu_mesh, sig_mesh = np.meshgrid(mu_vals, sig_vals)
#
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(sig_mesh, mu_mesh, lnlik_vals, rstride=8,
cstride=1, cmap=cmap1)
ax.set_title('Log likelihood for values of mu and sigma')
ax.set_xlabel(r'$$\sigma$$')
ax.set_ylabel(r'$$\mu$$')
ax.set_zlabel(r'log likelihood')
plt.show()
# ----------------------------------------------------------------------
# 9. Variance-covariance matrix of ML estimates
# ----------------------------------------------------------------------
# Obtain hessian from minimize function
vcv_mle = results_cstr.hess_inv.todense()
# Standard error of mu
stderr_mu_mle = np.sqrt(vcv_mle[0,0])
# Standard error of sigma
stderr_sig_mle = np.sqrt(vcv_mle[1,1])
#
print('VCV(MLE) = ', vcv_mle)
print('Standard error for mu estimate = ', stderr_mu_mle)
print('Standard error for sigma estimate = ', stderr_sig_mle)


-- EDIT --

However, as I run the following piece of code I get two errors when python tries to obtain: hessian_(theta_autograd)

1. KeyError: < class 'pandas.core.series.Series' >
2. TypeError: Can't differentiate w.r.t. type < class 'pandas.core.series.Series' >

The two errors are driving me crazy and I hope somebody might be of help.

The csv-data can be found here: csv-data

import warnings
warnings.filterwarnings('ignore')
from scipy.optimize import minimize
import statsmodels.api as sm
import pandas as pd
import matplotlib.pyplot as plt

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

df.set_index('ID', inplace=True) # set ID as index

# # Plot
plt.scatter(df.X, df.Y, color='red', marker='+')
plt.title('X vs Y', fontsize=15)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

# ----------------------------------------------------------------------
# 0. Data x, y
# ----------------------------------------------------------------------

N = df.shape[0]             # number of observations
#
x_ = df.X                   # plug independent data!
x = np.c_[np.ones(N),x_]    # add constant to independent variables
y = df.Y                    # plug independent data!

# ----------------------------------------------------------------------
# 1. Parameters MLE
# ----------------------------------------------------------------------

# number of parameters
K = 2
# Parameter values
beta = np.random.randn(K)

# ----------------------------------------------------------------------
# 2. Negative Loglikelihood
# ----------------------------------------------------------------------
def neg_loglike(theta):
beta = theta[:-1]
# transform theta[-1] so that sigma > 0
sigma = np.exp(theta[-1])
mu = np.dot(x,beta)
ll = norm.logpdf(y,mu,sigma).sum()
return -1 * ll

# jacobian of neg_loglike
jacobian_  = jacobian(neg_loglike)
# hessian of neg_loglike
hessian_ = hessian(neg_loglike)

# ----------------------------------------------------------------------
# 2. Minimizing Negative Loglikelihood
# ----------------------------------------------------------------------
# initialization
theta_start = np.append(np.zeros(beta.shape[0]),0.0)
# minimization sci-py criterion
res1 = minimize(neg_loglike, theta_start, method = 'BFGS', jac = jacobian_)
print(res1)
# estimated parameters

# # ----------------------------------------------------------------------
# # 3. Results
# # ----------------------------------------------------------------------
# # std errors, calculate the information matrix and using the autograd hessian
information1 = np.transpose(hessian_(theta_autograd)) # <-- GIVES ERROR
# # se1 = np.sqrt(np.diagonal(np.linalg.inv(information1)))
# # #
# # print('se1 :', se1)
#
# # Put Results in a DataFrame
# names = ['beta_'+str(i) for i in range(K)]
# names.append('log(Sigma)')
# results_a['Variable'] = names
# #
# print(results_a)
# #

# Jacobian of ML estimation

# Loglikelihood of ML estimation
print("ML log-likelihood = ", -1*neg_loglike(res1.x))

# ----------------------------------------------------------------------
# 4. Check with Statsmodels
# ----------------------------------------------------------------------
model = sm.OLS(y, x)
results = model.fit()
print(results.summary())


EDIT 2--

def neg_loglike(theta):
beta = theta[:-1]
loc = np.dot(x, beta)
scale = np.exp(theta[-1])
ll = gumbel_r.logpdf(x, loc, scale).sum()
return -1 * ll

• For 1st question: probability density function (pdf) is not probability, so it can > 1 ===> log pdf > 0 and log pdf = log likelihood. Commented Jul 20, 2019 at 14:38
• Thanks! It makes sense. I noticed indeed that pdf > 1, especially when the sigma is very low. Basically squeezes the pdf to the top! Commented Jul 20, 2019 at 14:42

I can suggest a number of improvements that will make any minimization method more stable.

    # Generates values from a normal distributed pdf
pdf_vals = norm_pdf(xvals, mu, sigma)
# Take log of normal distributed pdf values
ln_pdf_vals = np.log(pdf_vals)


you should analytically compute the log of the pdf. The reason is that the exponential term in the pdf is going to very very large for values away from the true values, and hence could underflow. Taking the log of an underflow'd number gives gibberish results. Instead, write code for the log of the pdf.

1. Currently, you are not giving the solver any gradient (jacobian) information. This is extremely helpful for convergence stability and performance. Take a look at this blog post on how to automatically include gradients: https://rlhick.people.wm.edu/posts/mle-autograd.html Autograd also can compute hessian's exactly!

2. Bounds, like on $$\sigma$$ can be sometimes tricky for inference. You can instead choose a new value, $$\psi$$ that is defined as $$\sigma = \exp \psi$$. So now $$\psi$$ can vary over any negative and positive number, but after exponentiating, $$\sigma$$ is still positive. This opens up other minimization routines to use.

• Thank you for the insightful comments! Really helpful, I upvoted. I will have to take a look at the blog for more information regarding the jacobian information. Commented Jul 20, 2019 at 19:39
• @ Cam Davidson Pilon, I implemented the suggestions you made. Unfortunately, I get some errors that are added in the editted part of the question. Would you mind taking a look at that? Commented Aug 6, 2019 at 17:23
• Try replacing with numpy arrays. Example, X = df.X.values Commented Aug 6, 2019 at 18:56
• It worked! Thank you, much appreciated. It almost drove me nuts. You are such a great help! Commented Aug 6, 2019 at 19:40

I got what looks like a good fit of the data in the .csv file to scipy's right-sided Gumbel distribution:

with the following info:

Right Sided Gumbel distribution
http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gumbel_r.html

Fit Statistics for 994 data points:
Negative Two Log Likelihood = 6.3305082874095360E+03
AIC = 6.3345082874095360E+03
AICc (Burnham and Anderson) = 6.3345203963903632E+03

Parameters:
location = 2.5653678690984339E+01
scale = 4.9869185100533118E+00

The probability density function for gumbel_r is:

gumbel_r.pdf(x) = exp(-(x + exp(-x)))


Here is a graphical fitter for gumbel_r:

import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt

f = open('Data.csv')
f.close()

lines = text.split('\n')

data = []
for line in lines[1:-1]:
ID, lineData = line.split(',')
if lineData:
data.append(float(lineData))

P = ss.gumbel_r.fit(data)
rX = np.linspace(min(data), max(data), 50)
rP = ss.gumbel_r.pdf(rX, *P)

plt.hist(data,bins=25, normed=True, color='slategrey')

plt.plot(rX, rP, color='darkturquoise')
plt.show()

• Well done! Really interesting, I intended to take a logarithmic transformation of the data and then employ a normal distribution. I noticed from a QQ-plot, that a log. transformation went actually quite well. However, the Gumbel distribution appears to be a decent alternative. I'll definitely look into the Gumbel distribution. Thanks! Commented Jul 20, 2019 at 20:25
• this is the "right sided" Gumbel distribution, the scipy documentation link is in the answer. Isn't it curious what you can find on the internet these days... Commented Jul 20, 2019 at 22:10
• I'm trying to use the gumbel_r distribution in the log likelihood expression. But I get two errors implementing Autograd MLE(see edit). If you could take a look at it, or know why it is the case that would be much appreciated. Commented Aug 6, 2019 at 17:26
• I personally have never used Autograd_MLE. Commented Aug 6, 2019 at 19:54
• Bear my ignorance, but the problem is that I don't see exactly how the Gumbel_r distribution can be used in a linear regression. What parameters in the Gumbel_r distribution should be optimized for the loglikelihood? I tried something like depicted in EDIT2. Commented Aug 6, 2019 at 20:38