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 / statistics.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-10-09  |  6.2 KB  |  279 lines
  1. /* statistics - Basic statistical functions for XLISP-STAT.            */
  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 "xlisp.h"
  8. #include "osdef.h"
  9. #ifdef ANSI
  10. #include "xlproto.h"
  11. #include "xlsproto.h"
  12. #include "osproto.h"
  13. #else
  14. #include "xlfun.h"
  15. #include "xlsfun.h"
  16. #include "osfun.h"
  17. #endif ANSI
  18.  
  19. #ifdef ANSI
  20. LVAL datareduce1(LVAL (*)(),LVAL (*)(),LVAL,int),lastcdr(LVAL),
  21.      elementlist(LVAL),basequant(void),quant(void),base_ifelse(void);
  22. int all_simple(LVAL);
  23. void make_number(Number *,LVAL),base_mean(int *,Number *,LVAL);
  24. #else
  25. LVAL datareduce1),lastcdr(),
  26.      elementlist(),basequant(),quant(),base_ifelse();
  27. int all_simple();
  28. void make_number),base_mean();
  29. #endif ANSI
  30.  
  31. LVAL xssum()  { return(datareduce1(xssum, xsradd, cvfixnum((FIXTYPE) 0), FALSE)); }
  32. LVAL xsprod() { return(datareduce1(xsprod, xsrmul, cvfixnum((FIXTYPE) 1), FALSE)); }
  33. LVAL xsmin() { return(datareduce1(xsmin, xsrmin, NIL, FALSE)); }
  34. LVAL xsmax() { return(datareduce1(xsmax, xsrmax, NIL, FALSE)); }
  35. LVAL xscount()  { return(datareduce1(xscount, xsradd, cvfixnum((FIXTYPE) 0), TRUE)); }
  36.  
  37. static LVAL datareduce1(f, bf, nullval, count)
  38.     LVAL (*f)(), (*bf)(), nullval;
  39.     int count;
  40. {
  41.   LVAL fcn, x, result;
  42.  
  43.   switch (xlargc) {
  44.   case 0: result = nullval;
  45.   case 1: 
  46.     if (compoundp(peekarg(0))) {
  47.       xlstkcheck(2);
  48.       xlsave(x);
  49.       xlsave(fcn);
  50.       fcn = cvsubr(bf, SUBR, 0);
  51.       x = subr_map_elements(f);
  52.       x = compounddataseq(x);
  53.       result = reduce(fcn, x, FALSE, NIL);
  54.       xlpopn(2);
  55.     }
  56.     else result = (count) ? cvfixnum((FIXTYPE) 1) : xlgetarg();
  57.     break;
  58.   default:
  59.     xlsave1(x);
  60.     x = makearglist(xlargc, xlargv);
  61.     result = xscallsubr1(f, x);
  62.     xlpop();
  63.   }
  64.   return(result);
  65. }
  66.  
  67. /* check if all elements of the sequence x are simple (i. e. not compound) */
  68. LOCAL int all_simple(x)
  69.      LVAL x;
  70. {
  71.   int i, n;
  72.  
  73.   if (! sequencep(x)) xlerror("not a sequence", x);
  74.   if (listp(x)) {
  75.     for (; consp(x); x = cdr(x)) 
  76.       if (compoundp(car(x))) return(FALSE);
  77.   }
  78.   else {
  79.     n = getsize(x);
  80.     for (i = 0; i < n; i++) 
  81.       if (compoundp(getelement(x, i))) return(FALSE);
  82.   }
  83.   return(TRUE);
  84. }
  85.  
  86. static LVAL lastcdr(x)
  87.     LVAL x;
  88. {
  89.   LVAL last = NIL;
  90.   
  91.   for (; consp(x); x = cdr(x)) last = x;
  92.   
  93.   return(last);
  94. }
  95.  
  96. static LVAL elementlist(x)
  97.     LVAL x;
  98. {
  99.   LVAL next, last, result;
  100.   
  101.   if (!compoundp(x)) result = consa(x);
  102.   else {
  103.     xlprot1(x);
  104.     x = compounddataseq(x);
  105.     x = (vectorp(x)) ? coerce_to_list(x) : copylist(x);
  106.     if (all_simple(x)) result = x;
  107.     else {
  108.       for (next = x; consp(next); next = cdr(next))
  109.         rplaca(next, elementlist(car(next)));
  110.       result = car(x);
  111.       last = lastcdr(car(x));
  112.       for (next = cdr(x); consp(next); next = cdr(next)) {
  113.         rplacd(last, car(next));
  114.         last = lastcdr(car(next));
  115.       }
  116.     }
  117.     xlpop();
  118.   }
  119.   return(result);
  120. }
  121.       
  122. LVAL elementseq(x)
  123.     LVAL x;
  124. {
  125.   if (! compoundp(x)) xlerror("not a compound data item", x);
  126.   x = compounddataseq(x);
  127.   if (all_simple(x)) return(x);
  128.   else return(elementlist(x));
  129. }
  130.  
  131. LVAL xselement_seq() { return(elementseq(xlgetarg())); }
  132.  
  133. static LVAL data;
  134.  
  135. static LVAL basequant()
  136. {
  137.   LVAL prob;
  138.   double p, lowelem, highelem;
  139.   int lowk, highk, n;
  140.   
  141.   prob = xlgetarg();
  142.   p = makedouble(prob);
  143.   if ((0.0 > p) || (1.0 < p))
  144.     xlerror("not a probability", prob);
  145.   xllastarg(); 
  146.   
  147.   n = getsize(data);
  148.   lowk = (n - 1) * p;
  149.   highk = (n - 1) - floor((n - 1) * (1 - p));
  150.   lowelem = makedouble(getelement(data, lowk));
  151.   highelem = makedouble(getelement(data, highk));
  152.   return(cvflonum((FLOTYPE) (lowelem + highelem) / 2.0));
  153. }
  154.  
  155. static LVAL quant() { return(recursive_subr_map_elements(basequant, quant)); }
  156.  
  157. LVAL xsquantile()
  158. {
  159.   LVAL result;
  160.   
  161.   xlsave1(data);
  162.   data = xlgetarg();
  163.   data = xscallsubr1(xssortdata, data);
  164.   data = coerce_to_vector(data);
  165.   result = quant();
  166.   xlpop();
  167.   return(result);
  168. }
  169.  
  170. static LVAL base_ifelse()
  171. {
  172.   LVAL a, b, c;
  173.   
  174.   a = xlgetarg();
  175.   b = xlgetarg();
  176.   c = xlgetarg();
  177.   xllastarg();
  178.   
  179.   return((a != NIL) ? b : c);
  180. }
  181.  
  182. LVAL xsifelse() { return(subr_map_elements(base_ifelse)); }
  183. /* in xlsdef.h JKL
  184. typedef struct {
  185.   double real, imag;
  186.   int complex;
  187. } Number;
  188. */
  189. static void make_number(num, x)
  190.      Number *num;
  191.      LVAL x;
  192. {
  193.   if (fixp(x) || floatp(x)) {
  194.     num->real = makedouble(x);
  195.     num->imag = 0.0;
  196.     num->complex = FALSE;
  197.   }
  198.   else if (complexp(x)) {
  199.     num->real = makedouble(realpart(x));
  200.     num->imag = makedouble(imagpart(x));
  201.     num->complex = TRUE;
  202.   }
  203.   else xlerror("not a number", x);
  204. }
  205.  
  206. static void base_mean(count, mean, x)
  207.      int *count;
  208.      Number *mean;
  209.      LVAL x;
  210. {
  211.   LVAL y;
  212.   Number num;
  213.   double c, p, q;
  214.   int i, n;
  215.  
  216.   if (! compoundp(x)) {
  217.     c = *count; p = c / (c + 1.0); q = 1.0 - p;
  218.     make_number(&num, x);
  219.     mean->real = p * mean->real + q * num.real;
  220.     mean->complex = mean->complex || num.complex;
  221.     if (mean->complex) mean->imag = p * mean->imag + q * num.imag;
  222.     (*count)++;
  223.   }
  224.   else {
  225.     x = compounddataseq(x);
  226.     n = seqlen(x);
  227.     for (i = 0; i < n; i++) {
  228.       y = getnextelement(&x, i);
  229.       base_mean(count, mean, y);
  230.     }
  231.   }
  232. }
  233.  
  234. LVAL xsmean()
  235. {
  236.   Number mean;
  237.   int count;
  238.   LVAL x;
  239.  
  240.   x = xlgetarg();
  241.   xllastarg();
  242.  
  243.   mean.real = 0.0; mean.imag = 0.0; mean.complex = FALSE;
  244.   count = 0;
  245.   base_mean(&count, &mean, x);
  246.   if (mean.complex) return(newdcomplex(mean.real,mean.imag));
  247.   else return(cvflonum((FLOTYPE) mean.real));
  248. }
  249.  
  250. LVAL xssample()
  251. {
  252.   LVAL x, result, temp;
  253.   int n, N, replace, i, j;
  254.   
  255.   x = xsgetsequence();
  256.   n = getfixnum(xlgafixnum());
  257.   N = seqlen(x);
  258.   replace = (moreargs()) ? (xlgetarg() != NIL) : FALSE;
  259.   
  260.   xlstkcheck(2);
  261.   xlprotect(x);
  262.   xlsave(result);
  263.   x = (listp(x)) ? coerce_to_vector(x) : copyvector(x);
  264.   /* result = NIL; not needed JKL */
  265.   if (N > 0 && n > 0) {
  266.     for (i = 0; i < n; i++) {
  267.       j = (replace) ? osrand(N) : i + osrand(N - i);
  268.       result = cons(getelement(x, j), result);
  269.       if (! replace) {           /* swap elements i and j */
  270.         temp = getelement(x, i);
  271.         setelement(x, i, getelement(x, j));
  272.         setelement(x, j, temp);
  273.       }
  274.     }
  275.   }
  276.   xlpopn(2);
  277.   return(result);
  278. }
  279.