Created
June 22, 2015 18:33
-
-
Save saketkc/dab5f3d8907f62232062 to your computer and use it in GitHub Desktop.
n_obs=200000
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 356 @profile | |
| 357 def _smw_solver(s, A, AtA, B, BI): | |
| 358 """ | |
| 359 Solves the system (s*I + A*B*A') * x = rhs for an arbitrary rhs. | |
| 360 | |
| 361 Parameters | |
| 362 ---------- | |
| 363 s : scalar | |
| 364 See above for usage | |
| 365 A : square symmetric ndarray | |
| 366 See above for usage | |
| 367 AtA : square ndarray | |
| 368 A.T * A | |
| 369 B : square symmetric ndarray | |
| 370 See above for usage | |
| 371 BI : square symmetric ndarray | |
| 372 The inverse of `B`. | |
| 373 | |
| 374 Returns | |
| 375 ------- | |
| 376 A function that takes `rhs` as an input argument and returns a | |
| 377 solution to the linear system defined above. | |
| 378 """ | |
| 379 | |
| 380 # Use SMW identity | |
| 381 5880 49233 8.4 1.4 qmat = BI + AtA / s | |
| 382 5880 3297098 560.7 95.4 qmati = np.linalg.solve(qmat, A.T) | |
| 383 | |
| 384 5880 105434 17.9 3.1 @profile | |
| 385 def solver(rhs): | |
| 386 ql = np.dot(qmati, rhs) | |
| 387 ql = np.dot(A, ql) | |
| 388 rslt = rhs / s - ql / s**2 | |
| 389 return rslt | |
| 390 | |
| 391 5880 4084 0.7 0.1 return solver | |
| Total time: 18.2266 s | |
| File: statsmodels/regression/mixed_linear_model.orig.py | |
| Function: solver at line 384 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 384 @profile | |
| 385 def solver(rhs): | |
| 386 13440 8058252 599.6 44.2 ql = np.dot(qmati, rhs) | |
| 387 13440 9291510 691.3 51.0 ql = np.dot(A, ql) | |
| 388 13440 868379 64.6 4.8 rslt = rhs / s - ql / s**2 | |
| 389 13440 8449 0.6 0.0 return rslt | |
| Total time: 0.18264 s | |
| File: statsmodels/regression/mixed_linear_model.orig.py | |
| Function: _smw_logdet at line 394 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 394 @profile | |
| 395 def _smw_logdet(s, A, AtA, B, BI, B_logdet): | |
| 396 """ | |
| 397 Use the matrix determinant lemma to accelerate the calculation of | |
| 398 the log determinant of s*I + A*B*A'. | |
| 399 | |
| 400 Parameters | |
| 401 ---------- | |
| 402 s : scalar | |
| 403 See above for usage | |
| 404 A : square symmetric ndarray | |
| 405 See above for usage | |
| 406 AtA : square matrix | |
| 407 A.T * A | |
| 408 B : square symmetric ndarray | |
| 409 See above for usage | |
| 410 BI : square symmetric ndarray | |
| 411 The inverse of `B`; can be None if B is singular. | |
| 412 B_logdet : real | |
| 413 The log determinant of B | |
| 414 | |
| 415 Returns | |
| 416 ------- | |
| 417 The log determinant of s*I + A*B*A'. | |
| 418 """ | |
| 419 | |
| 420 1500 1868 1.2 1.0 p = A.shape[0] | |
| 421 1500 8406 5.6 4.6 ld = p * np.log(s) | |
| 422 1500 11205 7.5 6.1 qmat = BI + AtA / s | |
| 423 1500 159363 106.2 87.3 _, ld1 = np.linalg.slogdet(qmat) | |
| 424 1500 1798 1.2 1.0 return B_logdet + ld + ld1 | |
| Total time: 14.2719 s | |
| File: statsmodels/regression/mixed_linear_model.orig.py | |
| Function: _gen_dV_dPar at line 1313 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 1313 @profile | |
| 1314 def _gen_dV_dPar(self, ex_r, solver, group, max_ix=None): | |
| 1315 """ | |
| 1316 A generator that yields the element-wise derivative of the | |
| 1317 marginal covariance matrix with respect to the random effects | |
| 1318 variance and covariance parameters. | |
| 1319 | |
| 1320 ex_r : array-like | |
| 1321 The random effects design matrix | |
| 1322 solver : function | |
| 1323 A function that given x returns V^{-1}x, where V | |
| 1324 is the group's marginal covariance matrix. | |
| 1325 group : scalar | |
| 1326 The group label | |
| 1327 max_ix : integer or None | |
| 1328 If not None, the generator ends when this index | |
| 1329 is reached. | |
| 1330 """ | |
| 1331 | |
| 1332 1560 7171154 4596.9 50.2 axr = solver(ex_r) | |
| 1333 | |
| 1334 # Regular random effects | |
| 1335 1560 937 0.6 0.0 jj = 0 | |
| 1336 1560 4707 3.0 0.0 for j1 in range(self.k_re): | |
| 1337 for j2 in range(j1 + 1): | |
| 1338 if max_ix is not None and jj > max_ix: | |
| 1339 return | |
| 1340 mat_l, mat_r = ex_r[:,j1:j1+1], ex_r[:,j2:j2+1] # Need 2d | |
| 1341 vsl, vsr = axr[:,j1:j1+1], axr[:,j2:j2+1] | |
| 1342 yield jj, mat_l, mat_r, vsl, vsr, j1 == j2 | |
| 1343 jj += 1 | |
| 1344 | |
| 1345 # Variance components | |
| 1346 4620 3834 0.8 0.0 for ky in self._vc_names: | |
| 1347 3120 4155 1.3 0.0 if group in self.exog_vc[ky]: | |
| 1348 3120 1796 0.6 0.0 if max_ix is not None and jj > max_ix: | |
| 1349 60 71 1.2 0.0 return | |
| 1350 3060 2188 0.7 0.0 mat = self.exog_vc[ky][group] | |
| 1351 3060 7078682 2313.3 49.6 axmat = solver(mat) | |
| 1352 3060 2299 0.8 0.0 yield jj, mat, mat, axmat, axmat, True | |
| 1353 3060 2074 0.7 0.0 jj += 1 | |
| Total time: 17.5665 s | |
| File: statsmodels/regression/mixed_linear_model.orig.py | |
| Function: score at line 1356 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 1356 @profile | |
| 1357 def score(self, params, profile_fe=True): | |
| 1358 """ | |
| 1359 Returns the score vector of the profile log-likelihood. | |
| 1360 | |
| 1361 Notes | |
| 1362 ----- | |
| 1363 The score vector that is returned is computed with respect to | |
| 1364 the parameterization defined by this model instance's | |
| 1365 `use_sqrt` attribute. | |
| 1366 """ | |
| 1367 | |
| 1368 23 29 1.3 0.0 if type(params) is not MixedLMParams: | |
| 1369 23 33 1.4 0.0 params = MixedLMParams.from_packed(params, self.k_fe, | |
| 1370 23 18 0.8 0.0 self.k_re, self.use_sqrt, | |
| 1371 23 1530 66.5 0.0 with_fe=False) | |
| 1372 | |
| 1373 23 20 0.9 0.0 if profile_fe: | |
| 1374 23 1894457 82367.7 10.8 params.fe_params = self.get_fe_params(params.cov_re, params.vcomp) | |
| 1375 | |
| 1376 23 27 1.2 0.0 if self.use_sqrt: | |
| 1377 23 15670289 681316.9 89.2 score_fe, score_re, score_vc = self.score_sqrt(params, calc_fe=not profile_fe) | |
| 1378 else: | |
| 1379 score_fe, score_re, score_vc = self.score_full(params, calc_fe=not profile_fe) | |
| 1380 | |
| 1381 23 18 0.8 0.0 if self._freepat is not None: | |
| 1382 score_fe *= self._freepat.fe_params | |
| 1383 score_re *= self._freepat.cov_re[self._freepat._ix] | |
| 1384 score_vc *= self._freepat.vcomp | |
| 1385 | |
| 1386 23 14 0.6 0.0 if profile_fe: | |
| 1387 23 60 2.6 0.0 return np.concatenate((score_re, score_vc)) | |
| 1388 else: | |
| 1389 return np.concatenate((score_fe, score_re, score_vc)) | |
| Total time: 15.5742 s | |
| File: statsmodels/regression/mixed_linear_model.orig.py | |
| Function: score_full at line 1392 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 1392 @profile | |
| 1393 def score_full(self, params, calc_fe): | |
| 1394 """ | |
| 1395 Returns the score with respect to untransformed parameters. | |
| 1396 | |
| 1397 Calculates the score vector for the profiled log-likelihood of | |
| 1398 the mixed effects model with respect to the parameterization | |
| 1399 in which the random effects covariance matrix is represented | |
| 1400 in its full form (not using the Cholesky factor). | |
| 1401 | |
| 1402 Parameters | |
| 1403 ---------- | |
| 1404 params : MixedLMParams or array-like | |
| 1405 The parameter at which the score function is evaluated. | |
| 1406 If array-like, must contain the packed random effects | |
| 1407 parameters (cov_re and vcomp) without fe_params. | |
| 1408 calc_fe : boolean | |
| 1409 If True, calculate the score vector for the fixed effects | |
| 1410 parameters. If False, this vector is not calculated, and | |
| 1411 a vector of zeros is returned in its place. | |
| 1412 | |
| 1413 Returns | |
| 1414 ------- | |
| 1415 score_fe : array-like | |
| 1416 The score vector with respect to the fixed effects | |
| 1417 parameters. | |
| 1418 score_re : array-like | |
| 1419 The score vector with respect to the random effects | |
| 1420 parameters (excluding variance components parameters). | |
| 1421 score_vc : array-like | |
| 1422 The score vector with respect to variance components | |
| 1423 parameters. | |
| 1424 | |
| 1425 Notes | |
| 1426 ----- | |
| 1427 `score_re` is taken with respect to the parameterization in | |
| 1428 which `cov_re` is represented through its lower triangle | |
| 1429 (without taking the Cholesky square root). | |
| 1430 """ | |
| 1431 | |
| 1432 23 27 1.2 0.0 fe_params = params.fe_params | |
| 1433 23 28 1.2 0.0 cov_re = params.cov_re | |
| 1434 23 26 1.1 0.0 vcomp = params.vcomp | |
| 1435 | |
| 1436 23 23 1.0 0.0 try: | |
| 1437 23 426 18.5 0.0 cov_re_inv = np.linalg.inv(cov_re) | |
| 1438 except np.linalg.LinAlgError: | |
| 1439 cov_re_inv = None | |
| 1440 | |
| 1441 23 110 4.8 0.0 score_fe = np.zeros(self.k_fe) | |
| 1442 23 63 2.7 0.0 score_re = np.zeros(self.k_re2) | |
| 1443 23 44 1.9 0.0 score_vc = np.zeros(self.k_vc) | |
| 1444 | |
| 1445 # Handle the covariance penalty. | |
| 1446 23 31 1.3 0.0 if self.cov_pen is not None: | |
| 1447 score_re -= self.cov_pen.grad(cov_re, cov_re_inv) | |
| 1448 | |
| 1449 # Handle the fixed effects penalty. | |
| 1450 23 23 1.0 0.0 if calc_fe and (self.fe_pen is not None): | |
| 1451 score_fe -= self.fe_pen.grad(fe_params) | |
| 1452 | |
| 1453 # resid' V^{-1} resid, summed over the groups (a scalar) | |
| 1454 23 23 1.0 0.0 rvir = 0. | |
| 1455 | |
| 1456 # exog' V^{-1} resid, summed over the groups (a k_fe | |
| 1457 # dimensional vector) | |
| 1458 23 23 1.0 0.0 xtvir = 0. | |
| 1459 | |
| 1460 # exog' V^{_1} exog, summed over the groups (a k_fe x k_fe | |
| 1461 # matrix) | |
| 1462 23 23 1.0 0.0 xtvix = 0. | |
| 1463 | |
| 1464 # V^{-1} exog' dV/dQ_jj exog V^{-1}, where Q_jj is the jj^th | |
| 1465 # covariance parameter. | |
| 1466 23 42 1.8 0.0 xtax = [0.,] * (self.k_re2 + self.k_vc) | |
| 1467 | |
| 1468 # Temporary related to the gradient of log |V| | |
| 1469 23 46 2.0 0.0 dlv = np.zeros(self.k_re2 + self.k_vc) | |
| 1470 | |
| 1471 # resid' V^{-1} dV/dQ_jj V^{-1} resid (a scalar) | |
| 1472 23 89 3.9 0.0 rvavr = np.zeros(self.k_re2 + self.k_vc) | |
| 1473 | |
| 1474 1403 2096 1.5 0.0 for group_ix, group in enumerate(self.group_labels): | |
| 1475 | |
| 1476 1380 91349 66.2 0.6 cov_aug, cov_aug_inv, _ = self._augment_cov_re(cov_re, cov_re_inv, vcomp, group) | |
| 1477 | |
| 1478 1380 2454 1.8 0.0 exog = self.exog_li[group_ix] | |
| 1479 1380 2252 1.6 0.0 ex_r, ex2_r = self._aex_r[group_ix], self._aex_r2[group_ix] | |
| 1480 1380 867894 628.9 5.6 solver = _smw_solver(1., ex_r, ex2_r, cov_aug, cov_aug_inv) | |
| 1481 | |
| 1482 # The residuals | |
| 1483 1380 2550 1.8 0.0 resid = self.endog_li[group_ix] | |
| 1484 1380 1656 1.2 0.0 if self.k_fe > 0: | |
| 1485 1380 46926 34.0 0.3 expval = np.dot(exog, fe_params) | |
| 1486 1380 9635 7.0 0.1 resid = resid - expval | |
| 1487 | |
| 1488 1380 1692 1.2 0.0 if self.reml: | |
| 1489 1380 766214 555.2 4.9 viexog = solver(exog) | |
| 1490 1380 45652 33.1 0.3 xtvix += np.dot(exog.T, viexog) | |
| 1491 | |
| 1492 # Contributions to the covariance parameter gradient | |
| 1493 1380 205120 148.6 1.3 vir = solver(resid) | |
| 1494 4140 12798991 3091.5 82.2 for jj, matl, matr, vsl, vsr, sym in self._gen_dV_dPar(ex_r, solver, group): | |
| 1495 2760 139621 50.6 0.9 dlv[jj] = np.sum(matr * vsl) # trace dot | |
| 1496 2760 3831 1.4 0.0 if not sym: | |
| 1497 dlv[jj] += np.sum(matl * vsr) # trace dot | |
| 1498 | |
| 1499 2760 95452 34.6 0.6 ul = np.dot(vir, matl) | |
| 1500 2760 6467 2.3 0.0 ur = ul.T if sym else np.dot(matr.T, vir) | |
| 1501 2760 8057 2.9 0.1 ulr = np.dot(ul, ur) | |
| 1502 2760 6870 2.5 0.0 rvavr[jj] += ulr | |
| 1503 2760 3043 1.1 0.0 if not sym: | |
| 1504 rvavr[jj] += ulr.T | |
| 1505 | |
| 1506 2760 3121 1.1 0.0 if self.reml: | |
| 1507 2760 411924 149.2 2.6 ul = np.dot(viexog.T, matl) | |
| 1508 2760 4835 1.8 0.0 ur = ul.T if sym else np.dot(matr.T, viexog) | |
| 1509 2760 7859 2.8 0.1 ulr = np.dot(ul, ur) | |
| 1510 2760 9793 3.5 0.1 xtax[jj] += ulr | |
| 1511 2760 3491 1.3 0.0 if not sym: | |
| 1512 xtax[jj] += ulr.T | |
| 1513 | |
| 1514 # Contribution of log|V| to the covariance parameter | |
| 1515 # gradient. | |
| 1516 1380 1815 1.3 0.0 if self.k_re > 0: | |
| 1517 score_re -= 0.5 * dlv[0:self.k_re2] | |
| 1518 1380 1599 1.2 0.0 if self.k_vc > 0: | |
| 1519 1380 8991 6.5 0.1 score_vc -= 0.5 * dlv[self.k_re2:] | |
| 1520 | |
| 1521 1380 7580 5.5 0.0 rvir += np.dot(resid, vir) | |
| 1522 | |
| 1523 1380 1593 1.2 0.0 if calc_fe: | |
| 1524 xtvir += np.dot(exog.T, vir) | |
| 1525 | |
| 1526 23 32 1.4 0.0 fac = self.n_totobs | |
| 1527 23 24 1.0 0.0 if self.reml: | |
| 1528 23 30 1.3 0.0 fac -= self.k_fe | |
| 1529 | |
| 1530 23 23 1.0 0.0 if calc_fe and self.k_fe > 0: | |
| 1531 score_fe += fac * xtvir / rvir | |
| 1532 | |
| 1533 23 25 1.1 0.0 if self.k_re > 0: | |
| 1534 score_re += 0.5 * fac * rvavr[0:self.k_re2] / rvir | |
| 1535 23 26 1.1 0.0 if self.k_vc > 0: | |
| 1536 23 234 10.2 0.0 score_vc += 0.5 * fac * rvavr[self.k_re2:] / rvir | |
| 1537 | |
| 1538 23 30 1.3 0.0 if self.reml: | |
| 1539 23 1327 57.7 0.0 xtvixi = np.linalg.inv(xtvix) | |
| 1540 23 65 2.8 0.0 for j in range(self.k_re2): | |
| 1541 score_re[j] += 0.5 * np.sum(xtvixi.T * xtax[j]) # trace dot | |
| 1542 69 100 1.4 0.0 for j in range(self.k_vc): | |
| 1543 46 762 16.6 0.0 score_vc[j] += 0.5 * np.sum(xtvixi.T * xtax[self.k_re2 + j]) # trace dot | |
| 1544 | |
| 1545 23 25 1.1 0.0 return score_fe, score_re, score_vc | |
| Total time: 15.6697 s | |
| File: statsmodels/regression/mixed_linear_model.orig.py | |
| Function: score_sqrt at line 1548 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 1548 @profile | |
| 1549 def score_sqrt(self, params, calc_fe=True): | |
| 1550 """ | |
| 1551 Returns the score with respect to transformed parameters. | |
| 1552 | |
| 1553 Calculates the score vector with respect to the | |
| 1554 parameterization in which the random effects covariance matrix | |
| 1555 is represented through its Cholesky square root. | |
| 1556 | |
| 1557 Parameters | |
| 1558 ---------- | |
| 1559 params : MixedLMParams or array-like | |
| 1560 The model parameters. If array-like must contain packed | |
| 1561 parameters that are compatible with this model instance. | |
| 1562 calc_fe : boolean | |
| 1563 If True, calculate the score vector for the fixed effects | |
| 1564 parameters. If False, this vector is not calculated, and | |
| 1565 a vector of zeros is returned in its place. | |
| 1566 | |
| 1567 Returns | |
| 1568 ------- | |
| 1569 score_fe : array-like | |
| 1570 The score vector with respect to the fixed effects | |
| 1571 parameters. | |
| 1572 score_re : array-like | |
| 1573 The score vector with respect to the random effects | |
| 1574 parameters (excluding variance components parameters). | |
| 1575 score_vc : array-like | |
| 1576 The score vector with respect to variance components | |
| 1577 parameters. | |
| 1578 """ | |
| 1579 | |
| 1580 23 15667517 681196.4 100.0 score_fe, score_re, score_vc = self.score_full(params, calc_fe=calc_fe) | |
| 1581 23 346 15.0 0.0 params_vec = params.get_packed(use_sqrt=True, with_fe=True) | |
| 1582 | |
| 1583 23 62 2.7 0.0 score_full = np.concatenate((score_fe, score_re, score_vc)) | |
| 1584 23 15 0.7 0.0 scr = 0. | |
| 1585 161 138 0.9 0.0 for i in range(len(params_vec)): | |
| 1586 138 928 6.7 0.0 v = self._lin[i] + 2 * np.dot(self._quad[i], params_vec) | |
| 1587 138 548 4.0 0.0 scr += score_full[i] * v | |
| 1588 23 29 1.3 0.0 score_fe = scr[0:self.k_fe] | |
| 1589 23 31 1.3 0.0 score_re = scr[self.k_fe:self.k_fe + self.k_re2] | |
| 1590 23 23 1.0 0.0 score_vc = scr[self.k_fe + self.k_re2:] | |
| 1591 | |
| 1592 23 14 0.6 0.0 return score_fe, score_re, score_vc | |
| Total time: 38.9766 s | |
| File: statsmodels/regression/mixed_linear_model.orig.py | |
| Function: hessian at line 1594 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 1594 @profile | |
| 1595 def hessian(self, params): | |
| 1596 """ | |
| 1597 Returns the model's Hessian matrix. | |
| 1598 | |
| 1599 Calculates the Hessian matrix for the linear mixed effects | |
| 1600 model with respect to the parameterization in which the | |
| 1601 covariance matrix is represented directly (without square-root | |
| 1602 transformation). | |
| 1603 | |
| 1604 Parameters | |
| 1605 ---------- | |
| 1606 params : MixedLMParams or array-like | |
| 1607 The model parameters at which the Hessian is calculated. | |
| 1608 If array-like, must contain the packed parameters in a | |
| 1609 form that is compatible with this model instance. | |
| 1610 | |
| 1611 Returns | |
| 1612 ------- | |
| 1613 hess : 2d ndarray | |
| 1614 The Hessian matrix, evaluated at `params`. | |
| 1615 """ | |
| 1616 | |
| 1617 1 1 1.0 0.0 if type(params) is not MixedLMParams: | |
| 1618 params = MixedLMParams.from_packed(params, self.k_fe, self.k_re, | |
| 1619 use_sqrt=self.use_sqrt, | |
| 1620 with_fe=True) | |
| 1621 | |
| 1622 1 1 1.0 0.0 fe_params = params.fe_params | |
| 1623 1 2 2.0 0.0 vcomp = params.vcomp | |
| 1624 1 1 1.0 0.0 cov_re = params.cov_re | |
| 1625 1 1 1.0 0.0 if self.k_re > 0: | |
| 1626 cov_re_inv = np.linalg.inv(cov_re) | |
| 1627 else: | |
| 1628 1 4 4.0 0.0 cov_re_inv = np.empty((0, 0)) | |
| 1629 | |
| 1630 # Blocks for the fixed and random effects parameters. | |
| 1631 1 2 2.0 0.0 hess_fe = 0. | |
| 1632 1 3 3.0 0.0 hess_re = np.zeros((self.k_re2 + self.k_vc, self.k_re2 + self.k_vc)) | |
| 1633 1 2 2.0 0.0 hess_fere = np.zeros((self.k_re2 + self.k_vc, self.k_fe)) | |
| 1634 | |
| 1635 1 1 1.0 0.0 fac = self.n_totobs | |
| 1636 1 2 2.0 0.0 if self.reml: | |
| 1637 1 3 3.0 0.0 fac -= self.exog.shape[1] | |
| 1638 | |
| 1639 1 1 1.0 0.0 rvir = 0. | |
| 1640 1 1 1.0 0.0 xtvix = 0. | |
| 1641 1 3 3.0 0.0 xtax = [0.,] * (self.k_re2 + self.k_vc) | |
| 1642 1 2 2.0 0.0 m = self.k_re2 + self.k_vc | |
| 1643 1 4 4.0 0.0 B = np.zeros(m) | |
| 1644 1 2 2.0 0.0 D = np.zeros((m, m)) | |
| 1645 3 8 2.7 0.0 F = [[0.] * m for k in range(m)] | |
| 1646 61 134 2.2 0.0 for k, group in enumerate(self.group_labels): | |
| 1647 | |
| 1648 60 6386 106.4 0.0 cov_aug, cov_aug_inv, _ = self._augment_cov_re(cov_re, cov_re_inv, vcomp, group) | |
| 1649 | |
| 1650 60 128 2.1 0.0 exog = self.exog_li[k] | |
| 1651 60 124 2.1 0.0 ex_r, ex2_r = self._aex_r[k], self._aex_r2[k] | |
| 1652 60 43165 719.4 0.1 solver = _smw_solver(1., ex_r, ex2_r, cov_aug, cov_aug_inv) | |
| 1653 | |
| 1654 # The residuals | |
| 1655 60 150 2.5 0.0 resid = self.endog_li[k] | |
| 1656 60 98 1.6 0.0 if self.k_fe > 0: | |
| 1657 60 1819 30.3 0.0 expval = np.dot(exog, fe_params) | |
| 1658 60 478 8.0 0.0 resid = resid - expval | |
| 1659 | |
| 1660 60 33386 556.4 0.1 viexog = solver(exog) | |
| 1661 60 2019 33.6 0.0 xtvix += np.dot(exog.T, viexog) | |
| 1662 60 9100 151.7 0.0 vir = solver(resid) | |
| 1663 60 403 6.7 0.0 rvir += np.dot(resid, vir) | |
| 1664 | |
| 1665 180 557280 3096.0 1.4 for jj1, matl1, matr1, vsl1, vsr1, sym1 in self._gen_dV_dPar(ex_r, solver, group): | |
| 1666 | |
| 1667 120 17454 145.4 0.0 ul = np.dot(viexog.T, matl1) | |
| 1668 120 4104 34.2 0.0 ur = np.dot(matr1.T, vir) | |
| 1669 120 1149 9.6 0.0 hess_fere[jj1, :] += np.dot(ul, ur) | |
| 1670 120 196 1.6 0.0 if not sym1: | |
| 1671 ul = np.dot(viexog.T, matr1) | |
| 1672 ur = np.dot(matl1.T, vir) | |
| 1673 hess_fere[jj1, :] += np.dot(ul, ur) | |
| 1674 | |
| 1675 120 185 1.5 0.0 if self.reml: | |
| 1676 120 17082 142.3 0.0 ul = np.dot(viexog.T, matl1) | |
| 1677 120 221 1.8 0.0 ur = ul if sym1 else np.dot(viexog.T, matr1) | |
| 1678 120 466 3.9 0.0 ulr = np.dot(ul, ur.T) | |
| 1679 120 446 3.7 0.0 xtax[jj1] += ulr | |
| 1680 120 174 1.4 0.0 if not sym1: | |
| 1681 xtax[jj1] += ulr.T | |
| 1682 | |
| 1683 120 3857 32.1 0.0 ul = np.dot(vir, matl1) | |
| 1684 120 212 1.8 0.0 ur = ul if sym1 else np.dot(vir, matr1) | |
| 1685 120 876 7.3 0.0 B[jj1] += np.dot(ul, ur) * (1 if sym1 else 2) | |
| 1686 | |
| 1687 # V^{-1} * dV/d_theta | |
| 1688 120 6276293 52302.4 16.1 E = np.dot(vsl1, matr1.T) | |
| 1689 120 404 3.4 0.0 if not sym1: | |
| 1690 E += np.dot(vsr1, matl1.T) | |
| 1691 | |
| 1692 300 968825 3229.4 2.5 for jj2, matl2, matr2, vsl2, vsr2, sym2 in self._gen_dV_dPar(ex_r, solver, group, jj1): | |
| 1693 | |
| 1694 180 28761228 159784.6 73.8 re = np.dot(matr2.T, E) | |
| 1695 180 11733 65.2 0.0 rev = np.dot(re, vir[:, None]) | |
| 1696 180 11308 62.8 0.0 vl = np.dot(vir[None, :], matl2) | |
| 1697 180 3808 21.2 0.0 vt = 2*np.dot(vl, rev) | |
| 1698 | |
| 1699 180 312 1.7 0.0 if not sym2: | |
| 1700 le = np.dot(matl2.T, E) | |
| 1701 lev = np.dot(le, vir[:, None]) | |
| 1702 vr = np.dot(vir[None, :], matr2) | |
| 1703 vt += 2*np.dot(vr, lev) | |
| 1704 | |
| 1705 180 2127 11.8 0.0 D[jj1, jj2] += vt | |
| 1706 180 339 1.9 0.0 if jj1 != jj2: | |
| 1707 60 224 3.7 0.0 D[jj2, jj1] += vt | |
| 1708 | |
| 1709 180 24104 133.9 0.1 rt = np.sum(vsl2 * re.T) / 2 # trace dot | |
| 1710 180 301 1.7 0.0 if not sym2: | |
| 1711 rt += np.sum(vsr2 * le.T) / 2 # trace dot | |
| 1712 | |
| 1713 180 494 2.7 0.0 hess_re[jj1, jj2] += rt | |
| 1714 180 248 1.4 0.0 if jj1 != jj2: | |
| 1715 60 106 1.8 0.0 hess_re[jj2, jj1] += rt | |
| 1716 | |
| 1717 180 367 2.0 0.0 if self.reml: | |
| 1718 180 2148143 11934.1 5.5 ev = np.dot(E, viexog) | |
| 1719 180 33557 186.4 0.1 u1 = np.dot(viexog.T, matl2) | |
| 1720 180 25086 139.4 0.1 u2 = np.dot(matr2.T, ev) | |
| 1721 180 666 3.7 0.0 um = np.dot(u1, u2) | |
| 1722 180 4853 27.0 0.0 F[jj1][jj2] += um + um.T | |
| 1723 180 351 1.9 0.0 if not sym2: | |
| 1724 u1 = np.dot(viexog.T, matr2) | |
| 1725 u2 = np.dot(matl2.T, ev) | |
| 1726 um = np.dot(u1, u2) | |
| 1727 F[jj1][jj2] += um + um.T | |
| 1728 | |
| 1729 1 26 26.0 0.0 hess_fe -= fac * xtvix / rvir | |
| 1730 1 47 47.0 0.0 hess_re = hess_re - 0.5 * fac * (D/rvir - np.outer(B, B) / rvir**2) | |
| 1731 1 9 9.0 0.0 hess_fere = -fac * hess_fere / rvir | |
| 1732 | |
| 1733 1 1 1.0 0.0 if self.reml: | |
| 1734 3 184 61.3 0.0 QL = [np.linalg.solve(xtvix, x) for x in xtax] | |
| 1735 3 9 3.0 0.0 for j1 in range(self.k_re2 + self.k_vc): | |
| 1736 5 8 1.6 0.0 for j2 in range(j1 + 1): | |
| 1737 3 55 18.3 0.0 a = np.sum(QL[j1].T * QL[j2]) # trace dot | |
| 1738 3 162 54.0 0.0 a -= np.trace(np.linalg.solve(xtvix, F[j1][j2])) | |
| 1739 3 6 2.0 0.0 a *= 0.5 | |
| 1740 3 7 2.3 0.0 hess_re[j1, j2] += a | |
| 1741 3 6 2.0 0.0 if j1 > j2: | |
| 1742 1 2 2.0 0.0 hess_re[j2, j1] += a | |
| 1743 | |
| 1744 # Put the blocks together to get the Hessian. | |
| 1745 1 2 2.0 0.0 m = self.k_fe + self.k_re2 + self.k_vc | |
| 1746 1 4 4.0 0.0 hess = np.zeros((m, m)) | |
| 1747 1 4 4.0 0.0 hess[0:self.k_fe, 0:self.k_fe] = hess_fe | |
| 1748 1 3 3.0 0.0 hess[0:self.k_fe, self.k_fe:] = hess_fere.T | |
| 1749 1 3 3.0 0.0 hess[self.k_fe:, 0:self.k_fe] = hess_fere | |
| 1750 1 2 2.0 0.0 hess[self.k_fe:, self.k_fe:] = hess_re | |
| 1751 | |
| 1752 1 2 2.0 0.0 return hess | |
| Total time: 61.0446 s | |
| File: statsmodels/regression/mixed_linear_model.orig.py | |
| Function: fit at line 1807 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 1807 @profile | |
| 1808 def fit(self, start_params=None, reml=True, | |
| 1809 fe_pen=None, cov_pen=None, free=None, | |
| 1810 full_output=False, method='bfgs', **kwargs): | |
| 1811 """ | |
| 1812 Fit a linear mixed model to the data. | |
| 1813 | |
| 1814 Parameters | |
| 1815 ---------- | |
| 1816 start_params: array-like or MixedLMParams | |
| 1817 Starting values for the profile log-likeihood. If not a | |
| 1818 `MixedLMParams` instance, this should be an array | |
| 1819 containing the packed parameters for the profile | |
| 1820 log-likelihood, including the fixed effects | |
| 1821 parameters. | |
| 1822 reml : bool | |
| 1823 If true, fit according to the REML likelihood, else | |
| 1824 fit the standard likelihood using ML. | |
| 1825 cov_pen : CovariancePenalty object | |
| 1826 A penalty for the random effects covariance matrix | |
| 1827 fe_pen : Penalty object | |
| 1828 A penalty on the fixed effects | |
| 1829 free : MixedLMParams object | |
| 1830 If not `None`, this is a mask that allows parameters to be | |
| 1831 held fixed at specified values. A 1 indicates that the | |
| 1832 correspondinig parameter is estimated, a 0 indicates that | |
| 1833 it is fixed at its starting value. Setting the `cov_re` | |
| 1834 component to the identity matrix fits a model with | |
| 1835 independent random effects. Note that some optimization | |
| 1836 methods do not respect this contraint (bfgs and lbfgs both | |
| 1837 work). | |
| 1838 full_output : bool | |
| 1839 If true, attach iteration history to results | |
| 1840 method : string | |
| 1841 Optimization method. | |
| 1842 | |
| 1843 Returns | |
| 1844 ------- | |
| 1845 A MixedLMResults instance. | |
| 1846 """ | |
| 1847 | |
| 1848 1 3 3.0 0.0 _allowed_kwargs = ['gtol', 'maxiter'] | |
| 1849 1 3 3.0 0.0 for x in kwargs.keys(): | |
| 1850 if x not in _allowed_kwargs: | |
| 1851 raise ValueError("Argument %s not allowed for MixedLM.fit" % x) | |
| 1852 | |
| 1853 1 2 2.0 0.0 if method.lower() in ["newton", "ncg"]: | |
| 1854 raise ValueError("method %s not available for MixedLM" % method) | |
| 1855 | |
| 1856 1 2 2.0 0.0 self.reml = reml | |
| 1857 1 1 1.0 0.0 self.cov_pen = cov_pen | |
| 1858 1 2 2.0 0.0 self.fe_pen = fe_pen | |
| 1859 | |
| 1860 1 1 1.0 0.0 self._freepat = free | |
| 1861 | |
| 1862 1 2 2.0 0.0 if full_output: | |
| 1863 hist = [] | |
| 1864 else: | |
| 1865 1 2 2.0 0.0 hist = None | |
| 1866 | |
| 1867 1 2 2.0 0.0 if start_params is None: | |
| 1868 1 43 43.0 0.0 params = MixedLMParams(self.k_fe, self.k_re, self.k_vc) | |
| 1869 1 4 4.0 0.0 params.fe_params = np.zeros(self.k_fe) | |
| 1870 1 8 8.0 0.0 params.cov_re = np.eye(self.k_re) | |
| 1871 1 14 14.0 0.0 params.vcomp = np.ones(self.k_vc) | |
| 1872 else: | |
| 1873 if isinstance(start_params, MixedLMParams): | |
| 1874 params = start_params | |
| 1875 else: | |
| 1876 params = MixedLMParams.from_packed(start_params, self.k_fe, | |
| 1877 self.k_re, self.use_sqrt, | |
| 1878 with_fe=True) | |
| 1879 | |
| 1880 1 2 2.0 0.0 kwargs["retall"] = hist is not None | |
| 1881 1 1 1.0 0.0 if "disp" not in kwargs: | |
| 1882 1 2 2.0 0.0 kwargs["disp"] = False | |
| 1883 1 14 14.0 0.0 packed = params.get_packed(use_sqrt=self.use_sqrt, with_fe=False) | |
| 1884 1 3 3.0 0.0 rslt = super(MixedLM, self).fit(start_params=packed, | |
| 1885 1 1 1.0 0.0 skip_hessian=True, | |
| 1886 1 1 1.0 0.0 method=method, | |
| 1887 1 21922116 21922116.0 35.9 **kwargs) | |
| 1888 | |
| 1889 # The optimization succeeded | |
| 1890 1 22 22.0 0.0 params = np.atleast_1d(rslt.params) | |
| 1891 1 2 2.0 0.0 if hist is not None: | |
| 1892 hist.append(rslt.mle_retvals) | |
| 1893 | |
| 1894 1 1 1.0 0.0 converged = rslt.mle_retvals['converged'] | |
| 1895 1 2 2.0 0.0 if not converged: | |
| 1896 msg = "Gradient optimization failed." | |
| 1897 warnings.warn(msg, ConvergenceWarning) | |
| 1898 | |
| 1899 # Convert to the final parameterization (i.e. undo the square | |
| 1900 # root transform of the covariance matrix, and the profiling | |
| 1901 # over the error variance). | |
| 1902 1 3 3.0 0.0 params = MixedLMParams.from_packed(params, self.k_fe, self.k_re, | |
| 1903 1 66 66.0 0.0 use_sqrt=self.use_sqrt, with_fe=False) | |
| 1904 1 2 2.0 0.0 cov_re_unscaled = params.cov_re | |
| 1905 1 1 1.0 0.0 vcomp_unscaled = params.vcomp | |
| 1906 1 80065 80065.0 0.1 fe_params = self.get_fe_params(cov_re_unscaled, vcomp_unscaled) | |
| 1907 1 3 3.0 0.0 params.fe_params = fe_params | |
| 1908 1 47347 47347.0 0.1 scale = self.get_scale(fe_params, cov_re_unscaled, vcomp_unscaled) | |
| 1909 1 6 6.0 0.0 cov_re = scale * cov_re_unscaled | |
| 1910 1 3 3.0 0.0 vcomp = scale * vcomp_unscaled | |
| 1911 | |
| 1912 1 1 1.0 0.0 if (((self.k_re > 0) and (np.min(np.abs(np.diag(cov_re))) < 0.01)) or | |
| 1913 1 20 20.0 0.0 ((self.k_vc > 0) and (np.min(np.abs(vcomp)) < 0.01))): | |
| 1914 msg = "The MLE may be on the boundary of the parameter space." | |
| 1915 warnings.warn(msg, ConvergenceWarning) | |
| 1916 | |
| 1917 # Compute the Hessian at the MLE. Note that this is the | |
| 1918 # Hessian with respect to the random effects covariance matrix | |
| 1919 # (not its square root). It is used for obtaining standard | |
| 1920 # errors, not for optimization. | |
| 1921 1 38990445 38990445.0 63.9 hess = self.hessian(params) | |
| 1922 1 15 15.0 0.0 hess_diag = np.diag(hess) | |
| 1923 1 2 2.0 0.0 if free is not None: | |
| 1924 pcov = np.zeros_like(hess) | |
| 1925 pat = self._freepat.get_packed(use_sqrt=False, with_fe=True) | |
| 1926 ii = np.flatnonzero(pat) | |
| 1927 hess_diag = hess_diag[ii] | |
| 1928 if len(ii) > 0: | |
| 1929 hess1 = hess[np.ix_(ii, ii)] | |
| 1930 pcov[np.ix_(ii, ii)] = np.linalg.inv(-hess1) | |
| 1931 else: | |
| 1932 1 52 52.0 0.0 pcov = np.linalg.inv(-hess) | |
| 1933 1 28 28.0 0.0 if np.any(hess_diag >= 0): | |
| 1934 msg = "The Hessian matrix at the estimated parameter values is not positive definite." | |
| 1935 warnings.warn(msg, ConvergenceWarning) | |
| 1936 | |
| 1937 # Prepare a results class instance | |
| 1938 1 14 14.0 0.0 params_packed = params.get_packed(use_sqrt=False, with_fe=True) | |
| 1939 1 4250 4250.0 0.0 results = MixedLMResults(self, params_packed, pcov / scale) | |
| 1940 1 2 2.0 0.0 results.params_object = params | |
| 1941 1 2 2.0 0.0 results.fe_params = fe_params | |
| 1942 1 2 2.0 0.0 results.cov_re = cov_re | |
| 1943 1 1 1.0 0.0 results.vcomp = vcomp | |
| 1944 1 2 2.0 0.0 results.scale = scale | |
| 1945 1 1 1.0 0.0 results.cov_re_unscaled = cov_re_unscaled | |
| 1946 1 3 3.0 0.0 results.method = "REML" if self.reml else "ML" | |
| 1947 1 1 1.0 0.0 results.converged = converged | |
| 1948 1 2 2.0 0.0 results.hist = hist | |
| 1949 1 1 1.0 0.0 results.reml = self.reml | |
| 1950 1 2 2.0 0.0 results.cov_pen = self.cov_pen | |
| 1951 1 2 2.0 0.0 results.k_fe = self.k_fe | |
| 1952 1 2 2.0 0.0 results.k_re = self.k_re | |
| 1953 1 2 2.0 0.0 results.k_re2 = self.k_re2 | |
| 1954 1 3 3.0 0.0 results.k_vc = self.k_vc | |
| 1955 1 2 2.0 0.0 results.use_sqrt = self.use_sqrt | |
| 1956 1 3 3.0 0.0 results.freepat = self._freepat | |
| 1957 | |
| 1958 1 7 7.0 0.0 return MixedLMResultsWrapper(results) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment