sascha@2645: package de.intevation.flys.artifacts.math;
sascha@2645: 
sascha@2646: import org.apache.commons.math.MathException;
sascha@2646: 
sascha@2645: import org.apache.commons.math.stat.descriptive.moment.Mean;
sascha@2645: import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
sascha@2645: 
sascha@2645: import org.apache.commons.math.distribution.TDistributionImpl;
sascha@2645: 
sascha@3011: import java.util.Collections;
sascha@2645: import java.util.List;
sascha@2645: import java.util.ArrayList;
sascha@2645: 
sascha@2646: import org.apache.log4j.Logger;
sascha@2646: 
sascha@2645: public class Outlier
sascha@2645: {
sascha@3011:     public static final double DEFAULT_ALPHA = 0.05;
sascha@3011: 
sascha@2646:     private static Logger log = Logger.getLogger(Outlier.class);
sascha@2646: 
sascha@3011:     public static class IndexedValue
sascha@3011:     implements          Comparable<IndexedValue>
sascha@3011:     {
sascha@2645:         protected int    index;
sascha@2645:         protected double value;
sascha@2645: 
sascha@2645:         public IndexedValue() {
sascha@2645:         }
sascha@2645: 
sascha@2645:         public IndexedValue(int index, double value) {
sascha@2645:             this.index = index;
sascha@2645:             this.value = value;
sascha@2645:         }
sascha@2645: 
sascha@2645:         public int getIndex() {
sascha@2645:             return index;
sascha@2645:         }
sascha@2645: 
sascha@2645:         public void setIndex(int index) {
sascha@2645:             this.index = index;
sascha@2645:         }
sascha@2645: 
sascha@2645:         public double getValue() {
sascha@2645:             return value;
sascha@2645:         }
sascha@2645: 
sascha@2645:         public void setValue(double value) {
sascha@2645:             this.value = value;
sascha@2645:         }
sascha@3011: 
sascha@3011:         @Override
sascha@3011:         public int compareTo(IndexedValue other) {
sascha@3011:             int diff = index - other.index;
sascha@3011:             if (index < 0) return -1;
sascha@3011:             return index > 0 ? +1 : 0;
sascha@3011:         }
sascha@2645:     } // class IndexedValue
sascha@2645: 
sascha@3011:     public static class Outliers {
sascha@3011: 
sascha@3011:         protected List<IndexedValue> retained;
sascha@3011:         protected List<IndexedValue> removed;
sascha@3011: 
sascha@3011:         public Outliers() {
sascha@3011:         }
sascha@3011: 
sascha@3011:         public Outliers(
sascha@3011:             List<IndexedValue> retained,
sascha@3011:             List<IndexedValue> removed
sascha@3011:         ) {
sascha@3011:             this.retained = retained;
sascha@3011:             this.removed  = removed;
sascha@3011:         }
sascha@3011: 
sascha@3011:         public boolean hasOutliers() {
sascha@3011:             return !removed.isEmpty();
sascha@3011:         }
sascha@3011: 
sascha@3011:         public List<IndexedValue> getRetained() {
sascha@3011:             return retained;
sascha@3011:         }
sascha@3011: 
sascha@3011:         public void setRetained(List<IndexedValue> retained) {
sascha@3011:             this.retained = retained;
sascha@3011:         }
sascha@3011: 
sascha@3011:         public List<IndexedValue> getRemoved() {
sascha@3011:             return removed;
sascha@3011:         }
sascha@3011: 
sascha@3011:         public void setRemoved(List<IndexedValue> removed) {
sascha@3011:             this.removed = removed;
sascha@3011:         }
sascha@3011:     } // class Outliers
sascha@3011: 
sascha@2645:     public Outlier() {
sascha@2645:     }
sascha@2645: 
sascha@3011:     public static Outliers findOutliers(List<IndexedValue> inputValues) {
sascha@3011:         return findOutliers(inputValues, DEFAULT_ALPHA);
sascha@3011:     }
sascha@3011: 
sascha@3011:     public static Outliers findOutliers(
sascha@2645:         List<IndexedValue> inputValues,
sascha@2645:         double             alpha
sascha@2645:     ) {
sascha@2645:         ArrayList<IndexedValue> outliers = new ArrayList<IndexedValue>();
sascha@2645: 
sascha@2645:         ArrayList<IndexedValue> values =
sascha@2645:             new ArrayList<IndexedValue>(inputValues);
sascha@2645: 
sascha@2645:         for (;;) {
sascha@2645:             int N = values.size();
sascha@2645: 
sascha@2645:             if (N < 4) {
sascha@2645:                 break;
sascha@2645:             }
sascha@2645: 
sascha@2645:             Mean mean = new Mean();
sascha@2645:             StandardDeviation std = new StandardDeviation();
sascha@2645: 
sascha@2645:             for (IndexedValue value: values) {
sascha@2645:                 mean.increment(value.getValue());
sascha@3011:                 std .increment(value.getValue());
sascha@2645:             }
sascha@2645: 
sascha@2645:             double m = mean.getResult();
sascha@2645:             double s = std.getResult();
sascha@2645: 
sascha@2645:             double maxZ = -Double.MAX_VALUE;
sascha@2645:             int iv = -1;
sascha@2646:             for (int i = N-1; i >= 0; --i) {
sascha@2645:                 IndexedValue v = values.get(i);
sascha@2645:                 double z = Math.abs(m - v.getValue())/s;
sascha@2645:                 if (z > maxZ) {
sascha@2645:                     maxZ = z;
sascha@2645:                     iv = i;
sascha@2645:                 }
sascha@2645:             }
sascha@2645: 
sascha@2646:             double t = Math.sqrt((N*(N-2)*maxZ*maxZ)
sascha@2646:                 /((N-1)*(N-1) - N*maxZ*maxZ));
sascha@2645: 
sascha@2645:             TDistributionImpl tdist = new TDistributionImpl(N-2);
sascha@2645: 
sascha@2646:             try {
sascha@2646:                 double p = tdist.cumulativeProbability(t);
sascha@2646: 
sascha@2646:                 if (p < alpha) {
sascha@2646:                     outliers.add(values.get(iv));
sascha@2646:                     values.remove(iv);
sascha@2646:                 }
sascha@3011:                 else {
sascha@3011:                     break;
sascha@3011:                 }
sascha@2646:             }
sascha@2646:             catch (MathException me) {
sascha@2646:                 log.error(me);
sascha@2646:             }
sascha@2645:         }
sascha@2645: 
sascha@3011:         Collections.sort(outliers);
sascha@2645: 
sascha@3011:         return new Outliers(values, outliers);
sascha@2645:     }
sascha@2645: }
sascha@2645: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :