lib/derivatives.c

   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.                       */
   6
   7 # include "xmath.h"
   8 extern double macheps();
   9
  10 # define nil 0L
  11 # define FALSE 0
  12 # define TRUE 1
  13
  14 void
  15 numergrad(int n, double *x, double *grad, double *fsum,
  16           int ffun(double *, double *, double *, double **),
  17           double h, double *typx)
  18      /*     int (*ffun)();*/
  19 {
  20   int i;
  21   double old_xi, f1, f2, hi;
  22
  23   for (i = 0; i < n; i++) {
  24     old_xi = x[i];
  25     hi = (typx != nil) ? typx[i] * h : h;
  26     x[i] = old_xi + hi;
  27     (*ffun)(x, &f1, nil, nil);
  28     x[i] = old_xi - hi;
  29     (*ffun)(x, &f2, nil, nil);
  30     x[i] = old_xi;
  31     grad[i] = (f1 - f2) / (2.0 * hi);
  32     fsum[i] = f1 + f2;
  33   }
  34 }
  35
  36 void
  37 numerhess(int n, double *x, double **hess, double f, double *fsum,
  38           int *ffun(double *, double *, double *, double *),
  39           double h, double *typx)
  40 {
  41   int i, j;
  42   double old_xi, old_xj, f1, f2, hi, hj;
  43
  44   for (i = 0; i < n; i++) {
  45     hi = (typx != nil) ? typx[i] * h : h;
  46     hess[i][i] = (fsum[i] - 2 * f) / (hi * hi);
  47     for (j = i + 1; j < n; j++) {
  48       hj = (typx != nil) ? typx[j] * h : h;
  49       old_xi = x[i];
  50       old_xj = x[j];
  51       x[i] = old_xi + hi;
  52       x[j] = old_xj + hj;
  53       (*ffun)(x, &f1, nil, nil);
  54       x[i] = old_xi - hi;
  55       x[j] = old_xj - hj;
  56       (*ffun)(x, &f2, nil, nil);
  57       x[i] = old_xi;
  58       x[j] = old_xj;
  59       hess[i][j] = (2 * f + f1 + f2 - fsum[i] - fsum[j]) / (2.0 * hi * hj);
  60       hess[j][i] = hess[i][j];
  61     }
  62   }
  63 }