correlation.h

/***************************************************************************
                          correlation.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2000-2021 by Thies Jochimsen
    email                : thies@jochimsen.de
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef CORRELATION_H
#define CORRELATION_H
 
#include<odindata/statistics.h>
 
struct correlationResult {
 
  correlationResult() : r(0.0), p(0.0), z(0.0) {}
 
  double r;
 
  double p;
 
  double z;
};
 
 
 
 
template <typename T, int N_rank>
correlationResult correlation(const Array<T,N_rank>& x, const Array<T,N_rank>& y) {
  Log<OdinData> odinlog("","correlation");
 
  correlationResult result;
 
  if(x.shape()!=y.shape()) {
    ODINLOG(odinlog,errorLog) << "Shape mismatch" << STD_endl;
    return result;
  }
 
  int n=x.numElements();
 
  if(n<2) return result;
 
  statisticResult xstat=statistics(x);
  statisticResult ystat=statistics(y);
 
  result.r=secureDivision(
    sum( (x-xstat.mean)*(y-ystat.mean) ) ,
    sqrt( sum( (x-xstat.mean) * (x-xstat.mean) ) ) * sqrt ( sum( (y-ystat.mean) * (y-ystat.mean) ) )
  );
 
#ifdef HAVE_ERFC
  result.p=0.5*erfc(fabs(result.r)*sqrt(double(n)/2.0)); // Eq 14.5.2 in NRC
#else
#error erfc() is missing
#endif
 
  double argument=secureDivision(1.0+result.r,1.0-result.r);
  if(argument) {
//    result.z=0.5*log(argument);  // Eq 14.5.6 in NRC
    if(n>3) result.z=0.5*sqrt(double(n-3))*log(argument);   // http://www.fmrib.ox.ac.uk/~stuart/thesis/chapter_6/section6_4.html
  }
 
  return result;
}
 
 
 
template <typename T, int N_rank>
correlationResult kendall(const Array<T,N_rank>& x, const Array<T,N_rank>& y) {
  Log<OdinData> odinlog("","kendall");
 
  correlationResult result;
 
  if(x.shape()!=y.shape()) {
    ODINLOG(odinlog,errorLog) << "Shape mismatch" << STD_endl;
    return result;
  }
 
  int n=x.size();
 
  int nx=0;
  int ny=0;
  int diff=0;
  for(int j=0; j<n; j++) {
    for(int i=j+1; i<n; i++) {
      TinyVector<int,N_rank> ii=index2extent(x.shape(),i);
      TinyVector<int,N_rank> jj=index2extent(x.shape(),j);
      float dx=x(jj)-x(ii);
      float dy=y(jj)-y(ii);
      if(dx) nx++;
      if(dy) ny++;
      float dd=dx*dy;
      if(dd>0.0) diff++;
      if(dd<0.0) diff--;
    }
  }
 
  double tau=secureDivision(diff,sqrt(double(nx))*sqrt(double(ny)));
 
  result.r=tau;
 
  result.z=3.0*tau*sqrt( secureDivision(n*(n-1),2*(2*n+5) ) );  // Eq. 30.4 in Handbook of Parametric and Nonparametric Statistical Procedures
 
  return result;
}
 
 
 
struct fmriResult {
 
  fmriResult() : Sbaseline(0.0), Srest(0.0), Sstim(0.0), rel_diff(0.0), rel_err(0.0) {}
 
  float Sbaseline;
 
  float Srest;
 
  float Sstim;
  
  float rel_diff;
 
  float rel_err;
};
 
 
fmriResult fmri_eval(const Data<float,1>& timecourse, const Data<float,1>& designvec);
 
 
#endif