home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d3xx / d386 / xlispstat.lha / XLispStat / src2.lzh / XLisp-Stat / kernel.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-10-05  |  2.2 KB  |  101 lines
  1. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  2. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  3. /* You may give out copies of this software; for conditions see the    */
  4. /* file COPYING included with this distribution.                       */
  5.  
  6. #include "xlisp.h"
  7. #include "osdef.h"
  8. #ifdef ANSI
  9. #include "xlsproto.h"
  10. #else
  11. #include "xlsfun.h"
  12. #endif ANSI
  13.  
  14. #ifdef ANSI
  15. double kernel(double,double,double,int);
  16. #else
  17. double kernel();
  18. #endif ANSI
  19.  
  20. #ifndef ROOT2PI
  21. #define ROOT2PI 2.50662827463100050241
  22. #endif  ROOT2PI
  23.  
  24. #ifndef nil
  25. #define nil 0L
  26. #endif  nil
  27.  
  28. static double kernel(x, y, w, type)
  29.      double x, y, w;
  30.      int type;
  31. {
  32.   double z, k;
  33.   
  34.   if (w > 0.0) {
  35.     z = (x - y) / w;
  36.     switch (type) {
  37.     case 'B':
  38.       w = 2.0 * w;
  39.       z = 0.5 * z;
  40.       if (-0.5 < z && z < 0.5) {
  41.         z = (1.0 - 4 * z * z);
  42.         k = 15.0 * z * z / 8.0;
  43.       }
  44.       else k = 0.0;
  45.       break;
  46.     case 'G':
  47.       w = 0.25 * w;
  48.       z = 4.0 * z;
  49.       k = exp(- 0.5 * z * z) / ROOT2PI; 
  50.       break;
  51.     case 'U':
  52.       w = 1.5 * w;
  53.       z = .75 * z;
  54.       k = (fabs(z) < 0.5) ? 1.0 : 0.0;
  55.       break;
  56.     case 'T': 
  57.       if (-1.0 <= z && z < 0.0) k = 1.0 + z;
  58.       else if (0.0 <= z && z < 1.0) k = 1.0 - z;
  59.       else k = 0.0;
  60.       break;
  61.     default:  k = 0.0; break;
  62.     }
  63.     k = k / w;
  64.   }
  65.   else k = 0.0;
  66.  
  67.   return(k);
  68. }
  69.  
  70. int kernel_smooth(x, y, n, width, wts, wds, xs, ys, ns, ktype)
  71.      double /* *x, *y,*/ width/*,  *wts, *wds, *xs, *ys */;/* changed JKL */
  72.      RVector x, y, wts, wds, xs, ys;
  73.      int n, ns, ktype;
  74. {
  75.   int i, j;
  76.   double wsum, ysum, lwidth, lwt, xmin, xmax;
  77.  
  78.   if (n < 1) return(1);
  79.   if (width <= 0.0) {
  80.     if (n < 2) return(1);
  81.     for (xmin = xmax = x[0], i = 1; i < n; i++) {
  82.       if (xmin > x[i]) xmin = x[i];
  83.       if (xmax < x[i]) xmax = x[i];
  84.     }
  85.     width = (xmax - xmin) / (1 + log((double) n));
  86.   }
  87.  
  88.   for (i = 0; i < ns; i++) {
  89.     for (j = 0, wsum = 0.0, ysum = 0.0; j < n; j++) {
  90.       lwidth = (wds != nil) ? width * wds[j] : width;
  91.       lwt = kernel(xs[i], x[j], lwidth, ktype);
  92.       if (wts != nil) lwt *= wts[j];
  93.       wsum += lwt;
  94.       if (y != nil) ysum += lwt * y[j];
  95.     }
  96.     if (y != nil) ys[i] = (wsum > 0.0) ? ysum / wsum : 0.0;
  97.     else ys[i] = wsum / n;
  98.   }
  99.   return(0);
  100. }
  101.