correlation.h

/***************************************************************************
                          correlation.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2001 by Thies Jochimsen
    email                : jochimse@cns.mpg.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

---