/*-----------------------------------------------------------------------------

Copyright (C) 2017

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-------------------------------------------------------------------------------

Function      simple - computes simple statistics for each column of a realmat

Syntax        #include "libscl.h"
              realmat simple(const realmat& X); 
              realmat simple(const trealmat& X); 

Prototype in  libscl.h

Description   For each column of X simple statistics are computed

Remark        If X is a trealmat, then simple statistics are for each 
              column of T(X).

Return value  If S = simple(X), x_j is the jth column of X, and s_j is
              the jth column of S, then s_j contains
	      1. mean of x_j
	      2. standard deviation of x_j
	      3. variance of x_j
	      4. skewness of x_j
	      5. kurtosis of x_j
	      
Functions     Library: (none)
called        libscl:  realmat

------------------------------------------------------------------------------*/

#include "libscl.h"

namespace scl {

  realmat simple(const realmat& X) 
  {
    INTEGER n = X.nrow();
    INTEGER k = X.ncol();

    if (n <= 1) {
      error("Error, simple, not enough data; need two or more observations");
    }

    realmat S(5,k);

    for (INTEGER j = 1; j<=k; ++j) {

      REAL ave = 0.0;

      for (INTEGER i=1; i<=n; ++i) {
        ave += X(i,j);  
      }
      
      ave /= REAL(n);

      REAL adev, sdev, var, skew, curt;
      
      adev = sdev = var = skew = curt = 0.0;

      REAL r, p;

      for (INTEGER i=1; i<=n; ++i) {
        adev += fabs(r = X(i,j) - ave);
	var  += (p = r*r);
	skew += (p *= r);
	curt += (p *= r);
      }
      
      REAL fn = REAL(n);

      adev /= fn;
      var  /= (fn-1.0);
      sdev = sqrt(var);
      if (var) {
        skew /= (fn * var * sdev);
        curt =  curt/(fn * var * var) - 3.0;
      }
      else {
        skew = curt = REAL_MAX;
      }
    
      S(1,j) = ave;
      S(2,j) = sdev;
      S(3,j) = var;
      S(4,j) = skew;
      S(5,j) = curt;
    }

    return S;
  }

  realmat simple(const trealmat& X) 
  {
    INTEGER n = X.nrow();
    INTEGER k = X.ncol();

    if (n <= 1) {
      error("Error, simple, not enough data; need two or more observations");
    }

    realmat S(5,k);

    for (INTEGER j = 1; j<=k; ++j) {

      REAL ave = 0.0;

      for (INTEGER i=1; i<=n; ++i) {
        ave += X(i,j);  
      }
      
      ave /= REAL(n);

      REAL adev, sdev, var, skew, curt;
      
      adev = sdev = var = skew = curt = 0.0;

      REAL r, p;

      for (INTEGER i=1; i<=n; ++i) {
        adev += fabs(r = X(i,j) - ave);
	var  += (p = r*r);
	skew += (p *= r);
	curt += (p *= r);
      }
      
      REAL fn = REAL(n);

      adev /= fn;
      var  /= (fn-1.0);
      sdev = sqrt(var);
      if (var) {
        skew /= (fn * var * sdev);
        curt =  curt/(fn * var * var) - 3.0;
      }
      else {
        skew = curt = REAL_MAX;
      }
    
      S(1,j) = ave;
      S(2,j) = sdev;
      S(3,j) = var;
      S(4,j) = skew;
      S(5,j) = curt;
    }

    return S;
  }

}
