File Coverage

blib/lib/Statistics/RankCorrelation.pm

Criterion	Covered	Total	%
statement	148	150	98.6
branch	49	58	84.4
condition	14	21	66.6
subroutine	21	21	100.0
pod	16	16	100.0
total	248	266	93.2

package Statistics::RankCorrelation;
BEGIN {
  $Statistics::RankCorrelation::AUTHORITY = 'cpan:GENE';
}
# ABSTRACT: Compute the rank correlation between two vectors

use strict;
use warnings;

our $VERSION = '0.1204';

use Carp;

sub new {
    my $proto = shift;
    my $class = ref($proto) || $proto;
    my $self = {};
    bless $self, $class;
    $self->_init(@_);
    return $self;
}

sub _init {
    my $self = shift;

    # Handle vector and named parameters.
    while( my $arg = shift ) {
        if( ref $arg eq 'ARRAY' ) {
               if( !defined $self->x_data ) { $self->x_data( $arg ) }
            elsif( !defined $self->y_data ) { $self->y_data( $arg ) }
        }
        elsif( !ref $arg ) {
            my $v = shift;
            $self->{$arg} = defined $v ? $v : $arg;
        }
    }

    # Automatically compute the ranks if given data.
    if( $self->x_data && $self->y_data &&
        @{ $self->x_data } && @{ $self->y_data }
    ) {
        # "Co-normalize" the vectors if they are of unequal size.
        my( $x, $y ) = pad_vectors( $self->x_data, $self->y_data );

        # "Co-sort" the bivariate data set by the first one.
        ( $x, $y ) = co_sort( $x, $y ) if $self->{sorted};

        # Set the massaged data.
        $self->x_data( $x );
        $self->y_data( $y );

        # Set the size of the data vector.
        $self->size( scalar @{ $self->x_data } );

        # Set the ranks and ties of the vectors.
        ( $x, $y ) = rank( $self->x_data );
        $self->x_rank( $x );
        $self->x_ties( $y );
        ( $x, $y ) = rank( $self->y_data );
        $self->y_rank( $x );
        $self->y_ties( $y );
    }
}

sub size {
    my $self = shift;
    $self->{size} = shift if @_;
    return $self->{size};
}

sub x_data {
    my $self = shift;
    $self->{x_data} = shift if @_;
    return $self->{x_data};
}

sub y_data {
    my $self = shift;
    $self->{y_data} = shift if @_;
    return $self->{y_data};
}

sub x_rank {
    my $self = shift;
    $self->{x_rank} = shift if @_;
    return $self->{x_rank};
}

sub y_rank {
    my $self = shift;
    $self->{y_rank} = shift if @_;
    return $self->{y_rank};
}

sub x_ties {
    my $self = shift;
    $self->{x_ties} = shift if @_;
    return $self->{x_ties};
}

sub y_ties {
    my $self = shift;
    $self->{y_ties} = shift if @_;
    return $self->{y_ties};
}

sub spearman {
    my $self = shift;
    # Algorithm contributed by Jon Schutz :
    my($x_sum, $y_sum) = (0, 0);
    $x_sum += $_ for @{$self->{x_rank}};
    $y_sum += $_ for @{$self->{y_rank}};
    my $n = $self->size;
    my $x_mean = $x_sum / $n;
    my $y_mean = $y_sum / $n;
    # Compute the sum of the difference of the squared ranks.
    my($x_sum2, $y_sum2, $xy_sum) = (0, 0, 0);
    for( 0 .. $self->size - 1 ) {
        $x_sum2 += ($self->{x_rank}[$_] - $x_mean) ** 2;
        $y_sum2 += ($self->{y_rank}[$_] - $y_mean) ** 2;
        $xy_sum += ($self->{x_rank}[$_] - $x_mean) * ($self->{y_rank}[$_] - $y_mean);
    }
    return 1 if $x_sum2 == 0 || $y_sum2 == 0;
    return $xy_sum / sqrt($x_sum2 * $y_sum2);
}


