correlation.h
1 /***************************************************************************
2  correlation.h - description
3  -------------------
4  begin : Fri Apr 6 2001
5  copyright : (C) 2000-2015 by Thies Jochimsen
6  email : thies@jochimsen.de
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #ifndef CORRELATION_H
19 #define CORRELATION_H
20 
21 #include<odindata/statistics.h>
22 
33 
34  correlationResult() : r(0.0), p(0.0), z(0.0) {}
35 
39  double r;
40 
44  double p;
45 
49  double z;
50 };
51 
52 
54 
55 
60 template <typename T, int N_rank>
61 correlationResult correlation(const Array<T,N_rank>& x, const Array<T,N_rank>& y) {
62  Log<OdinData> odinlog("","correlation");
63 
64  correlationResult result;
65 
66  if(x.shape()!=y.shape()) {
67  ODINLOG(odinlog,errorLog) << "Shape mismatch" << STD_endl;
68  return result;
69  }
70 
71  int n=x.numElements();
72 
73  if(n<2) return result;
74 
75  statisticResult xstat=statistics(x);
76  statisticResult ystat=statistics(y);
77 
78  result.r=secureDivision(
79  sum( (x-xstat.mean)*(y-ystat.mean) ) ,
80  sqrt( sum( (x-xstat.mean) * (x-xstat.mean) ) ) * sqrt ( sum( (y-ystat.mean) * (y-ystat.mean) ) )
81  );
82 
83 #ifdef HAVE_ERFC
84  result.p=0.5*erfc(fabs(result.r)*sqrt(double(n)/2.0)); // Eq 14.5.2 in NRC
85 #else
86 #error erfc() is missing
87 #endif
88 
89  double argument=secureDivision(1.0+result.r,1.0-result.r);
90  if(argument) {
91 // result.z=0.5*log(argument); // Eq 14.5.6 in NRC
92  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
93  }
94 
95  return result;
96 }
97 
99 
100 
105 template <typename T, int N_rank>
106 correlationResult kendall(const Array<T,N_rank>& x, const Array<T,N_rank>& y) {
107  Log<OdinData> odinlog("","kendall");
108 
109  correlationResult result;
110 
111  if(x.shape()!=y.shape()) {
112  ODINLOG(odinlog,errorLog) << "Shape mismatch" << STD_endl;
113  return result;
114  }
115 
116  int n=x.size();
117 
118  int nx=0;
119  int ny=0;
120  int diff=0;
121  for(int j=0; j<n; j++) {
122  for(int i=j+1; i<n; i++) {
123  TinyVector<int,N_rank> ii=index2extent(x.shape(),i);
124  TinyVector<int,N_rank> jj=index2extent(x.shape(),j);
125  float dx=x(jj)-x(ii);
126  float dy=y(jj)-y(ii);
127  if(dx) nx++;
128  if(dy) ny++;
129  float dd=dx*dy;
130  if(dd>0.0) diff++;
131  if(dd<0.0) diff--;
132  }
133  }
134 
135  double tau=secureDivision(diff,sqrt(double(nx))*sqrt(double(ny)));
136 
137  result.r=tau;
138 
139  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
140 
141  return result;
142 }
143 
145 
146 
151 struct fmriResult {
152 
153  fmriResult() : Sbaseline(0.0), Srest(0.0), Sstim(0.0), rel_diff(0.0), rel_err(0.0) {}
154 
158  float Sbaseline;
159 
163  float Srest;
164 
168  float Sstim;
169 
173  float rel_diff;
174 
178  float rel_err;
179 };
180 
182 
187 fmriResult fmri_eval(const Data<float,1>& timecourse, const Data<float,1>& designvec);
188 
189 
193 #endif
194 
correlationResult correlation(const Array< T, N_rank > &x, const Array< T, N_rank > &y)
Definition: correlation.h:61
fmriResult fmri_eval(const Data< float, 1 > &timecourse, const Data< float, 1 > &designvec)
Definition: tjlog.h:218
correlationResult kendall(const Array< T, N_rank > &x, const Array< T, N_rank > &y)
Definition: correlation.h:106
float Srest
Definition: correlation.h:163
float Sstim
Definition: correlation.h:168
float Sbaseline
Definition: correlation.h:158
float rel_err
Definition: correlation.h:178
statisticResult statistics(const Array< T, N_rank > &ensemble, const Array< T, N_rank > *mask=0)
Definition: statistics.h:79
double secureDivision(double numerator, double denominator)
float rel_diff
Definition: correlation.h:173