Skip to content

Commit

Permalink
Merge 0f24c10 into 0529da7
Browse files Browse the repository at this point in the history
  • Loading branch information
jseabold committed Sep 21, 2014
2 parents 0529da7 + 0f24c10 commit e59d31a
Show file tree
Hide file tree
Showing 5 changed files with 40 additions and 19 deletions.
24 changes: 15 additions & 9 deletions statsmodels/base/model.py
Expand Up @@ -596,7 +596,7 @@ def score(self, params):
kwds.setdefault('centered', True)
return approx_fprime(params, self.loglike, **kwds).ravel()

def jac(self, params, **kwds):
def score_obs(self, params, **kwds):
'''
Jacobian/Gradient of log-likelihood evaluated at params for each
observation.
Expand All @@ -605,6 +605,9 @@ def jac(self, params, **kwds):
kwds.setdefault('centered', True)
return approx_fprime(params, self.loglikeobs, **kwds)

jac = np.deprecate(score_obs, 'jac', 'score_obs', "Use score_obs method."
" jac will be removed in 0.7.")

def hessian(self, params):
'''
Hessian of log-likelihood evaluated at params
Expand Down Expand Up @@ -1609,10 +1612,14 @@ def bic(self):
return -2 * self.llf + np.log(self.nobs) * (self.df_modelwc)

@cache_readonly
def jacv(self):
def score_obsv(self):
'''cached Jacobian of log-likelihood
'''
return self.model.jac(self.params)
return self.model.score_obs(self.params)

jacv = np.deprecate(score_obsv, 'jacv', 'score_obsv',
"Use score_obsv attribute."
" jacv will be removed in 0.7.")

@cache_readonly
def hessv(self):
Expand All @@ -1631,7 +1638,7 @@ def covjac(self):
## raise ValueError('need to call fit first')
## #self.fit()
## self.jacv = jacv = self.jac(self._results.params)
jacv = self.jacv
jacv = self.score_obsv
return np.linalg.inv(np.dot(jacv.T, jacv))

