line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
package Statistics::RankCorrelation; |
2
|
|
|
|
|
|
|
our $AUTHORITY = 'cpan:GENE'; |
3
|
|
|
|
|
|
|
# ABSTRACT: Compute the rank correlation between two vectors |
4
|
|
|
|
|
|
|
|
5
|
1
|
|
|
1
|
|
72378
|
use strict; |
|
1
|
|
|
|
|
12
|
|
|
1
|
|
|
|
|
32
|
|
6
|
1
|
|
|
1
|
|
6
|
use warnings; |
|
1
|
|
|
|
|
2
|
|
|
1
|
|
|
|
|
1805
|
|
7
|
|
|
|
|
|
|
|
8
|
|
|
|
|
|
|
our $VERSION = '0.1206'; |
9
|
|
|
|
|
|
|
|
10
|
|
|
|
|
|
|
sub new { |
11
|
15
|
|
|
15
|
1
|
806
|
my $proto = shift; |
12
|
15
|
|
33
|
|
|
60
|
my $class = ref($proto) || $proto; |
13
|
15
|
|
|
|
|
28
|
my $self = {}; |
14
|
15
|
|
|
|
|
28
|
bless $self, $class; |
15
|
15
|
|
|
|
|
47
|
$self->_init(@_); |
16
|
15
|
|
|
|
|
84
|
return $self; |
17
|
|
|
|
|
|
|
} |
18
|
|
|
|
|
|
|
|
19
|
|
|
|
|
|
|
sub _init { |
20
|
15
|
|
|
15
|
|
27
|
my $self = shift; |
21
|
|
|
|
|
|
|
|
22
|
|
|
|
|
|
|
# Handle vector and named parameters. |
23
|
15
|
|
|
|
|
36
|
while( my $arg = shift ) { |
24
|
29
|
100
|
|
|
|
77
|
if( ref $arg eq 'ARRAY' ) { |
|
|
50
|
|
|
|
|
|
25
|
28
|
100
|
|
|
|
50
|
if( !defined $self->x_data ) { $self->x_data( $arg ) } |
|
14
|
50
|
|
|
|
28
|
|
26
|
14
|
|
|
|
|
28
|
elsif( !defined $self->y_data ) { $self->y_data( $arg ) } |
27
|
|
|
|
|
|
|
} |
28
|
|
|
|
|
|
|
elsif( !ref $arg ) { |
29
|
1
|
|
|
|
|
13
|
my $v = shift; |
30
|
1
|
50
|
|
|
|
7
|
$self->{$arg} = defined $v ? $v : $arg; |
31
|
|
|
|
|
|
|
} |
32
|
|
|
|
|
|
|
} |
33
|
|
|
|
|
|
|
|
34
|
|
|
|
|
|
|
# Automatically compute the ranks if given data. |
35
|
15
|
50
|
66
|
|
|
32
|
if( $self->x_data && $self->y_data && |
|
|
|
66
|
|
|
|
|
|
|
|
33
|
|
|
|
|
36
|
14
|
|
|
|
|
26
|
@{ $self->x_data } && @{ $self->y_data } |
|
14
|
|
|
|
|
25
|
|
37
|
|
|
|
|
|
|
) { |
38
|
|
|
|
|
|
|
# "Co-normalize" the vectors if they are of unequal size. |
39
|
14
|
|
|
|
|
25
|
my( $x, $y ) = pad_vectors( $self->x_data, $self->y_data ); |
40
|
|
|
|
|
|
|
|
41
|
|
|
|
|
|
|
# "Co-sort" the bivariate data set by the first one. |
42
|
14
|
100
|
|
|
|
36
|
( $x, $y ) = co_sort( $x, $y ) if $self->{sorted}; |
43
|
|
|
|
|
|
|
|
44
|
|
|
|
|
|
|
# Set the massaged data. |
45
|
14
|
|
|
|
|
36
|
$self->x_data( $x ); |
46
|
14
|
|
|
|
|
33
|
$self->y_data( $y ); |
47
|
|
|
|
|
|
|
|
48
|
|
|
|
|
|
|
# Set the size of the data vector. |
49
|
14
|
|
|
|
|
43
|
$self->size( scalar @{ $self->x_data } ); |
|
14
|
|
|
|
|
24
|
|
50
|
|
|
|
|
|
|
|
51
|
|
|
|
|
|
|
# Set the ranks and ties of the vectors. |
52
|
14
|
|
|
|
|
27
|
( $x, $y ) = rank( $self->x_data ); |
53
|
14
|
|
|
|
|
42
|
$self->x_rank( $x ); |
54
|
14
|
|
|
|
|
35
|
$self->x_ties( $y ); |
55
|
14
|
|
|
|
|
28
|
( $x, $y ) = rank( $self->y_data ); |
56
|
14
|
|
|
|
|
42
|
$self->y_rank( $x ); |
57
|
14
|
|
|
|
|
30
|
$self->y_ties( $y ); |
58
|
|
|
|
|
|
|
} |
59
|
|
|
|
|
|
|
} |
60
|
|
|
|
|
|
|
|
61
|
|
|
|
|
|
|
sub size { |
62
|
196
|
|
|
196
|
1
|
279
|
my $self = shift; |
63
|
196
|
100
|
|
|
|
370
|
$self->{size} = shift if @_; |
64
|
196
|
|
|
|
|
403
|
return $self->{size}; |
65
|
|
|
|
|
|
|
} |
66
|
|
|
|
|
|
|
|
67
|
|
|
|
|
|
|
sub x_data { |
68
|
130
|
|
|
130
|
1
|
926
|
my $self = shift; |
69
|
130
|
100
|
|
|
|
244
|
$self->{x_data} = shift if @_; |
70
|
130
|
|
|
|
|
319
|
return $self->{x_data}; |
71
|
|
|
|
|
|
|
} |
72
|
|
|
|
|
|
|
|
73
|
|
|
|
|
|
|
sub y_data { |
74
|
101
|
|
|
101
|
1
|
152
|
my $self = shift; |
75
|
101
|
100
|
|
|
|
183
|
$self->{y_data} = shift if @_; |
76
|
101
|
|
|
|
|
240
|
return $self->{y_data}; |
77
|
|
|
|
|
|
|
} |
78
|
|
|
|
|
|
|
|
79
|
|
|
|
|
|
|
sub x_rank { |
80
|
18
|
|
|
18
|
1
|
33
|
my $self = shift; |
81
|
18
|
100
|
|
|
|
43
|
$self->{x_rank} = shift if @_; |
82
|
18
|
|
|
|
|
42
|
return $self->{x_rank}; |
83
|
|
|
|
|
|
|
} |
84
|
|
|
|
|
|
|
|
85
|
|
|
|
|
|
|
sub y_rank { |
86
|
19
|
|
|
19
|
1
|
32
|
my $self = shift; |
87
|
19
|
100
|
|
|
|
52
|
$self->{y_rank} = shift if @_; |
88
|
19
|
|
|
|
|
45
|
return $self->{y_rank}; |
89
|
|
|
|
|
|
|
} |
90
|
|
|
|
|
|
|
|
91
|
|
|
|
|
|
|
sub x_ties { |
92
|
46
|
|
|
46
|
1
|
119
|
my $self = shift; |
93
|
46
|
100
|
|
|
|
94
|
$self->{x_ties} = shift if @_; |
94
|
46
|
|
|
|
|
171
|
return $self->{x_ties}; |
95
|
|
|
|
|
|
|
} |
96
|
|
|
|
|
|
|
|
97
|
|
|
|
|
|
|
sub y_ties { |
98
|
44
|
|
|
44
|
1
|
77
|
my $self = shift; |
99
|
44
|
100
|
|
|
|
98
|
$self->{y_ties} = shift if @_; |
100
|
44
|
|
|
|
|
180
|
return $self->{y_ties}; |
101
|
|
|
|
|
|
|
} |
102
|
|
|
|
|
|
|
|
103
|
|
|
|
|
|
|
sub spearman { |
104
|
14
|
|
|
14
|
1
|
29
|
my $self = shift; |
105
|
|
|
|
|
|
|
# Algorithm contributed by Jon Schutz : |
106
|
14
|
|
|
|
|
30
|
my($x_sum, $y_sum) = (0, 0); |
107
|
14
|
|
|
|
|
24
|
$x_sum += $_ for @{$self->{x_rank}}; |
|
14
|
|
|
|
|
61
|
|
108
|
14
|
|
|
|
|
24
|
$y_sum += $_ for @{$self->{y_rank}}; |
|
14
|
|
|
|
|
50
|
|
109
|
14
|
|
|
|
|
33
|
my $n = $self->size; |
110
|
14
|
|
|
|
|
32
|
my $x_mean = $x_sum / $n; |
111
|
14
|
|
|
|
|
24
|
my $y_mean = $y_sum / $n; |
112
|
|
|
|
|
|
|
# Compute the sum of the difference of the squared ranks. |
113
|
14
|
|
|
|
|
29
|
my($x_sum2, $y_sum2, $xy_sum) = (0, 0, 0); |
114
|
14
|
|
|
|
|
31
|
for( 0 .. $self->size - 1 ) { |
115
|
119
|
|
|
|
|
216
|
$x_sum2 += ($self->{x_rank}[$_] - $x_mean) ** 2; |
116
|
119
|
|
|
|
|
186
|
$y_sum2 += ($self->{y_rank}[$_] - $y_mean) ** 2; |
117
|
119
|
|
|
|
|
218
|
$xy_sum += ($self->{x_rank}[$_] - $x_mean) * ($self->{y_rank}[$_] - $y_mean); |
118
|
|
|
|
|
|
|
} |
119
|
14
|
100
|
66
|
|
|
67
|
return 1 if $x_sum2 == 0 || $y_sum2 == 0; |
120
|
13
|
|
|
|
|
172
|
return $xy_sum / sqrt($x_sum2 * $y_sum2); |
121
|
|
|
|
|
|
|
} |
122
|
|
|
|
|
|
|
|
123
|
|
|
|
|
|
|
|
124
|
|
|
|
|
|
|
sub rank { |
125
|
28
|
|
|
28
|
1
|
41
|
my $u = shift; |
126
|
|
|
|
|
|
|
|
127
|
|
|
|
|
|
|
# Make a list of tied ranks for each datum. |
128
|
28
|
|
|
|
|
45
|
my %ties; |
129
|
28
|
|
|
|
|
74
|
push @{ $ties{ $u->[$_] } }, $_ for 0 .. @$u - 1; |
|
238
|
|
|
|
|
609
|
|
130
|
|
|
|
|
|
|
|
131
|
28
|
|
|
|
|
53
|
my ($old, $cur) = (0, 0); |
132
|
|
|
|
|
|
|
|
133
|
|
|
|
|
|
|
# Set the averaged ranks. |
134
|
28
|
|
|
|
|
42
|
my @ranks; |
135
|
28
|
|
|
|
|
128
|
for my $x (sort { $a <=> $b } keys %ties) { |
|
381
|
|
|
|
|
603
|
|
136
|
|
|
|
|
|
|
# Get the number of ties. |
137
|
197
|
|
|
|
|
271
|
my $ties = @{ $ties{$x} }; |
|
197
|
|
|
|
|
321
|
|
138
|
197
|
|
|
|
|
264
|
$cur += $ties; |
139
|
|
|
|
|
|
|
|
140
|
197
|
100
|
|
|
|
344
|
if ($ties > 1) { |
141
|
|
|
|
|
|
|
# Average the tied data. |
142
|
15
|
|
|
|
|
30
|
my $average = $old + ($ties + 1) / 2; |
143
|
15
|
|
|
|
|
23
|
$ranks[$_] = $average for @{ $ties{$x} }; |
|
15
|
|
|
|
|
50
|
|
144
|
|
|
|
|
|
|
} |
145
|
|
|
|
|
|
|
else { |
146
|
|
|
|
|
|
|
# Add the single rank to the list of ranks. |
147
|
182
|
|
|
|
|
290
|
$ranks[ $ties{$x}[0] ] = $cur; |
148
|
|
|
|
|
|
|
} |
149
|
|
|
|
|
|
|
|
150
|
197
|
|
|
|
|
307
|
$old = $cur; |
151
|
|
|
|
|
|
|
} |
152
|
|
|
|
|
|
|
|
153
|
|
|
|
|
|
|
# Remove the non-tied ranks. |
154
|
28
|
|
|
|
|
79
|
delete @ties{ grep @{ $ties{$_} } <= 1, keys %ties }; |
|
197
|
|
|
|
|
345
|
|
155
|
|
|
|
|
|
|
|
156
|
|
|
|
|
|
|
# Return the ranks arrayref in a scalar context and include ties |
157
|
|
|
|
|
|
|
# if called in a list context. |
158
|
28
|
50
|
|
|
|
108
|
return wantarray ? (\@ranks, \%ties) : \@ranks; |
159
|
|
|
|
|
|
|
} |
160
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
sub co_sort { |
162
|
1
|
|
|
1
|
1
|
4
|
my( $u, $v ) = @_; |
163
|
1
|
50
|
|
|
|
4
|
return unless @$u == @$v; |
164
|
|
|
|
|
|
|
# Ye olde Schwartzian Transforme: |
165
|
|
|
|
|
|
|
$v = [ |
166
|
10
|
|
|
|
|
18
|
map { $_->[1] } |
167
|
23
|
50
|
|
|
|
49
|
sort { $a->[0] <=> $b->[0] || $a->[1] <=> $b->[1] } |
168
|
1
|
|
|
|
|
6
|
map { [ $u->[$_], $v->[$_] ] } |
|
10
|
|
|
|
|
34
|
|
169
|
|
|
|
|
|
|
0 .. @$u - 1 |
170
|
|
|
|
|
|
|
]; |
171
|
|
|
|
|
|
|
# Sort the independent vector last. |
172
|
1
|
|
|
|
|
7
|
$u = [ sort { $a <=> $b } @$u ]; |
|
23
|
|
|
|
|
37
|
|
173
|
1
|
|
|
|
|
4
|
return $u, $v; |
174
|
|
|
|
|
|
|
} |
175
|
|
|
|
|
|
|
|
176
|
|
|
|
|
|
|
sub csim { |
177
|
3
|
|
|
3
|
1
|
7
|
my $self = shift; |
178
|
|
|
|
|
|
|
|
179
|
|
|
|
|
|
|
# Get the pitch matrices for each vector. |
180
|
3
|
|
|
|
|
12
|
my $m1 = correlation_matrix($self->{x_data}); |
181
|
|
|
|
|
|
|
#warn map { "@$_\n" } @$m1; |
182
|
3
|
|
|
|
|
9
|
my $m2 = correlation_matrix($self->{y_data}); |
183
|
|
|
|
|
|
|
#warn map { "@$_\n" } @$m2; |
184
|
|
|
|
|
|
|
|
185
|
|
|
|
|
|
|
# Compute the rank correlation. |
186
|
3
|
|
|
|
|
7
|
my $k = 0; |
187
|
3
|
|
|
|
|
8
|
for my $i (0 .. @$m1 - 1) { |
188
|
16
|
|
|
|
|
29
|
for my $j (0 .. @$m1 - 1) { |
189
|
96
|
100
|
|
|
|
185
|
$k++ if $m1->[$i][$j] == $m2->[$i][$j]; |
190
|
|
|
|
|
|
|
} |
191
|
|
|
|
|
|
|
} |
192
|
|
|
|
|
|
|
|
193
|
|
|
|
|
|
|
# Return the rank correlation normalized by the number of rows in |
194
|
|
|
|
|
|
|
# the pitch matrices. |
195
|
3
|
|
|
|
|
21
|
return $k / (@$m1 * @$m1); |
196
|
|
|
|
|
|
|
} |
197
|
|
|
|
|
|
|
|
198
|
|
|
|
|
|
|
sub pad_vectors { |
199
|
14
|
|
|
14
|
1
|
29
|
my ($u, $v) = @_; |
200
|
|
|
|
|
|
|
|
201
|
14
|
50
|
|
|
|
41
|
if (@$u > @$v) { |
|
|
50
|
|
|
|
|
|
202
|
0
|
|
|
|
|
0
|
$v = [ @$v, (0) x (@$u - @$v) ]; |
203
|
|
|
|
|
|
|
} |
204
|
|
|
|
|
|
|
elsif (@$u < @$v) { |
205
|
0
|
|
|
|
|
0
|
$u = [ @$u, (0) x (@$v - @$u) ]; |
206
|
|
|
|
|
|
|
} |
207
|
|
|
|
|
|
|
|
208
|
14
|
|
|
|
|
32
|
return $u, $v; |
209
|
|
|
|
|
|
|
} |
210
|
|
|
|
|
|
|
|
211
|
|
|
|
|
|
|
sub correlation_matrix { |
212
|
6
|
|
|
6
|
1
|
12
|
my $u = shift; |
213
|
6
|
|
|
|
|
9
|
my $c; |
214
|
|
|
|
|
|
|
|
215
|
|
|
|
|
|
|
# Is a row value (i) lower than a column value (j)? |
216
|
6
|
|
|
|
|
15
|
for my $i (0 .. @$u - 1) { |
217
|
32
|
|
|
|
|
54
|
for my $j (0 .. @$u - 1) { |
218
|
192
|
100
|
|
|
|
379
|
$c->[$i][$j] = $u->[$i] < $u->[$j] ? 1 : 0; |
219
|
|
|
|
|
|
|
} |
220
|
|
|
|
|
|
|
} |
221
|
|
|
|
|
|
|
|
222
|
6
|
|
|
|
|
17
|
return $c; |
223
|
|
|
|
|
|
|
} |
224
|
|
|
|
|
|
|
|
225
|
|
|
|
|
|
|
sub kendall { |
226
|
13
|
|
|
13
|
1
|
30
|
my $self = shift; |
227
|
|
|
|
|
|
|
|
228
|
|
|
|
|
|
|
# Calculate number of concordant and discordant pairs. |
229
|
13
|
|
|
|
|
27
|
my( $concordant, $discordant ) = ( 0, 0 ); |
230
|
13
|
|
|
|
|
29
|
for my $i ( 0 .. $self->size - 1 ) { |
231
|
115
|
|
|
|
|
213
|
for my $j ( $i + 1 .. $self->size - 1 ) { |
232
|
574
|
|
|
|
|
1074
|
my $x_sign = sign( $self->{x_data}[$j] - $self->{x_data}[$i] ); |
233
|
574
|
|
|
|
|
1179
|
my $y_sign = sign( $self->{y_data}[$j] - $self->{y_data}[$i] ); |
234
|
574
|
100
|
100
|
|
|
1593
|
if (not($x_sign and $y_sign)) {} |
|
|
100
|
|
|
|
|
|
235
|
288
|
|
|
|
|
464
|
elsif ($x_sign == $y_sign) { $concordant++ } |
236
|
178
|
|
|
|
|
293
|
else { $discordant++ } |
237
|
|
|
|
|
|
|
} |
238
|
|
|
|
|
|
|
} |
239
|
|
|
|
|
|
|
|
240
|
|
|
|
|
|
|
# Set the indirect relationship. |
241
|
13
|
|
|
|
|
23
|
my $d = $self->size * ($self->size - 1) / 2; |
242
|
13
|
100
|
100
|
|
|
22
|
if( keys %{ $self->x_ties } || keys %{ $self->y_ties } ) { |
|
13
|
|
|
|
|
26
|
|
|
11
|
|
|
|
|
23
|
|
243
|
5
|
|
|
|
|
10
|
my $x = 0; |
244
|
5
|
|
|
|
|
12
|
$x += @$_ * (@$_ - 1) for values %{ $self->x_ties }; |
|
5
|
|
|
|
|
10
|
|
245
|
5
|
|
|
|
|
13
|
$x = $d - $x / 2; |
246
|
5
|
|
|
|
|
7
|
my $y = 0; |
247
|
5
|
|
|
|
|
8
|
$y += @$_ * (@$_ - 1) for values %{ $self->y_ties }; |
|
5
|
|
|
|
|
11
|
|
248
|
5
|
|
|
|
|
11
|
$y = $d - $y / 2; |
249
|
5
|
|
|
|
|
10
|
$d = sqrt($x * $y); |
250
|
|
|
|
|
|
|
} |
251
|
|
|
|
|
|
|
|
252
|
13
|
|
|
|
|
166
|
return ($concordant - $discordant) / $d; |
253
|
|
|
|
|
|
|
} |
254
|
|
|
|
|
|
|
|
255
|
|
|
|
|
|
|
sub sign { |
256
|
1148
|
|
|
1148
|
1
|
1599
|
my $x = shift; |
257
|
1148
|
100
|
|
|
|
1987
|
return 0 if $x == 0; |
258
|
1024
|
100
|
|
|
|
1752
|
return $x > 0 ? 1 : -1; |
259
|
|
|
|
|
|
|
} |
260
|
|
|
|
|
|
|
|
261
|
|
|
|
|
|
|
1; |
262
|
|
|
|
|
|
|
|
263
|
|
|
|
|
|
|
__END__ |