sub rank {
    my $u = shift;

    # Make a list of tied ranks for each datum.
    my %ties;
    push @{ $ties{ $u->[$_] } }, $_ for 0 .. @$u - 1;

    my ($old, $cur) = (0, 0);

    # Set the averaged ranks.
    my @ranks;
    for my $x (sort { $a <=> $b } keys %ties) {
        # Get the number of ties.
        my $ties = @{ $ties{$x} };
        $cur += $ties;

        if ($ties > 1) {
            # Average the tied data.
            my $average = $old + ($ties + 1) / 2;
            $ranks[$_] = $average for @{ $ties{$x} };
        }
        else {
            # Add the single rank to the list of ranks.
            $ranks[ $ties{$x}[0] ] = $cur;
        }

        $old = $cur;
    }

    # Remove the non-tied ranks.
    while( my( $k, $v ) = each %ties ) {
        delete $ties{$k} unless @$v > 1;
    }

    # Return the ranks arrayref in a scalar context and include ties
    # if called in a list context.
    return wantarray ? (\@ranks, \%ties) : \@ranks;
}

sub co_sort {
    my( $u, $v ) = @_;
    return unless @$u == @$v;
    # Ye olde Schwartzian Transforme:
    $v = [
        map { $_->[1] }
            sort { $a->[0] <=> $b->[0] || $a->[1] <=> $b->[1] }
                map { [ $u->[$_], $v->[$_] ] }
                    0 .. @$u - 1
    ];
    # Sort the independent vector last.
    $u = [ sort { $a <=> $b } @$u ];
    return $u, $v;
}

sub csim {
    my $self = shift;

    # Get the pitch matrices for each vector.
    my $m1 = correlation_matrix($self->{x_data});
#warn map { "@$_\n" } @$m1;
    my $m2 = correlation_matrix($self->{y_data});
#warn map { "@$_\n" } @$m2;

    # Compute the rank correlation.
    my $k = 0;
    for my $i (0 .. @$m1 - 1) {
        for my $j (0 .. @$m1 - 1) {
            $k++ if $m1->[$i][$j] == $m2->[$i][$j];
        }
    }

    # Return the rank correlation normalized by the number of rows in
    # the pitch matrices.
    return $k / (@$m1 * @$m1);
}

sub pad_vectors {
    my ($u, $v) = @_;

    if (@$u > @$v) {
        $v = [ @$v, (0) x (@$u - @$v) ];
    }
    elsif (@$u < @$v) {
        $u = [ @$u, (0) x (@$v - @$u) ];
    }

    return $u, $v;
}

sub correlation_matrix {
    my $u = shift;
    my $c;

    # Is a row value (i) lower than a column value (j)?
    for my $i (0 .. @$u - 1) {
        for my $j (0 .. @$u - 1) {
            $c->[$i][$j] = $u->[$i] < $u->[$j] ? 1 : 0;
        }
    }

    return $c;
}

sub kendall {
    my $self = shift;

    # Calculate number of concordant and discordant pairs.
    my( $concordant, $discordant ) = ( 0, 0 );
    for my $i ( 0 .. $self->size - 1 ) {
        for my $j ( $i + 1 .. $self->size - 1 ) {
            my $x_sign = sign( $self->{x_data}[$j] - $self->{x_data}[$i] );
            my $y_sign = sign( $self->{y_data}[$j] - $self->{y_data}[$i] );
            if (not($x_sign and $y_sign)) {}
            elsif ($x_sign == $y_sign) { $concordant++ }
            else { $discordant++ }
        }
    }

    # Set the indirect relationship.
    my $d = $self->size * ($self->size - 1) / 2;
    if( keys %{ $self->x_ties } || keys %{ $self->y_ties } ) {
        my $x = 0;
        $x += @$_ * (@$_ - 1) for values %{ $self->x_ties };
        $x = $d - $x / 2;
        my $y = 0;
        $y += @$_ * (@$_ - 1) for values %{ $self->y_ties };
        $y = $d - $y / 2;
        $d = sqrt($x * $y);
    }

    return ($concordant - $discordant) / $d;
}

sub sign {
    my $x = shift;
    return 0 if $x == 0;
    return $x > 0 ? 1 : -1;
}

1;

__END__