MultivariateStatistics.cpp

Go to the documentation of this file.

#include <float.h>
#include <string>
#include <iostream>
#include "MultivariateStatistics.h"
#include "IException.h"


namespace Isis {
  MultivariateStatistics::MultivariateStatistics() {
    Reset();
  }

  void MultivariateStatistics::Reset() {
    p_x.Reset();
    p_y.Reset();
    p_sumxy = 0.0;

    p_validPixels = 0;
    p_invalidPixels = 0;
    p_totalPixels = 0;
  }

  MultivariateStatistics::~MultivariateStatistics() {};
  void MultivariateStatistics::AddData(const double *x, const double *y,
                                       const unsigned int count) {
    for(unsigned int i = 0; i < count; i++) {
      double yVal = y[i];
      double xVal = x[i];
      p_totalPixels++;

      if(Isis::IsValidPixel(xVal) && Isis::IsValidPixel(yVal)) {
        p_x.AddData(xVal);
        p_y.AddData(yVal);
        p_sumxy += xVal * yVal;
        p_validPixels++;
      }
      else {
        p_invalidPixels++;
      }
    }
  }

  void MultivariateStatistics::AddData(double x, double y, unsigned int count) {
    p_totalPixels += count;

    if(IsValidPixel(x) && IsValidPixel(y)) {
      p_sumxy += x * y * count;
      p_validPixels += count;

      for (unsigned int i = 0; i < count; i++) {
        p_x.AddData(x);
        p_y.AddData(y);
      }
    }
    else {
      p_invalidPixels += count;
    }
  }

  void MultivariateStatistics::RemoveData(const double *x, const double *y,
                                          const unsigned int count) {
    for(unsigned int i = 0; i < count; i++) {
      p_totalPixels--;

      if(Isis::IsValidPixel(x[i]) && Isis::IsValidPixel(y[i])) {
        p_x.RemoveData(&x[i], 1);
        p_y.RemoveData(&y[i], 1);
        p_sumxy -= x[i] * y[i];
        p_validPixels--;
      }
      else {
        p_invalidPixels--;
      }
    }

    if(p_totalPixels < 0) {
      std::string m = "You are removing non-existant data in [MultivariateStatistics::RemoveData]";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
  }

  double MultivariateStatistics::Covariance() const {
    if(p_validPixels <= 1) return Isis::NULL8;
    double covar = p_sumxy - p_y.Average() * p_x.Sum() - p_x.Average() * p_y.Sum() +
                   p_x.Average() * p_y.Average() * (double)p_validPixels;
    return covar / (double)(p_validPixels - 1);
  }

  double MultivariateStatistics::Correlation() const {

    if(p_validPixels <= 1) return Isis::NULL8;
    double covar = Covariance();
    double stdX = p_x.StandardDeviation();
    double stdY = p_y.StandardDeviation();
    if(stdX == 0.0 || stdX == Isis::NULL8) return Isis::NULL8;
    if(stdY == 0.0 || stdY == Isis::NULL8) return Isis::NULL8;
    if(covar == Isis::NULL8) return Isis::NULL8;
    return covar / (stdX * stdY);
  }

  BigInt MultivariateStatistics::TotalPixels() const {
    return p_totalPixels;
  }

  BigInt MultivariateStatistics::ValidPixels() const {
    return p_validPixels;
  }

  BigInt MultivariateStatistics::InvalidPixels() const {
    return p_invalidPixels;
  }

  void MultivariateStatistics::LinearRegression(double &a, double &b) const {
    // From Modern Elementary Statistics - 5th edition, Freund, pp 367
    double denom = (double)p_validPixels * p_x.SumSquare() - p_x.Sum() * p_x.Sum();
    if(denom == 0.0) {
      std::string msg = "Unable to compute linear regression in Multivariate Statistics";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    a = p_y.Sum() * p_x.SumSquare() - p_x.Sum() * p_sumxy;
    a = a / denom;

    b = (double)p_validPixels * p_sumxy - p_x.Sum() * p_y.Sum();
    b = b / denom;
  }

  double MultivariateStatistics::SumXY() const {
    return p_sumxy;
  }

  Isis::Statistics MultivariateStatistics::X() const {
    return p_x;
  };

  Isis::Statistics MultivariateStatistics::Y() const {
    return p_y;
  };
}