1 2from functools import partial 3 4import numpy as np 5from scipy import stats 6 7import statsmodels.api as sm 8from statsmodels.base.model import GenericLikelihoodModel 9from statsmodels.tools.numdiff import approx_fprime, approx_hess 10 11data = sm.datasets.spector.load() 12data.exog = sm.add_constant(data.exog, prepend=False) 13# in this dir 14 15probit_mod = sm.Probit(data.endog, data.exog) 16probit_res = probit_mod.fit() 17loglike = probit_mod.loglike 18score = probit_mod.score 19mod = GenericLikelihoodModel(data.endog, data.exog*2, loglike, score) 20res = mod.fit(method="nm", maxiter = 500) 21 22def probitloglike(params, endog, exog): 23 """ 24 Log likelihood for the probit 25 """ 26 q = 2*endog - 1 27 X = exog 28 return np.add.reduce(stats.norm.logcdf(q*np.dot(X,params))) 29 30 31model_loglike = partial(probitloglike, endog=data.endog, exog=data.exog) 32mod = GenericLikelihoodModel(data.endog, data.exog, loglike=model_loglike) 33res = mod.fit(method="nm", maxiter=500) 34print(res) 35 36 37np.allclose(res.params, probit_res.params, rtol=1e-4) 38print(res.params, probit_res.params) 39 40#datal = sm.datasets.longley.load() 41datal = sm.datasets.ccard.load() 42datal.exog = sm.add_constant(datal.exog, prepend=False) 43# Instance of GenericLikelihood model does not work directly, because loglike 44# cannot get access to data in self.endog, self.exog 45 46nobs = 5000 47rvs = np.random.randn(nobs,6) 48datal.exog = rvs[:,:-1] 49datal.exog = sm.add_constant(datal.exog, prepend=False) 50datal.endog = 1 + rvs.sum(1) 51 52show_error = False 53show_error2 = 1#False 54if show_error: 55 def loglike_norm_xb(self, params): 56 beta = params[:-1] 57 sigma = params[-1] 58 xb = np.dot(self.exog, beta) 59 return stats.norm.logpdf(self.endog, loc=xb, scale=sigma) 60 61 mod_norm = GenericLikelihoodModel(datal.endog, datal.exog, loglike_norm_xb) 62 res_norm = mod_norm.fit(method="nm", maxiter = 500) 63 64 print(res_norm.params) 65 66if show_error2: 67 def loglike_norm_xb(params, endog, exog): 68 beta = params[:-1] 69 sigma = params[-1] 70 #print exog.shape, beta.shape 71 xb = np.dot(exog, beta) 72 #print xb.shape, stats.norm.logpdf(endog, loc=xb, scale=sigma).shape 73 return stats.norm.logpdf(endog, loc=xb, scale=sigma).sum() 74 75 model_loglike3 = partial(loglike_norm_xb, 76 endog=datal.endog, exog=datal.exog) 77 mod_norm = GenericLikelihoodModel(datal.endog, datal.exog, model_loglike3) 78 res_norm = mod_norm.fit(start_params=np.ones(datal.exog.shape[1]+1), 79 method="nm", maxiter = 5000) 80 81 print(res_norm.params) 82 83class MygMLE(GenericLikelihoodModel): 84 # just for testing 85 def loglike(self, params): 86 beta = params[:-1] 87 sigma = params[-1] 88 xb = np.dot(self.exog, beta) 89 return stats.norm.logpdf(self.endog, loc=xb, scale=sigma).sum() 90 91 def loglikeobs(self, params): 92 beta = params[:-1] 93 sigma = params[-1] 94 xb = np.dot(self.exog, beta) 95 return stats.norm.logpdf(self.endog, loc=xb, scale=sigma) 96 97mod_norm2 = MygMLE(datal.endog, datal.exog) 98#res_norm = mod_norm.fit(start_params=np.ones(datal.exog.shape[1]+1), method="nm", maxiter = 500) 99res_norm2 = mod_norm2.fit(start_params=[1.]*datal.exog.shape[1]+[1], method="nm", maxiter = 500) 100np.allclose(res_norm.params, res_norm2.params) 101print(res_norm2.params) 102 103res2 = sm.OLS(datal.endog, datal.exog).fit() 104start_params = np.hstack((res2.params, np.sqrt(res2.mse_resid))) 105res_norm3 = mod_norm2.fit(start_params=start_params, method="nm", maxiter = 500, 106 retall=0) 107print(start_params) 108print(res_norm3.params) 109print(res2.bse) 110print(res_norm3.bse) 111print('llf', res2.llf, res_norm3.llf) 112 113bse = np.sqrt(np.diag(np.linalg.inv(res_norm3.model.hessian(res_norm3.params)))) 114res_norm3.model.score(res_norm3.params) 115 116#fprime in fit option cannot be overwritten, set to None, when score is defined 117# exception is fixed, but I do not think score was supposed to be called 118 119res_bfgs = mod_norm2.fit(start_params=start_params, method="bfgs", fprime=None, 120 maxiter=500, retall=0) 121 122hb=-approx_hess(res_norm3.params, mod_norm2.loglike, epsilon=-1e-4) 123hf=-approx_hess(res_norm3.params, mod_norm2.loglike, epsilon=1e-4) 124hh = (hf+hb)/2. 125print(np.linalg.eigh(hh)) 126 127grad = -approx_fprime(res_norm3.params, mod_norm2.loglike, epsilon=-1e-4) 128print(grad) 129gradb = -approx_fprime(res_norm3.params, mod_norm2.loglike, epsilon=-1e-4) 130gradf = -approx_fprime(res_norm3.params, mod_norm2.loglike, epsilon=1e-4) 131print((gradb+gradf)/2.) 132 133print(res_norm3.model.score(res_norm3.params)) 134print(res_norm3.model.score(start_params)) 135mod_norm2.loglike(start_params/2.) 136print(np.linalg.inv(-1*mod_norm2.hessian(res_norm3.params))) 137print(np.sqrt(np.diag(res_bfgs.cov_params()))) 138print(res_norm3.bse) 139 140print("MLE - OLS parameter estimates") 141print(res_norm3.params[:-1] - res2.params) 142print("bse diff in percent") 143print((res_norm3.bse[:-1] / res2.bse)*100. - 100) 144 145''' 146Optimization terminated successfully. 147 Current function value: 12.818804 148 Iterations 6 149Optimization terminated successfully. 150 Current function value: 12.818804 151 Iterations: 439 152 Function evaluations: 735 153Optimization terminated successfully. 154 Current function value: 12.818804 155 Iterations: 439 156 Function evaluations: 735 157<statsmodels.model.LikelihoodModelResults object at 0x02131290> 158[ 1.6258006 0.05172931 1.42632252 -7.45229732] [ 1.62581004 0.05172895 1.42633234 -7.45231965] 159Warning: Maximum number of function evaluations has been exceeded. 160[ -1.18109149 246.94438535 -16.21235536 24.05282629 -324.80867176 161 274.07378453] 162Warning: Maximum number of iterations has been exceeded 163[ 17.57107 -149.87528787 19.89079376 -72.49810777 -50.06067953 164 306.14170418] 165Optimization terminated successfully. 166 Current function value: 506.488765 167 Iterations: 339 168 Function evaluations: 550 169[ -3.08181404 234.34702702 -14.99684418 27.94090839 -237.1465136 170 284.75079529] 171[ -3.08181304 234.34701361 -14.99684381 27.94088692 -237.14649571 172 274.6857294 ] 173[ 5.51471653 80.36595035 7.46933695 82.92232357 199.35166485] 174llf -506.488764864 -506.488764864 175Optimization terminated successfully. 176 Current function value: 506.488765 177 Iterations: 9 178 Function evaluations: 13 179 Gradient evaluations: 13 180(array([ 2.41772580e-05, 1.62492628e-04, 2.79438138e-04, 181 1.90996240e-03, 2.07117946e-01, 1.28747174e+00]), array([[ 1.52225754e-02, 2.01838216e-02, 6.90127235e-02, 182 -2.57002471e-04, -5.25941060e-01, -8.47339404e-01], 183 [ 2.39797491e-01, -2.32325602e-01, -9.36235262e-01, 184 3.02434938e-03, 3.95614029e-02, -1.02035585e-01], 185 [ -2.11381471e-02, 3.01074776e-02, 7.97208277e-02, 186 -2.94955832e-04, 8.49402362e-01, -5.20391053e-01], 187 [ -1.55821981e-01, -9.66926643e-01, 2.01517298e-01, 188 1.52397702e-03, 4.13805882e-03, -1.19878714e-02], 189 [ -9.57881586e-01, 9.87911166e-02, -2.67819451e-01, 190 1.55192932e-03, -1.78717579e-02, -2.55757014e-02], 191 [ -9.96486655e-04, -2.03697290e-03, -2.98130314e-03, 192 -9.99992985e-01, -1.71500426e-05, 4.70854949e-06]])) 193[[ -4.91007768e-05 -7.28732630e-07 -2.51941401e-05 -2.50111043e-08 194 -4.77484718e-08 -9.72022463e-08]] 195[[ -1.64845915e-08 -2.87059265e-08 -2.88764568e-07 -6.82121026e-09 196 2.84217094e-10 -1.70530257e-09]] 197[ -4.90678076e-05 -6.71320777e-07 -2.46166110e-05 -1.13686838e-08 198 -4.83169060e-08 -9.37916411e-08] 199[ -4.56753924e-05 -6.50857146e-07 -2.31756303e-05 -1.70530257e-08 200 -4.43378667e-08 -1.75592936e-02] 201[[ 2.99386348e+01 -1.24442928e+02 9.67254672e+00 -1.58968536e+02 202 -5.91960010e+02 -2.48738183e+00] 203 [ -1.24442928e+02 5.62972166e+03 -5.00079203e+02 -7.13057475e+02 204 -7.82440674e+03 -1.05126925e+01] 205 [ 9.67254672e+00 -5.00079203e+02 4.87472259e+01 3.37373299e+00 206 6.96960872e+02 7.69866589e-01] 207 [ -1.58968536e+02 -7.13057475e+02 3.37373299e+00 6.82417837e+03 208 4.84485862e+03 3.21440021e+01] 209 [ -5.91960010e+02 -7.82440674e+03 6.96960872e+02 4.84485862e+03 210 3.43753691e+04 9.37524459e+01] 211 [ -2.48738183e+00 -1.05126925e+01 7.69866589e-01 3.21440021e+01 212 9.37524459e+01 5.23915258e+02]] 213>>> res_norm3.bse 214array([ 5.47162086, 75.03147114, 6.98192136, 82.60858536, 215 185.40595756, 22.88919522]) 216>>> print res_norm3.model.score(res_norm3.params) 217[ -4.90678076e-05 -6.71320777e-07 -2.46166110e-05 -1.13686838e-08 218 -4.83169060e-08 -9.37916411e-08] 219>>> print res_norm3.model.score(start_params) 220[ -4.56753924e-05 -6.50857146e-07 -2.31756303e-05 -1.70530257e-08 221 -4.43378667e-08 -1.75592936e-02] 222>>> mod_norm2.loglike(start_params/2.) 223-598.56178102781314 224>>> print np.linalg.inv(-1*mod_norm2.hessian(res_norm3.params)) 225[[ 2.99386348e+01 -1.24442928e+02 9.67254672e+00 -1.58968536e+02 226 -5.91960010e+02 -2.48738183e+00] 227 [ -1.24442928e+02 5.62972166e+03 -5.00079203e+02 -7.13057475e+02 228 -7.82440674e+03 -1.05126925e+01] 229 [ 9.67254672e+00 -5.00079203e+02 4.87472259e+01 3.37373299e+00 230 6.96960872e+02 7.69866589e-01] 231 [ -1.58968536e+02 -7.13057475e+02 3.37373299e+00 6.82417837e+03 232 4.84485862e+03 3.21440021e+01] 233 [ -5.91960010e+02 -7.82440674e+03 6.96960872e+02 4.84485862e+03 234 3.43753691e+04 9.37524459e+01] 235 [ -2.48738183e+00 -1.05126925e+01 7.69866589e-01 3.21440021e+01 236 9.37524459e+01 5.23915258e+02]] 237>>> print np.sqrt(np.diag(res_bfgs.cov_params())) 238[ 5.10032831 74.34988912 6.96522122 76.7091604 169.8117832 239 22.91695494] 240>>> print res_norm3.bse 241[ 5.47162086 75.03147114 6.98192136 82.60858536 185.40595756 242 22.88919522] 243>>> res_norm3.conf_int 244<bound method LikelihoodModelResults.conf_int of <statsmodels.model.LikelihoodModelResults object at 0x021317F0>> 245>>> res_norm3.conf_int() 246array([[0.96421437, 1.01999835], 247 [0.99251725, 1.04863332], 248 [0.95721328, 1.01246222], 249 [0.97134549, 1.02695393], 250 [0.97050081, 1.02660988], 251 [0.97773434, 1.03290028], 252 [0.97529207, 1.01428874]]) 253 254>>> res_norm3.params 255array([ -3.08181304, 234.34701361, -14.99684381, 27.94088692, 256 -237.14649571, 274.6857294 ]) 257>>> res2.params 258array([ -3.08181404, 234.34702702, -14.99684418, 27.94090839, 259 -237.1465136 ]) 260>>> 261>>> res_norm3.params - res2.params 262Traceback (most recent call last): 263 File "<stdin>", line 1, in <module> 264ValueError: shape mismatch: objects cannot be broadcast to a single shape 265 266>>> res_norm3.params[:-1] - res2.params 267array([ 9.96859735e-07, -1.34122981e-05, 3.72278400e-07, 268 -2.14645839e-05, 1.78919019e-05]) 269>>> 270>>> res_norm3.bse[:-1] - res2.bse 271array([ -0.04309567, -5.33447922, -0.48741559, -0.31373822, -13.94570729]) 272>>> (res_norm3.bse[:-1] / res2.bse) - 1 273array([-0.00781467, -0.06637735, -0.06525554, -0.00378352, -0.06995531]) 274>>> (res_norm3.bse[:-1] / res2.bse)*100. - 100 275array([-0.7814667 , -6.6377355 , -6.52555369, -0.37835193, -6.99553089]) 276>>> np.sqrt(np.diag(np.linalg.inv(res_norm3.model.hessian(res_bfgs.params)))) 277array([ NaN, NaN, NaN, NaN, NaN, NaN]) 278>>> np.sqrt(np.diag(np.linalg.inv(-res_norm3.model.hessian(res_bfgs.params)))) 279array([ 5.10032831, 74.34988912, 6.96522122, 76.7091604 , 280 169.8117832 , 22.91695494]) 281>>> res_norm3.bse 282array([ 5.47162086, 75.03147114, 6.98192136, 82.60858536, 283 185.40595756, 22.88919522]) 284>>> res2.bse 285array([ 5.51471653, 80.36595035, 7.46933695, 82.92232357, 286 199.35166485]) 287>>> 288>>> bse_bfgs = np.sqrt(np.diag(np.linalg.inv(-res_norm3.model.hessian(res_bfgs.params)))) 289>>> (bse_bfgs[:-1] / res2.bse)*100. - 100 290array([ -7.51422527, -7.4858335 , -6.74913633, -7.49275094, -14.8179759 ]) 291>>> hb=-approx_hess(res_bfgs.params, mod_norm2.loglike, epsilon=-1e-4) 292>>> hf=-approx_hess(res_bfgs.params, mod_norm2.loglike, epsilon=1e-4) 293>>> hh = (hf+hb)/2. 294>>> bse_bfgs = np.sqrt(np.diag(np.linalg.inv(-hh))) 295>>> bse_bfgs 296array([ NaN, NaN, NaN, NaN, NaN, NaN]) 297>>> bse_bfgs = np.sqrt(np.diag(np.linalg.inv(hh))) 298>>> np.diag(hh) 299array([ 9.81680159e-01, 1.39920076e-02, 4.98101826e-01, 300 3.60955710e-04, 9.57811608e-04, 1.90709670e-03]) 301>>> np.diag(np.inv(hh)) 302Traceback (most recent call last): 303 File "<stdin>", line 1, in <module> 304AttributeError: 'module' object has no attribute 'inv' 305 306>>> np.diag(np.linalg.inv(hh)) 307array([ 2.64875153e+01, 5.91578496e+03, 5.13279911e+01, 308 6.11533345e+03, 3.33775960e+04, 5.24357391e+02]) 309>>> res2.bse**2 310array([ 3.04120984e+01, 6.45868598e+03, 5.57909945e+01, 311 6.87611175e+03, 3.97410863e+04]) 312>>> bse_bfgs 313array([ 5.14660231, 76.91414015, 7.1643556 , 78.20059751, 314 182.69536402, 22.89885131]) 315>>> bse_bfgs - res_norm3.bse 316array([-0.32501855, 1.88266901, 0.18243424, -4.40798785, -2.71059354, 317 0.00965609]) 318>>> (bse_bfgs[:-1] / res2.bse)*100. - 100 319array([-6.67512508, -4.29511526, -4.0831115 , -5.69415552, -8.35523538]) 320>>> (res_norm3.bse[:-1] / res2.bse)*100. - 100 321array([-0.7814667 , -6.6377355 , -6.52555369, -0.37835193, -6.99553089]) 322>>> (bse_bfgs / res_norm3.bse)*100. - 100 323array([-5.94007812, 2.50917247, 2.61295176, -5.33599242, -1.46197759, 324 0.04218624]) 325>>> bse_bfgs 326array([ 5.14660231, 76.91414015, 7.1643556 , 78.20059751, 327 182.69536402, 22.89885131]) 328>>> res_norm3.bse 329array([ 5.47162086, 75.03147114, 6.98192136, 82.60858536, 330 185.40595756, 22.88919522]) 331>>> res2.bse 332array([ 5.51471653, 80.36595035, 7.46933695, 82.92232357, 333 199.35166485]) 334>>> dir(res_bfgs) 335['__class__', '__delattr__', '__dict__', '__doc__', '__getattribute__', '__hash__', '__init__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__str__', '__weakref__', 'bse', 'conf_int', 'cov_params', 'f_test', 'initialize', 'llf', 'mle_retvals', 'mle_settings', 'model', 'normalized_cov_params', 'params', 'scale', 't', 't_test'] 336>>> res_bfgs.scale 3371.0 338>>> res2.scale 33981083.015420213851 340>>> res2.mse_resid 34181083.015420213851 342>>> print np.sqrt(np.diag(np.linalg.inv(-1*mod_norm2.hessian(res_bfgs.params)))) 343[ 5.10032831 74.34988912 6.96522122 76.7091604 169.8117832 344 22.91695494] 345>>> print np.sqrt(np.diag(np.linalg.inv(-1*res_bfgs.model.hessian(res_bfgs.params)))) 346[ 5.10032831 74.34988912 6.96522122 76.7091604 169.8117832 347 22.91695494] 348 349Is scale a misnomer, actually scale squared, i.e. variance of error term ? 350''' 351 352print(res_norm3.model.score_obs(res_norm3.params).shape) 353 354jac = res_norm3.model.score_obs(res_norm3.params) 355print(np.sqrt(np.diag(np.dot(jac.T, jac)))/start_params) 356jac2 = res_norm3.model.score_obs(res_norm3.params, centered=True) 357 358print(np.sqrt(np.diag(np.linalg.inv(np.dot(jac.T, jac))))) 359print(res_norm3.bse) 360print(res2.bse) 361