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