Created
June 22, 2015 18:35
-
-
Save saketkc/23fd682d4ece6b57c895 to your computer and use it in GitHub Desktop.
overwrites
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
| Total time: 3.02405 s | |
| File: statsmodels/regression/mixed_linear_model.mod.py | |
| Function: _smw_solver at line 356 | |
| 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 49134 8.4 1.6 qmat = BI + AtA / s | |
| 382 5880 2868186 487.8 94.8 qmati = scipy.linalg.solve(qmat, A.T, overwrite_a=True, check_finite=False) #np.linalg.solve(qmat, A.T) | |
| 383 | |
| 384 5880 102470 17.4 3.4 @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 4264 0.7 0.1 return solver | |
| Total time: 24.0576 s | |
| File: statsmodels/regression/mixed_linear_model.mod.py | |
| Function: solver at line 384 | |
| Line # Hits Time Per Hit % Time Line Contents | |
| ============================================================== | |
| 384 @profile | |
| 385 def solver(rhs): | |
| 386 13440 13665521 1016.8 56.8 ql = np.dot(qmati, rhs) | |
| 387 13440 9524465 708.7 39.6 ql = np.dot(A, ql) | |
| 388 13440 859182 63.9 3.6 rslt = rhs / s - ql / s**2 | |
| 389 13440 8410 0.6 0.0 return rslt | |
| Total time: 0.183025 s | |
| File: statsmodels/regression/mixed_linear_model.mod.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 2077 1.4 1.1 p = A.shape[0] | |
| 421 1500 8634 5.8 4.7 ld = p * np.log(s) | |
| 422 1500 11393 7.6 6.2 qmat = BI + AtA / s | |
| 423 1500 159054 106.0 86.9 _, ld1 = np.linalg.slogdet(qmat) | |
| 424 1500 1867 1.2 1.0 return B_logdet + ld + ld1 | |
| Total time: 19.2829 s | |
| File: statsmodels/regression/mixed_linear_model.mod.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 9635596 6176.7 50.0 axr = solver(ex_r) | |
| 1333 | |
| 1334 # Regular random effects | |
| 1335 1560 928 0.6 0.0 jj = 0 | |
| 1336 1560 4133 2.6 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 3931 0.9 0.0 for ky in self._vc_names: | |
| 1347 3120 4291 1.4 0.0 if group in self.exog_vc[ky]: | |
| 1348 3120 1743 0.6 0.0 if max_ix is not None and jj > max_ix: | |
| 1349 60 68 1.1 0.0 return | |
| 1350 3060 2189 0.7 0.0 mat = self.exog_vc[ky][group] | |
| 1351 3060 9625683 3145.6 49.9 axmat = solver(mat) | |
| 1352 3060 2294 0.7 0.0 yield jj, mat, mat, axmat, axmat, True | |
| 1353 3060 1995 0.7 0.0 jj += 1 | |
| Total time: 22.3189 s | |
| File: statsmodels/regression/mixed_linear_model.mod.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 28 1.2 0.0 if type(params) is not MixedLMParams: | |
| 1369 23 35 1.5 0.0 params = MixedLMParams.from_packed(params, self.k_fe, | |
| 1370 23 17 0.7 0.0 self.k_re, self.use_sqrt, | |
| 1371 23 1521 66.1 0.0 with_fe=False) | |
| 1372 | |
| 1373 23 18 0.8 0.0 if profile_fe: | |
| 1374 23 1978680 86029.6 8.9 params.fe_params = self.get_fe_params(params.cov_re, params.vcomp) | |
| 1375 | |
| 1376 23 26 1.1 0.0 if self.use_sqrt: | |
| 1377 23 20338446 884280.3 91.1 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 20 0.9 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 11 0.5 0.0 if profile_fe: | |
| 1387 23 52 2.3 0.0 return np.concatenate((score_re, score_vc)) | |
| 1388 else: | |
| 1389 return np.concatenate((score_fe, score_re, score_vc)) | |
| Total time: 20.2441 s | |
| File: statsmodels/regression/mixed_linear_model.mod.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 28 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 439 19.1 0.0 cov_re_inv = np.linalg.inv(cov_re) | |
| 1438 except np.linalg.LinAlgError: | |
| 1439 cov_re_inv = None | |
| 1440 | |
| 1441 23 92 4.0 0.0 score_fe = np.zeros(self.k_fe) | |
| 1442 23 65 2.8 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 27 1.2 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 22 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 22 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 43 1.9 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 48 2.1 0.0 rvavr = np.zeros(self.k_re2 + self.k_vc) | |
| 1473 | |
| 1474 1403 2122 1.5 0.0 for group_ix, group in enumerate(self.group_labels): | |
| 1475 | |
| 1476 1380 93150 67.5 0.5 cov_aug, cov_aug_inv, _ = self._augment_cov_re(cov_re, cov_re_inv, vcomp, group) | |
| 1477 | |
| 1478 1380 2394 1.7 0.0 exog = self.exog_li[group_ix] | |
| 1479 1380 2210 1.6 0.0 ex_r, ex2_r = self._aex_r[group_ix], self._aex_r2[group_ix] | |
| 1480 1380 821419 595.2 4.1 solver = _smw_solver(1., ex_r, ex2_r, cov_aug, cov_aug_inv) | |
| 1481 | |
| 1482 # The residuals | |
| 1483 1380 2724 2.0 0.0 resid = self.endog_li[group_ix] | |
| 1484 1380 1613 1.2 0.0 if self.k_fe > 0: | |
| 1485 1380 46352 33.6 0.2 expval = np.dot(exog, fe_params) | |
| 1486 1380 9717 7.0 0.0 resid = resid - expval | |
| 1487 | |
| 1488 1380 1715 1.2 0.0 if self.reml: | |
| 1489 1380 944717 684.6 4.7 viexog = solver(exog) | |
| 1490 1380 45748 33.2 0.2 xtvix += np.dot(exog.T, viexog) | |
| 1491 | |
| 1492 # Contributions to the covariance parameter gradient | |
| 1493 1380 224998 163.0 1.1 vir = solver(resid) | |
| 1494 4140 17297152 4178.1 85.4 for jj, matl, matr, vsl, vsr, sym in self._gen_dV_dPar(ex_r, solver, group): | |
| 1495 2760 152018 55.1 0.8 dlv[jj] = np.sum(matr * vsl) # trace dot | |
| 1496 2760 3889 1.4 0.0 if not sym: | |
| 1497 dlv[jj] += np.sum(matl * vsr) # trace dot | |
| 1498 | |
| 1499 2760 103546 37.5 0.5 ul = np.dot(vir, matl) | |
| 1500 2760 6389 2.3 0.0 ur = ul.T if sym else np.dot(matr.T, vir) | |
| 1501 2760 8025 2.9 0.0 ulr = np.dot(ul, ur) | |
| 1502 2760 6913 2.5 0.0 rvavr[jj] += ulr | |
| 1503 2760 2861 1.0 0.0 if not sym: | |
| 1504 rvavr[jj] += ulr.T | |
| 1505 | |
| 1506 2760 3132 1.1 0.0 if self.reml: | |
| 1507 2760 409631 148.4 2.0 ul = np.dot(viexog.T, matl) | |
| 1508 2760 4988 1.8 0.0 ur = ul.T if sym else np.dot(matr.T, viexog) | |
| 1509 2760 7883 2.9 0.0 ulr = np.dot(ul, ur) | |
| 1510 2760 10056 3.6 0.0 xtax[jj] += ulr | |
| 1511 2760 3431 1.2 0.0 if not sym: | |
| 1512 xtax[jj] += ulr.T | |
| 1513 | |
| 1514 # Contribution of log|V| to the covariance parameter | |
| 1515 # gradient. | |
| 1516 1380 1781 1.3 0.0 if self.k_re > 0: | |
| 1517 score_re -= 0.5 * dlv[0:self.k_re2] | |
| 1518 1380 1567 1.1 0.0 if self.k_vc > 0: | |
| 1519 1380 8917 6.5 0.0 score_vc -= 0.5 * dlv[self.k_re2:] | |
| 1520 | |
| 1521 1380 7548 5.5 0.0 rvir += np.dot(resid, vir) | |
| 1522 | |
| 1523 1380 1532 1.1 0.0 if calc_fe: | |
| 1524 xtvir += np.dot(exog.T, vir) | |
| 1525 | |
| 1526 23 28 1.2 0.0 fac = self.n_totobs | |
| 1527 23 24 1.0 0.0 if self.reml: | |
| 1528 23 28 1.2 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 25 1.1 0.0 if self.k_vc > 0: | |
| 1536 23 255 11.1 0.0 score_vc += 0.5 * fac * rvavr[self.k_re2:] / rvir | |
| 1537 | |
| 1538 23 32 1.4 0.0 if self.reml: | |
| 1539 23 1547 67.3 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 103 1.5 0.0 for j in range(self.k_vc): | |
| 1543 46 782 17.0 0.0 score_vc[j] += 0.5 * np.sum(xtvixi.T * xtax[self.k_re2 + j]) # trace dot | |
| 1544 | |
| 1545 23 23 1.0 0.0 return score_fe, score_re, score_vc | |
| Total time: 20.3378 s | |
| File: statsmodels/regression/mixed_linear_model.mod.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 20335705 884161.1 100.0 score_fe, score_re, score_vc = self.score_full(params, calc_fe=calc_fe) | |
| 1581 23 339 14.7 0.0 params_vec = params.get_packed(use_sqrt=True, with_fe=True) | |
| 1582 | |
| 1583 23 55 2.4 0.0 score_full = np.concatenate((score_fe, score_re, score_vc)) | |
| 1584 23 14 0.6 0.0 scr = 0. | |
| 1585 161 126 0.8 0.0 for i in range(len(params_vec)): | |
| 1586 138 930 6.7 0.0 v = self._lin[i] + 2 * np.dot(self._quad[i], params_vec) | |
| 1587 138 549 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 27 1.2 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 16 0.7 0.0 return score_fe, score_re, score_vc | |
| Total time: 41.2382 s | |
| File: statsmodels/regression/mixed_linear_model.mod.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 1 1.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 2 2.0 0.0 xtax = [0.,] * (self.k_re2 + self.k_vc) | |
| 1642 1 1 1.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 142 2.3 0.0 for k, group in enumerate(self.group_labels): | |
| 1647 | |
| 1648 60 6479 108.0 0.0 cov_aug, cov_aug_inv, _ = self._augment_cov_re(cov_re, cov_re_inv, vcomp, group) | |
| 1649 | |
| 1650 60 130 2.2 0.0 exog = self.exog_li[k] | |
| 1651 60 121 2.0 0.0 ex_r, ex2_r = self._aex_r[k], self._aex_r2[k] | |
| 1652 60 38955 649.2 0.1 solver = _smw_solver(1., ex_r, ex2_r, cov_aug, cov_aug_inv) | |
| 1653 | |
| 1654 # The residuals | |
| 1655 60 156 2.6 0.0 resid = self.endog_li[k] | |
| 1656 60 98 1.6 0.0 if self.k_fe > 0: | |
| 1657 60 1888 31.5 0.0 expval = np.dot(exog, fe_params) | |
| 1658 60 487 8.1 0.0 resid = resid - expval | |
| 1659 | |
| 1660 60 42226 703.8 0.1 viexog = solver(exog) | |
| 1661 60 2010 33.5 0.0 xtvix += np.dot(exog.T, viexog) | |
| 1662 60 9703 161.7 0.0 vir = solver(resid) | |
| 1663 60 396 6.6 0.0 rvir += np.dot(resid, vir) | |
| 1664 | |
| 1665 180 748306 4157.3 1.8 for jj1, matl1, matr1, vsl1, vsr1, sym1 in self._gen_dV_dPar(ex_r, solver, group): | |
| 1666 | |
| 1667 120 17543 146.2 0.0 ul = np.dot(viexog.T, matl1) | |
| 1668 120 5630 46.9 0.0 ur = np.dot(matr1.T, vir) | |
| 1669 120 1160 9.7 0.0 hess_fere[jj1, :] += np.dot(ul, ur) | |
| 1670 120 184 1.5 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 189 1.6 0.0 if self.reml: | |
| 1676 120 17181 143.2 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 473 3.9 0.0 ulr = np.dot(ul, ur.T) | |
| 1679 120 462 3.9 0.0 xtax[jj1] += ulr | |
| 1680 120 166 1.4 0.0 if not sym1: | |
| 1681 xtax[jj1] += ulr.T | |
| 1682 | |
| 1683 120 3838 32.0 0.0 ul = np.dot(vir, matl1) | |
| 1684 120 203 1.7 0.0 ur = ul if sym1 else np.dot(vir, matr1) | |
| 1685 120 926 7.7 0.0 B[jj1] += np.dot(ul, ur) * (1 if sym1 else 2) | |
| 1686 | |
| 1687 # V^{-1} * dV/d_theta | |
| 1688 120 6278210 52318.4 15.2 E = np.dot(vsl1, matr1.T) | |
| 1689 120 425 3.5 0.0 if not sym1: | |
| 1690 E += np.dot(vsr1, matl1.T) | |
| 1691 | |
| 1692 300 1293123 4310.4 3.1 for jj2, matl2, matr2, vsl2, vsr2, sym2 in self._gen_dV_dPar(ex_r, solver, group, jj1): | |
| 1693 | |
| 1694 180 30494830 169415.7 73.9 re = np.dot(matr2.T, E) | |
| 1695 180 11760 65.3 0.0 rev = np.dot(re, vir[:, None]) | |
| 1696 180 11262 62.6 0.0 vl = np.dot(vir[None, :], matl2) | |
| 1697 180 3949 21.9 0.0 vt = 2*np.dot(vl, rev) | |
| 1698 | |
| 1699 180 318 1.8 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 2209 12.3 0.0 D[jj1, jj2] += vt | |
| 1706 180 337 1.9 0.0 if jj1 != jj2: | |
| 1707 60 232 3.9 0.0 D[jj2, jj1] += vt | |
| 1708 | |
| 1709 180 23846 132.5 0.1 rt = np.sum(vsl2 * re.T) / 2 # trace dot | |
| 1710 180 309 1.7 0.0 if not sym2: | |
| 1711 rt += np.sum(vsr2 * le.T) / 2 # trace dot | |
| 1712 | |
| 1713 180 481 2.7 0.0 hess_re[jj1, jj2] += rt | |
| 1714 180 251 1.4 0.0 if jj1 != jj2: | |
| 1715 60 96 1.6 0.0 hess_re[jj2, jj1] += rt | |
| 1716 | |
| 1717 180 353 2.0 0.0 if self.reml: | |
| 1718 180 2151158 11950.9 5.2 ev = np.dot(E, viexog) | |
| 1719 180 33738 187.4 0.1 u1 = np.dot(viexog.T, matl2) | |
| 1720 180 25453 141.4 0.1 u2 = np.dot(matr2.T, ev) | |
| 1721 180 668 3.7 0.0 um = np.dot(u1, u2) | |
| 1722 180 4964 27.6 0.0 F[jj1][jj2] += um + um.T | |
| 1723 180 340 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 27 27.0 0.0 hess_fe -= fac * xtvix / rvir | |
| 1730 1 45 45.0 0.0 hess_re = hess_re - 0.5 * fac * (D/rvir - np.outer(B, B) / rvir**2) | |
| 1731 1 8 8.0 0.0 hess_fere = -fac * hess_fere / rvir | |
| 1732 | |
| 1733 1 2 2.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 11 2.2 0.0 for j2 in range(j1 + 1): | |
| 1737 3 56 18.7 0.0 a = np.sum(QL[j1].T * QL[j2]) # trace dot | |
| 1738 3 168 56.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 5 1.7 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 3 3.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: 68.2131 s | |
| File: statsmodels/regression/mixed_linear_model.mod.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 2 2.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 3 3.0 0.0 if method.lower() in ["newton", "ncg"]: | |
| 1854 raise ValueError("method %s not available for MixedLM" % method) | |
| 1855 | |
| 1856 1 1 1.0 0.0 self.reml = reml | |
| 1857 1 2 2.0 0.0 self.cov_pen = cov_pen | |
| 1858 1 2 2.0 0.0 self.fe_pen = fe_pen | |
| 1859 | |
| 1860 1 2 2.0 0.0 self._freepat = free | |
| 1861 | |
| 1862 1 1 1.0 0.0 if full_output: | |
| 1863 hist = [] | |
| 1864 else: | |
| 1865 1 1 1.0 0.0 hist = None | |
| 1866 | |
| 1867 1 1 1.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 9 9.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 1 1.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 26827200 26827200.0 39.3 **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 2 2.0 0.0 converged = rslt.mle_retvals['converged'] | |
| 1895 1 1 1.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 2 2.0 0.0 params = MixedLMParams.from_packed(params, self.k_fe, self.k_re, | |
| 1903 1 65 65.0 0.0 use_sqrt=self.use_sqrt, with_fe=False) | |
| 1904 1 1 1.0 0.0 cov_re_unscaled = params.cov_re | |
| 1905 1 1 1.0 0.0 vcomp_unscaled = params.vcomp | |
| 1906 1 84099 84099.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 45318 45318.0 0.1 scale = self.get_scale(fe_params, cov_re_unscaled, vcomp_unscaled) | |
| 1909 1 7 7.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 19 19.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 41252083 41252083.0 60.5 hess = self.hessian(params) | |
| 1922 1 16 16.0 0.0 hess_diag = np.diag(hess) | |
| 1923 1 1 1.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 53 53.0 0.0 pcov = np.linalg.inv(-hess) | |
| 1933 1 29 29.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 4008 4008.0 0.0 results = MixedLMResults(self, params_packed, pcov / scale) | |
| 1940 1 3 3.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 2 2.0 0.0 results.cov_re_unscaled = cov_re_unscaled | |
| 1946 1 2 2.0 0.0 results.method = "REML" if self.reml else "ML" | |
| 1947 1 2 2.0 0.0 results.converged = converged | |
| 1948 1 1 1.0 0.0 results.hist = hist | |
| 1949 1 2 2.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 1 1.0 0.0 results.k_re2 = self.k_re2 | |
| 1954 1 4 4.0 0.0 results.k_vc = self.k_vc | |
| 1955 1 1 1.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