@cache_readonly
Expand All @@ -1642,11 +1649,10 @@ def covjhj(self):
name should be covhjh
'''
jacv = self.jacv
## hessv = self.hessv
## hessinv = np.linalg.inv(hessv)
## self.hessinv = hessinv
hessinv = self.cov_params()
jacv = self.score_obsv
hessv = self.hessv
hessinv = np.linalg.inv(hessv)
## self.hessinv = hessin = self.cov_params()
return np.dot(hessinv, np.dot(np.dot(jacv.T, jacv), hessinv))

@cache_readonly
Expand Down
25 changes: 20 additions & 5 deletions statsmodels/discrete/discrete_model.py
Expand Up @@ -1056,7 +1056,7 @@ def score(self, params):
L = np.exp(np.dot(X,params) + offset + exposure)
return np.dot(self.endog - L, X)

def jac(self, params):
def score_obs(self, params):
"""
Poisson model Jacobian of the log-likelihood for each observation
Expand Down Expand Up @@ -1086,6 +1086,9 @@ def jac(self, params):
L = np.exp(np.dot(X,params) + offset + exposure)
return (self.endog - L)[:,None] * X

jac = np.deprecate(score_obs, 'jac', 'score_obs', "Use score_obs method."
" jac will be removed in 0.7")

def hessian(self, params):
"""
Poisson model Hessian matrix of the loglikelihood
Expand Down Expand Up @@ -1257,7 +1260,7 @@ def score(self, params):
L = self.cdf(np.dot(X,params))
return np.dot(y - L,X)

def jac(self, params):
def score_obs(self, params):
"""
Logit model Jacobian of the log-likelihood for each observation
Expand Down Expand Up @@ -1285,6 +1288,9 @@ def jac(self, params):
L = self.cdf(np.dot(X, params))
return (y - L)[:,None] * X

jac = np.deprecate(score_obs, 'jac', 'score_obs', "Use score_obs method."
" jac will be removed in 0.7")

def hessian(self, params):
"""
Logit model Hessian matrix of the log-likelihood
Expand Down Expand Up @@ -1465,7 +1471,7 @@ def score(self, params):
L = q*self.pdf(q*XB)/np.clip(self.cdf(q*XB), FLOAT_EPS, 1 - FLOAT_EPS)
return np.dot(L,X)

def jac(self, params):
def score_obs(self, params):
"""
Probit model Jacobian for each observation
Expand Down Expand Up @@ -1497,6 +1503,9 @@ def jac(self, params):
L = q*self.pdf(q*XB)/np.clip(self.cdf(q*XB), FLOAT_EPS, 1 - FLOAT_EPS)
return L[:,None] * X

jac = np.deprecate(score_obs, 'jac', 'score_obs', "Use score_obs method."
" jac will be removed in 0.7")

def hessian(self, params):
"""
Probit model Hessian matrix of the log-likelihood
Expand Down Expand Up @@ -1710,7 +1719,7 @@ def loglike_and_score(self, params):
score_array = np.dot(firstterm.T, self.exog).flatten()
return loglike_value, score_array

def jac(self, params):
def score_obs(self, params):
"""
Jacobian matrix for multinomial logit model log-likelihood
Expand Down Expand Up @@ -1741,6 +1750,9 @@ def jac(self, params):
#NOTE: might need to switch terms if params is reshaped
return (firstterm[:,:,None] * self.exog[:,None,:]).reshape(self.exog.shape[0], -1)

jac = np.deprecate(score_obs, 'jac', 'score_obs', "Use score_obs method."
" jac will be removed in 0.7")

def hessian(self, params):
"""
Multinomial logit Hessian matrix of the log-likelihood
Expand Down Expand Up @@ -2180,10 +2192,13 @@ def _hessian_nb2(self, params):
return hess_arr

#TODO: replace this with analytic where is it used?
def jac(self, params):
def score_obs(self, params):
sc = approx_fprime_cs(params, self.loglikeobs)
return sc

jac = np.deprecate(score_obs, 'jac', 'score_obs', "Use score_obs method."
" jac will be removed in 0.7")

def fit(self, start_params=None, method='bfgs', maxiter=35,
full_output=1, disp=1, callback=None,
cov_type='nonrobust', cov_kwds=None, use_t=None, **kwargs):
Expand Down
2 changes: 1 addition & 1 deletion statsmodels/discrete/tests/test_discrete.py
Expand Up @@ -111,7 +111,7 @@ def test_loglikeobs(self):

def test_jac(self):
#basic cross check
jacsum = self.res1.model.jac(self.res1.params).sum(0)
jacsum = self.res1.model.score_obs(self.res1.params).sum(0)
score = self.res1.model.score(self.res1.params)
assert_almost_equal(jacsum, score, DECIMAL_9) #Poisson has low precision ?

Expand Down
6 changes: 3 additions & 3 deletions statsmodels/examples/ex_generic_mle.py
Expand Up @@ -373,11 +373,11 @@ def loglikeobs(self, params):
Is scale a misnomer, actually scale squared, i.e. variance of error term ?
'''

print(res_norm3.model.jac(res_norm3.params).shape)
print(res_norm3.model.score_obs(res_norm3.params).shape)

jac = res_norm3.model.jac(res_norm3.params)
jac = res_norm3.model.score_obs(res_norm3.params)
print(np.sqrt(np.diag(np.dot(jac.T, jac)))/start_params)
jac2 = res_norm3.model.jac(res_norm3.params, centered=True)
jac2 = res_norm3.model.score_obs(res_norm3.params, centered=True)

print(np.sqrt(np.diag(np.linalg.inv(np.dot(jac.T, jac)))))
print(res_norm3.bse)
Expand Down
2 changes: 1 addition & 1 deletion statsmodels/genmod/tests/test_glm.py
Expand Up @@ -130,7 +130,7 @@ def test_compare_discrete(self):

assert_allclose(res1.llf, resd.llf, rtol=1e-10)
score_obs1 = res1.model.score_obs(res1.params)
score_obsd = resd.model.jac(resd.params)
score_obsd = resd.model.score_obs(resd.params)
assert_allclose(score_obs1, score_obsd, rtol=1e-10)

# score
Expand Down

0 comments on commit e59d31a

Please sign in to comment.