In [1]:
import warnings
import numpy as np
import pandas as pd
import matplotlib.pylab as plt
import seaborn as sns
from scipy.optimize import minimize, check_grad
from pandas.tools.plotting import scatter_matrix
np.set_printoptions(precision = 3, suppress = True)
sns.set_context('notebook')
%matplotlib inline
In [2]:
df = pd.read_csv('../data/attendance.csv').dropna()
df.sort_values(by='Attendance', inplace=True)
Y = df['Exam'].values
X = df['Attendance'].values
# Adjust by one to allow for log(X)
X += 1
print(df.head())
In [3]:
scatter_matrix(df, alpha=0.2, figsize=(6, 6), diagonal='kde')
plt.show()
In [4]:
class Model_Poly(object):
def __init__(self, Y, X, nterms):
self.nterms = nterms
self.Y, self.X = Y, X
self.N = Y.shape[0]
def xall(self, theta):
return np.atleast_2d(self.X).T ** list(range(len(theta)))
def yhat(self, theta):
return self.xall(theta).dot(theta)
def sse(self, theta):
return ((self.Y - self.yhat(theta))**2).mean()
def jac(self, theta):
return -2 * (self.Y - self.yhat(theta)).dot(self.xall(theta)) / self.N
def hess(self, theta):
return 2 * self.xall(theta).T.dot(self.xall(theta)) / self.N
In [5]:
class Model_Power(object):
def __init__(self, Y, X):
self.Y, self.X = Y, X
self.N = Y.shape[0]
def yhat(self, theta):
return theta[0] + theta[1] * X**theta[2]
def dyhat(self, theta):
return np.vstack([np.ones(self.N),
self.X**theta[2],
theta[1] * self.X**theta[2] * np.log(self.X)])
def sse(self, theta):
return ((self.Y - self.yhat(theta))**2).mean()
def jac(self, theta):
return -2 * (self.Y - self.yhat(theta)).dot(self.dyhat(theta).T) / self.N
In [6]:
methods = ['Nelder-Mead', 'Powell', 'CG', 'BFGS',
'Newton-CG', 'L-BFGS-B', 'TNC', 'COBYLA',
'SLSQP', 'dogleg', 'trust-ncg']
class NLS(object):
def __init__(self, model):
self.model = model
if hasattr(self.model, 'jac'):
self.jac = self.model.jac
else:
self.jac = None
if hasattr(self.model, 'hess'):
self.hess = self.model.hess
else:
self.hess = None
def estimate(self, theta_start, method):
try:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
self.res = minimize(self.model.sse, theta_start, method = method,
jac = self.jac, hess = self.hess)
self.theta_hat = self.res.x
except ValueError:
print('Hessian is required!')
In [7]:
nterms = 3
model = Model_Poly(Y, X, nterms)
plt.figure(figsize = (10, 5))
theta_start = np.ones(nterms)
print('The difference between analytic and numerical gradient =',
check_grad(model.sse, model.jac, theta_start))
estim = NLS(model)
for method in methods:
estim.estimate(theta_start, method)
print(method, estim.res.success, estim.theta_hat)
if estim.res.success:
Yhat = model.yhat(estim.theta_hat)
plt.plot(X, Yhat, label = method)
plt.scatter(X, Y)
plt.legend()
plt.show()
In [8]:
nterms = 3
model = Model_Power(Y, X)
plt.figure(figsize = (10, 5))
theta_start = [35, 5, .5]
print('The difference between analytic and numerical gradient =',
check_grad(model.sse, model.jac, theta_start))
estim = NLS(model)
for method in methods:
estim.estimate(theta_start, method)
print(method, estim.res.success, estim.theta_hat)
if estim.res.success:
Yhat = model.yhat(estim.theta_hat)
plt.plot(X, Yhat, label = method)
plt.scatter(X, Y)
plt.legend()
plt.show()
In [9]:
%%time
# Number of bootstrap samples
B = 1000
# Confidence level
alpha = .05
# Matrix of random indices
new_index = np.random.randint(Y.shape[0], size=[B, Y.shape[0]])
nterms = 3
model = Model_Poly(Y, X, nterms)
# Instantiate estimator
estim = NLS(model)
# Initial guess
theta_start = np.ones(nterms)
# Estimate parameters
estim.estimate(theta_start, method)
# Update initial guess for optimization speed
theta_start = estim.theta_hat
theta_b = []
for idx in new_index:
# Resample the data and pass it to the model
model = Model_Poly(Y[idx], X[idx], nterms)
# Instantiate estimator
estim = NLS(model)
# Estimate
estim.estimate(theta_start, method)
# Store the results
theta_b.append(estim.theta_hat)
# Convert the list to numpy array
theta_b = np.sort(np.vstack(theta_b), axis=0)
# Compute confidence intervals
lo = theta_b[np.floor(alpha * B)]
hi = theta_b[np.ceil((1-alpha) * B)]
print('NLS estimator = ', theta_start)
print('Lower bound = ', lo)
print('Upper bound = ', hi)