lisp-stat-model package (still need to add to ASDF).
[CommonLispStat.git] / lib / derivatives.c
blobf7ff8c4b8ee1757190f59ae989274189a9cdbec5
1 /* derivatives - for Bayes code in XLISP-STAT and S */
2 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
3 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
4 /* You may give out copies of this software; for conditions see the */
5 /* file COPYING included with this distribution. */
7 /*
8 #ifdef __APPLE__
9 #include <stdlib.h>
10 #endif
13 #include <stdlib.h>
15 # include "xmath.h"
16 extern double macheps();
18 # define nil 0L
19 # define FALSE 0
20 # define TRUE 1
22 void
23 numergrad(size_t n, double *x, double *grad, double *fsum,
24 int ffun(double *, double *, double *, double **),
25 double h, double *typx)
26 /* int (*ffun)();*/
28 size_t i;
29 double old_xi, f1, f2, hi;
31 for (i = 0; i < n; i++) {
32 old_xi = x[i];
33 hi = (typx != nil) ? typx[i] * h : h;
34 x[i] = old_xi + hi;
35 (*ffun)(x, &f1, nil, nil);
36 x[i] = old_xi - hi;
37 (*ffun)(x, &f2, nil, nil);
38 x[i] = old_xi;
39 grad[i] = (f1 - f2) / (2.0 * hi);
40 fsum[i] = f1 + f2;
44 void
45 numerhess(size_t n, double *x, double **hess, double f, double *fsum,
46 int *ffun(double *, double *, double *, double *),
47 double h, double *typx)
49 size_t i, j;
50 double old_xi, old_xj, f1, f2, hi, hj;
52 for (i = 0; i < n; i++) {
53 hi = (typx != nil) ? typx[i] * h : h;
54 hess[i][i] = (fsum[i] - 2 * f) / (hi * hi);
55 for (j = i + 1; j < n; j++) {
56 hj = (typx != nil) ? typx[j] * h : h;
57 old_xi = x[i];
58 old_xj = x[j];
59 x[i] = old_xi + hi;
60 x[j] = old_xj + hj;
61 (*ffun)(x, &f1, nil, nil);
62 x[i] = old_xi - hi;
63 x[j] = old_xj - hj;
64 (*ffun)(x, &f2, nil, nil);
65 x[i] = old_xi;
66 x[j] = old_xj;
67 hess[i][j] = (2 * f + f1 + f2 - fsum[i] - fsum[j]) / (2.0 * hi * hj);
68 hess[j][i] = hess[i][j];