SCM

SCM Repository

[matrix] Diff of /pkg/src/ssclme.c
ViewVC logotype

Diff of /pkg/src/ssclme.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 176, Fri May 28 13:23:44 2004 UTC revision 237, Thu Jul 1 03:20:33 2004 UTC
# Line 1310  Line 1310 
1310  }  }
1311    
1312  /**  /**
1313     * Return the Hessian of the ML or REML deviance.  This is a
1314     * placeholder until I work out the evaluation of the analytic
1315     * Hessian, which probably will involve several helper functions.
1316     *
1317     * @param x pointer to an ssclme object
1318     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1319     * @param Uncp pointer to a logical scalar indicating if the
1320     * unconstrained parameterization is to be used
1321     *
1322     * @return pointer to an approximate Hessian matrix
1323     */
1324    SEXP ssclme_Hessian(SEXP x, SEXP REMLp, SEXP Uncp)
1325    {
1326        int j, ncoef = coef_length(length(GET_SLOT(x, Matrix_OmegaSym)),
1327                                   INTEGER(GET_SLOT(x, Matrix_ncSym))),
1328            unc = asLogical(Uncp);
1329        SEXP ans = PROTECT(allocMatrix(REALSXP, ncoef, ncoef)),
1330            base = PROTECT(unc ? ssclme_coefUnc(x) : ssclme_coef(x)),
1331            current = PROTECT(duplicate(base)),
1332            gradient;
1333    
1334        for (j = 0; j < ncoef; j++) {
1335            double delta = (REAL(base)[j] ? 1.e-7 * REAL(base)[j] : 1.e-7);
1336            int i;
1337    
1338            for (i = 0; i < ncoef; i++) REAL(current)[i] = REAL(base)[i];
1339            REAL(current)[j] += delta/2.;
1340            if (unc) {
1341                ssclme_coefGetsUnc(x, current);
1342            } else {
1343                ssclme_coefGets(x, current);
1344            }
1345            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1346            for (i = 0; i < ncoef; i++) REAL(ans)[j * ncoef + i] = REAL(gradient)[i];
1347            UNPROTECT(1);
1348            REAL(current)[j] -= delta;
1349            if (unc) {
1350                ssclme_coefGetsUnc(x, current);
1351            } else {
1352                ssclme_coefGets(x, current);
1353            }
1354            PROTECT(gradient = ssclme_gradient(x, REMLp, Uncp));
1355            for (i = 0; i < ncoef; i++)
1356                REAL(ans)[j * ncoef + i] = (REAL(ans)[j * ncoef + i] - REAL(gradient)[i])/
1357                    delta;
1358            UNPROTECT(1);
1359        }
1360        UNPROTECT(3);
1361        return ans;
1362    }
1363    
1364    /**
1365   * Calculate and return the fitted values.   * Calculate and return the fitted values.
1366   *   *
1367   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object

Legend:
Removed from v.176  
changed lines
  Added in v.237

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge