home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / magazine / drdobbs / 1992 / 06 / noisy.asc < prev    next >
Text File  |  1992-05-18  |  4KB  |  151 lines
  1. _FINDING SIGNIFICANCE IN NOISY DATA_
  2. by Roy Kimbrell
  3.  
  4. [LISTING ONE]
  5.  
  6. #include <stdio.h>
  7. #include <stdlib.h>
  8. #include <math.h>
  9. #include <console.h>
  10.  
  11. void I_Lfilter(void);
  12. unsigned Lfilter(double Y, double confidence);
  13.  
  14. #define MAIN
  15.  
  16. #define T           15  /* number of historical values */
  17. #define S           30  /* number of stages */
  18. #define SIGNIFICANT     1
  19. #define NOT_SIGNIFICANT !SIGNIFICANT
  20.  
  21. /* The following macros define the circular buffer. There are T+2 values in 
  22. the buffer to allow exactly T historical values. The "last" and "plast" values
  23. are additional to the T values and are not used in the computations--they 
  24. exist in the buffer to allow them to be subtracted from current values.
  25. */
  26.  
  27. #define prior           (now == 0 ? T+1 : now-1)
  28. #define last            (now == T ? 0 : (now == T+1 ? 1 : now+2))
  29. #define plast           (now == T+1 ? 0 : now+1)
  30. #define cycle           (now == T+1 ? now = 0 : now++)
  31.  
  32. static int now; /* index of the current historical value */
  33.  
  34. static double Ef;           /* final error sum */
  35. static double priorEf, priorEff;    /* prior error sums */
  36. static double ef[S+1][T+2];     /* array of forward errors */
  37. static double eb[S+1][T+2];     /* array of backward errors */
  38. double Y_hat;               /* expected Y */
  39. static double Efb[S+1], Eff[S+1];
  40.  
  41. void I_Lfilter(void){
  42.     
  43.     int s, t;
  44.     now = 0;
  45.     Ef = priorEf = priorEff = Y_hat = 0.0;
  46.     for(s=0; s<=S; s++){
  47.         Efb[s] = Eff[s] = 0;
  48.         for(t=0; t<=T+1; t++) ef[s][t] = eb[s][t] = 0.0;
  49.         }
  50.     } /* I_Lfilter */
  51.     
  52. unsigned Lfilter(double Y, double confidence){
  53.  
  54.     double k, z=0.0, sd=0.0;
  55.     int s;
  56.     static unsigned iteration;
  57.     
  58.     ef[0][now] = eb[0][now] = Y;
  59.      for (s=1; s<=S; s++){
  60.       Efb[s] += ef[s-1][now] * eb[s-1][prior] - ef[s-1][last] * eb[s-1][plast];
  61.          Eff[s] += ef[s-1][now] * ef[s-1][now] - ef[s-1][last] * ef[s-1][last];
  62.         k = Efb[s] / Eff[s];
  63.         ef[s][now] = ef[s-1][now] - k * eb[s-1][prior];
  64.         eb[s][now] = eb[s-1][prior] - k * ef[s-1][now];
  65.         }
  66.     Y_hat = Y - ef[S][prior];
  67.     Ef += ef[S][now] - ef[S][last];
  68.     if (iteration != 0){
  69.         sd = sqrt((T * priorEff - priorEf * priorEf)/(T*(T-1)));
  70.         z = (ef[S][prior] - Ef / T) / sd;
  71.         }
  72.     priorEff = Eff[S];  /* use the output of the last stage */
  73.     priorEf = Ef;
  74.     iteration++;
  75.     cycle;
  76.     if (fabs(z) > confidence)
  77.         return SIGNIFICANT;
  78.     else
  79.         return NOT_SIGNIFICANT;
  80.     } /* Lfilter */
  81.  
  82. #ifdef MAIN
  83.  
  84. #define CONFIDENCE      2.0
  85.  
  86. void main(int ac, char **av){
  87.     
  88.     char buf[80];
  89.     unsigned i=0, significant;
  90.     double Y;
  91.     extern double Y_hat;
  92.     
  93.     while (fgets(buf,sizeof(buf),stdin)){
  94.         Y = atof(buf);
  95.         significant = Lfilter(Y,CONFIDENCE);
  96.         printf("%d\t%.4f\t%.4f",i++,Y,Y_hat);
  97.         if (significant)
  98.             printf("\tSIGNIFICANT\n");
  99.         else
  100.             printf("\n");
  101.         }
  102.     } /* main */
  103.  
  104. #endif
  105.  
  106.  
  107. Example 1: The value of to reflects trends in the relative 
  108. differences between forward and (prior) backward error values.
  109.  
  110.  
  111. At stage s,
  112.    for (t = now; t != last; t = prior(t)){
  113.     Efbs = eft s-1 * eb prior(t) s-1;
  114.     Effs = (eft s-1)2;
  115.     }
  116.    ks =  F(Efbs/Effs );
  117.  
  118.  
  119.  
  120. Example 2: Calculating the filtered value.
  121.  
  122. ef0 now = ef0 now = Y;
  123. for (s = 1; s # S; s++){
  124.     Efbs += efs-1 now * ebs-1 prior - efs-1 last * ebs-1 plast;
  125.     Effs += efs-1 now  * efs-1 now    - efs-1 last  * efs-1 last ;
  126.     k =  F(Efbs ,Effs);
  127.     efs now = efs-1 now - k * ebs-1 prior;
  128.     ebs now = ebs-1 prior - k * efs-1 now;
  129. /*      Comment The following happens to be one way of computing 
  130.         o(Y,^) but we actually do it a different way:
  131.      o(Y,^) += k * ebs-1 prior;
  132. */
  133.     }
  134.  o(Y,^) = Y - efS prior;
  135. Ef += efS now - efS last;
  136. if (iteration != 0){
  137.     sd =  R(, F(T * priorEff - (priorEf)2,T * (T-1)));
  138.     z =  F( B bc [(e S up3(f) S do3(S prior) -  F(Ef,T)),sd);
  139.     }
  140. priorEff = EffS;
  141. priorEf = Ef;
  142. iteration++;
  143. cycle;
  144. if ( B bc |(z) > confidence)
  145.     return SIGNIFICANT;
  146. else
  147.     return NOT_SIGNIFICANT;
  148.  
  149.  
  150.  
  151.