CODEWITHSUNDEEP
pythonscipyregression
user14241
user14241

Reputation: 748

how can I do a maximum likelihood regression using scipy.optimize.minimize

How can I do a maximum likelihood regression using scipy.optimize.minimize? I specifically want to use the minimize function here, because I have a complex model and need to add some constraints. I am currently trying a simple example using the following:

from scipy.optimize import minimize

def lik(parameters):
    m = parameters[0]
    b = parameters[1]
    sigma = parameters[2]
    for i in np.arange(0, len(x)):
        y_exp = m * x + b
    L = sum(np.log(sigma) + 0.5 * np.log(2 * np.pi) + (y - y_exp) ** 2 / (2 * sigma ** 2))
    return L

x = [1,2,3,4,5]
y = [2,3,4,5,6]
lik_model = minimize(lik, np.array([1,1,1]), method='L-BFGS-B', options={'disp': True})

When I run this, convergence fails. Does anyone know what is wrong with my code?

The message I get running this is 'ABNORMAL_TERMINATION_IN_LNSRCH'. I am using the same algorithm that I have working using optim in R.

Upvotes: 11

Views: 23496

Answers (2)

JeeyCi
JeeyCi

Reputation: 589

Better separate your model func. & objective func. to minimize SSE - sum of squared errors (or residuals) - in order to use with any other model. If you would like further to implement your stuff in the class for OOP - see S.O.L.I.D. principles here or here

Also, can use log-likelihood in the form of norm.logpdf(y, mu,sigma).sum() from scipy.stats rather than your manually defined Normal Law of residual's distribution.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Define a dictionary containing employee data
d = {'maturity':[0.75, 2.5, 3., 5., 7.5, 7.9, 8., 12., 16.],
        'price':[967.26, 834.21, 810.52, 769.30, 656.64, 639.71, 604.61, 502.11, 393.38]}

# Convert the dictionary into DataFrame
df = pd.DataFrame(d)
print(df)

x= df.maturity
y= df.price

from scipy.optimize import minimize
from scipy.stats import norm

x0= [0., 1., 1.]

# model
def mod(x, r, b):
    return 1000* np.exp(-1*r*x) + b

# loss-function
def ll(theta):
    r, b = theta[:-1]
##    sigma = theta[-1]
##    mu = np.dot(x,beta)
##    ll = -N/2 * np.log(2*np.pi*sigma**2) - (1/(2*sigma**2)) * np.sum((y - mu)**2)
##    return -1 * ll
    sigma = np.exp(theta[-1])
    for i in np.arange(0, len(x)):
        mu = mod(x, r, b)
    ll = norm.logpdf(y, mu,sigma).sum()     # work with this version of log-likelihood rather than the one above
    return -1 * ll


fit = minimize(ll, np.array(x0), method ='BFGS' )
print(fit.x)
##Formula: price ~ 1000 * exp(-r * maturity)

plt.scatter(x,y)
plt.plot(x, [mod(i, fit.x[0], fit.x[1]) for i in x])
plt.xlabel('time till maturity (maturity date)')
plt.ylabel('bond_price')
plt.show()

enter image description here

Nevertheless, scipy.optimize.minimize not always works good for regression -- depending on init_x0, method chosen and data itself, methods with jacobian and hessian can increase the speed of convergence

P.S. for exponential-fit with LS better linearize xy relations - e.g. last code - though there still can exist "log of zero" problem - pay attention to estimator

Upvotes: 0

user14241
user14241

Reputation: 748

Thank you Aleksander. You were correct that my likelihood function was wrong, not the code. Using a formula I found on wikipedia I adjusted the code to:

import numpy as np
from scipy.optimize import minimize

def lik(parameters):
    m = parameters[0]
    b = parameters[1]
    sigma = parameters[2]
    for i in np.arange(0, len(x)):
        y_exp = m * x + b
    L = (len(x)/2 * np.log(2 * np.pi) + len(x)/2 * np.log(sigma ** 2) + 1 /
         (2 * sigma ** 2) * sum((y - y_exp) ** 2))
    return L

x = np.array([1,2,3,4,5])
y = np.array([2,5,8,11,14])
lik_model = minimize(lik, np.array([1,1,1]), method='L-BFGS-B')
plt.scatter(x,y)
plt.plot(x, lik_model['x'][0] * x + lik_model['x'][1])
plt.show()

Now it seems to be working.

maximum likelihood regression

Thanks for the help!

Upvotes: 11

Related Questions