source: branches/DataAnalysis/HeuristicLab.Problems.DataAnalysis.MultiVariate.TimeSeriesPrognosis/3.3/Evaluators/OnlineTheilsUStatisticEvaluator.cs @ 6627

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Problems.DataAnalysis.Evaluators;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Parameters;

namespace HeuristicLab.Problems.DataAnalysis.MultiVariate.TimeSeriesPrognosis {
  public class OnlineTheilsUStatisticEvaluator : IOnlineTimeSeriesPrognosisEvaluator {
    private OnlineMeanAndVarianceCalculator squaredErrorMeanCalculator;
    private OnlineMeanAndVarianceCalculator unbiasedEstimatorMeanCalculator;
    private double prevOriginal;
    private int windowSize;
    private Queue<double> movingAverageWindow;

    public double TheilsUStatistic {
      get {
        return Math.Sqrt(squaredErrorMeanCalculator.Mean) / Math.Sqrt(unbiasedEstimatorMeanCalculator.Mean);
      }
    }

    public OnlineTheilsUStatisticEvaluator()
      : this(1) {
    }

    public OnlineTheilsUStatisticEvaluator(int movingAverageWindowSize) {
      this.windowSize = movingAverageWindowSize;
      movingAverageWindow = new Queue<double>(windowSize);
      squaredErrorMeanCalculator = new OnlineMeanAndVarianceCalculator();
      unbiasedEstimatorMeanCalculator = new OnlineMeanAndVarianceCalculator();
      Reset();
    }

    #region IOnlineEvaluator Members
    public double Value {
      get { return TheilsUStatistic; }
    }

    public void Add(double original, double estimated) {
      if (double.IsInfinity(original) || double.IsNaN(original) || double.IsInfinity(estimated) || double.IsNaN(estimated)) {
        throw new ArgumentException("Theil's U-statistic is not defined for series containing NaN or infinity values.");
      }

      if (!double.IsNaN(prevOriginal)) {
        // error of predicted change
        double errorEstimatedChange = (estimated - original);
        squaredErrorMeanCalculator.Add(errorEstimatedChange * errorEstimatedChange);

        // calculate trend observed in the MA window
        double d = CalculateTrend(movingAverageWindow);
        
        // shift window forward
        if (movingAverageWindow.Count == windowSize) {
          movingAverageWindow.Dequeue();
        }
        movingAverageWindow.Enqueue(original);

        double errorNoChange = (original - prevOriginal * (1+d));
        unbiasedEstimatorMeanCalculator.Add(errorNoChange * errorNoChange);
      }
    }

    private double CalculateTrend(Queue<double> movingAverageWindow) {
      double[] xs = movingAverageWindow.ToArray();
      double sum = 0.0;
      for (int i = 0; i < xs.Length - 1; i++) {
        sum += (xs[i + 1] - xs[i]) / xs[i];
      }
      return sum / xs.Length;
    }

    public void Reset() {
      prevOriginal = double.NaN;
      squaredErrorMeanCalculator.Reset();
      unbiasedEstimatorMeanCalculator.Reset();
      movingAverageWindow.Clear();
    }

    #endregion

    #region IOnlineTimeSeriesPrognosisEvaluator Members

    public void StartNewPredictionWindow(double referenceOriginalValue) {
      prevOriginal = referenceOriginalValue;
      movingAverageWindow.Enqueue(referenceOriginalValue);
    }

    #endregion
  }
}