/*--------------------------------ScoreContributions ------------------------------*/ /* This is a modification of Num1Derivative that allows the function to be vector valued, as with Gauss's gradp. The function returns n x 1, the parameter is p x 1, and the derivative is n x p. The matrix of derivatives is passed as the third argument, just as with Num1Derivative. Michael Creel, mcreel@volcano.uab.es 21 Aug. 1998 */ const decl SQRT_EPS =1E-8; /* appr. square root of machine precision */ const decl DIFF_EPS1=5E-6; /* Rice's formula: log(DIFF_EPS)=log(MACH_EPS)/3 */ static dFiniteDiff1(const x) { return max( (fabs(x) + SQRT_EPS) * SQRT_EPS, DIFF_EPS1); } ScoreContributions(const func, vP, const avScore) { decl i, cp = rows(vP), left, right, fknowf = FALSE, p, h, f, fm, fp, v; for (i = 0; i < cp; i++) /* get 1st derivative by central difference */ { p = double(vP[i][0]); h = dFiniteDiff1(p); vP[i][0] = p + h; right = func(vP, &fp, 0, 0); if(i==0) v = new matrix[rows(fp)][cp]; vP[i][0] = p - h; left = func(vP, &fm, 0, 0); vP[i][0] = p; /* restore original parameter */ if (left && right) v[][i] = (fp - fm) / (2 * h); /* take central difference */ else return FALSE; } avScore[0] = v; return TRUE; }