_FINDING SIGNIFICANCE IN NOISY DATA_
by Roy Kimbrell

[LISTING ONE]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <console.h>

void I_Lfilter(void);
unsigned Lfilter(double Y, double confidence);

#define MAIN

#define T           15  /* number of historical values */
#define S           30  /* number of stages */
#define SIGNIFICANT     1
#define NOT_SIGNIFICANT !SIGNIFICANT

/* The following macros define the circular buffer. There are T+2 values in 
the buffer to allow exactly T historical values. The "last" and "plast" values
are additional to the T values and are not used in the computations--they 
exist in the buffer to allow them to be subtracted from current values.
*/

#define prior           (now == 0 ? T+1 : now-1)
#define last            (now == T ? 0 : (now == T+1 ? 1 : now+2))
#define plast           (now == T+1 ? 0 : now+1)
#define cycle           (now == T+1 ? now = 0 : now++)

static int now; /* index of the current historical value */

static double Ef;           /* final error sum */
static double priorEf, priorEff;    /* prior error sums */
static double ef[S+1][T+2];     /* array of forward errors */
static double eb[S+1][T+2];     /* array of backward errors */
double Y_hat;               /* expected Y */
static double Efb[S+1], Eff[S+1];

void I_Lfilter(void){
    
    int s, t;
    now = 0;
    Ef = priorEf = priorEff = Y_hat = 0.0;
    for(s=0; s<=S; s++){
        Efb[s] = Eff[s] = 0;
        for(t=0; t<=T+1; t++) ef[s][t] = eb[s][t] = 0.0;
        }
    } /* I_Lfilter */
    
unsigned Lfilter(double Y, double confidence){

    double k, z=0.0, sd=0.0;
    int s;
    static unsigned iteration;
    
    ef[0][now] = eb[0][now] = Y;
     for (s=1; s<=S; s++){
      Efb[s] += ef[s-1][now] * eb[s-1][prior] - ef[s-1][last] * eb[s-1][plast];
         Eff[s] += ef[s-1][now] * ef[s-1][now] - ef[s-1][last] * ef[s-1][last];
        k = Efb[s] / Eff[s];
        ef[s][now] = ef[s-1][now] - k * eb[s-1][prior];
        eb[s][now] = eb[s-1][prior] - k * ef[s-1][now];
        }
    Y_hat = Y - ef[S][prior];
    Ef += ef[S][now] - ef[S][last];
    if (iteration != 0){
        sd = sqrt((T * priorEff - priorEf * priorEf)/(T*(T-1)));
        z = (ef[S][prior] - Ef / T) / sd;
        }
    priorEff = Eff[S];  /* use the output of the last stage */
    priorEf = Ef;
    iteration++;
    cycle;
    if (fabs(z) > confidence)
        return SIGNIFICANT;
    else
        return NOT_SIGNIFICANT;
    } /* Lfilter */

#ifdef MAIN

#define CONFIDENCE      2.0

void main(int ac, char **av){
    
    char buf[80];
    unsigned i=0, significant;
    double Y;
    extern double Y_hat;
    
    while (fgets(buf,sizeof(buf),stdin)){
        Y = atof(buf);
        significant = Lfilter(Y,CONFIDENCE);
        printf("%d\t%.4f\t%.4f",i++,Y,Y_hat);
        if (significant)
            printf("\tSIGNIFICANT\n");
        else
            printf("\n");
        }
    } /* main */

#endif


Example 1: The value of to reflects trends in the relative 
differences between forward and (prior) backward error values.


At stage s,
   for (t = now; t != last; t = prior(t)){
	Efbs = eft s-1 * eb prior(t) s-1;
	Effs = (eft s-1)2;
	}
   ks =  F(Efbs/Effs );



Example 2: Calculating the filtered value.

ef0 now = ef0 now = Y;
for (s = 1; s # S; s++){
	Efbs += efs-1 now * ebs-1 prior - efs-1 last * ebs-1 plast;
	Effs += efs-1 now  * efs-1 now	- efs-1 last  * efs-1 last ;
	k =  F(Efbs ,Effs);
	efs now = efs-1 now - k * ebs-1 prior;
	ebs now = ebs-1 prior - k * efs-1 now;
/*      Comment The following happens to be one way of computing 
        o(Y,^) but we actually do it a different way:
	 o(Y,^) += k * ebs-1 prior;
*/
	}
 o(Y,^) = Y - efS prior;
Ef += efS now - efS last;
if (iteration != 0){
	sd =  R(, F(T * priorEff - (priorEf)2,T * (T-1)));
	z =  F( B bc [(e S up3(f) S do3(S prior) -  F(Ef,T)),sd);
	}
priorEff = EffS;
priorEf = Ef;
iteration++;
cycle;
if ( B bc |(z) > confidence)
	return SIGNIFICANT;
else
	return NOT_SIGNIFICANT;