Statistics::Test::WilcoxonRankSum - perform the Wilcoxon (aka Mann-Whitney) rank sum test on two sets of numeric data.

    use Statistics::Test::WilcoxonRankSum;

    my $wilcox_test = Statistics::Test::WilcoxonRankSum->new();

    my @dataset_1 = (4.6, 4.7, 4.9, 5.1, 5.2, 5.5, 5.8, 6.1, 6.5, 6.5, 7.2);
    my @dataset_2 = (5.2, 5.3, 5.4, 5.6, 6.2, 6.3, 6.8, 7.7, 8.0, 8.1);

    $wilcox_test->load_data(\@dataset_1, \@dataset_2);
    my $prob = $wilcox_test->probability();

    my $pf = sprintf '%f', $prob; # prints 0.091022

    print $wilcox_test->probability_status();

    # prints something like:
    # Probability:   0.002797, exact
    # or
    # Probability:   0.511020, normal approx w. mean: 104.000000, std deviation:  41.840969, z:   0.657251

    my $pstatus = $wilcox_test->probability_status();
    # $pstatus is like the strings above

    $wilcox_test->summary();

    # prints something like:

    # ----------------------------------------------------------------
    # dataset |    n      | rank sum: observed / expected 
    # ----------------------------------------------------------------
    #   1    |     10    |               533      /    300
    # ----------------------------------------------------------------
    #   2    |     50    |              1296      /   1500
    # ----------------------------------------------------------------
    # N (size of both datasets):      60
    # Probability:   0.000006, normal approx w. mean: 305.000000, std deviation:  50.414945, z:   4.522468
    # Significant (at 0.05 level)
    # Ranks of dataset 1 are higher than expected

package Statistics::Test::WilcoxonRankSum;

use warnings;
use strict;
use Carp;
use Carp::Assert;

use version; our $VERSION = qv('0.0.7');

use Contextual::Return;
use List::Util qw(sum);
use Set::Partition;
use Math::BigFloat;
use Math::Counting ':big';
use Statistics::Distributions;

use Class::Std;

{
  ############ Data ######################################################################

  my %EXACT_UPTO  : ATTR( :init_arg<exact_upto> :default<20> );
  my %dataset1_of : ATTR( :get<dataset1> ); # array of numbers
  my %dataset2_of : ATTR( :get<dataset2> ); # array of numbers
  my %n1_of       : ATTR( :get<n1> );       # number of elements in dataset 1
  my %n2_of       : ATTR( :get<n2> );       # number of elements in dataset 2
  my %N_of        : ATTR( :get<N>  );       # overall number of elements (ranks)
  my %MaxSum_of   : ATTR( :get<max_rank_sum> ); # biggest possible ranksum
  my %ranks_of    : ATTR( :get<ranks>    :set<ranks>  ); # hash with ranked data
  my %rank_array_of : ATTR( :get<rank_array> ); # rank array from %ranks
  my %rankSum1_of : ATTR( :get<rankSum_dataset1> ); # rank sum for dataset 1
  my %expected_rank_sum_1_of : ATTR( :get<expected_rank_sum_dataset1>); # expected rank sum for dataset 1
  my %expected_rank_sum_2_of : ATTR( :get<expected_rank_sum_dataset2>); # expected rank sum for dataset 2
  my %rankSum2_of : ATTR( :get<rankSum_dataset2> ); # rank sum for dataset 2
  my %smaller_rank_sum_of : ATTR;
  my %smaller_ranks_count_of : ATTR;
  my %expected_rank_count_for_smaller_ranks_count_of : ATTR( :get<expected_rank_count_for_smaller_ranks_count>);
  my %smaller_rank_sums_count_of : ATTR;            # number of possible arrangements with lesser rank sum
                                                    # than the smaller rank sum
  my %rank_sums_other_than_expected_count_of : ATTR; # number of possible arrangements with rank sum
                                                    # other than the smaller rank sum
  my %probability_of : ATTR;                   # probability for the ranking with smaller rank sum
  my %probability_normal_approx_of : ATTR;
  
  ############ Utility subroutines #######################################################

  sub _check_dataset {
    my ($dataset_ref) = @_;

    croak "Need array ref to dataset\n"
      unless ($dataset_ref);

    croak "Datasets must be passed as array references\n"
      unless (ref($dataset_ref) eq 'ARRAY');
    
    my @dataset = grep { $_ > 0 } @{ $dataset_ref };
    croak "dataset has no element greater 0\n" unless (@dataset);

    return \@dataset;

  }

  sub _compute_N_MaxSum {
    my ($id) = @_;

    my $N;
    unless ($N_of{$id}) {
      $N = $n1_of{$id} + $n2_of{$id};
      $N_of{$id} = $N;
    }

    unless ($MaxSum_of{$id}) {
      $MaxSum_of{$id} = $N*($N+1)/2;
    }

    unless ($expected_rank_sum_1_of{$id}) {
      $expected_rank_sum_1_of{$id} = $n1_of{$id}*$N/2;
    }

    unless ($expected_rank_sum_2_of{$id}) {
      $expected_rank_sum_2_of{$id} = $n2_of{$id}*$N/2;
    }

    return;
  }

  

  sub _reset_dependant_datastructures {
    my ($id) = @_;

    delete $ranks_of{$id};
    delete $rank_array_of{$id};
    delete $rankSum1_of{$id};
    delete $rankSum2_of{$id};
    delete $N_of{$id};
    delete $MaxSum_of{$id};
    delete $smaller_rank_sum_of{$id};
    delete $smaller_rank_sums_count_of{$id};
    delete $probability_of{$id};
    delete $probability_normal_approx_of{$id};
    delete $expected_rank_sum_1_of{$id};
    delete $expected_rank_sum_2_of{$id};
    
    return;
  }

  sub _rank_sum_for {
    my ($self, $dataset) = @_;

    my $id = ident $self;
    
    my @rank_array;

    if ($rank_array_of{$id} and @{ $rank_array_of{$id} }) {
      @rank_array = @{ $rank_array_of{$id} };
    } else {
      @rank_array = $self->compute_rank_array();
    }
    
    return sum map { $_->[0] } grep { $_->[1] eq $dataset } @rank_array;
  }


  sub _set_smaller_rank_for {
    my ($id, $rank_sum_1, $rank_sum_2) = @_;
    if ($rank_sum_1 <= $rank_sum_2) {
      $smaller_rank_sum_of{$id} = $rank_sum_1;
      $smaller_ranks_count_of{$id} = $n1_of{$id};
      $expected_rank_count_for_smaller_ranks_count_of{$id} = $expected_rank_sum_1_of{$id};
    } else {
      $smaller_rank_sum_of{$id} = $rank_sum_2;
      $smaller_ranks_count_of{$id} = $n2_of{$id};
      $expected_rank_count_for_smaller_ranks_count_of{$id} = $expected_rank_sum_2_of{$id};
    };
    return;
  };

  sub _NormalZ {  # ($Z) -> $p
    my ($x) = @_;
    #
    # P(x) = 1 - Z(x)(b1*t+b2*t**2+b3*t**3+b4*t**4+b5*t**5)
    # Z(x) = exp(-$x*$x/2.0)/(sqrt(2*3.14159265358979323846))
    # t = 1/(1+p*x)
    #
    # Parameters
    my @b = (0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429);
    my $p = 0.2316419;
    my $t = 1/(1+$p*$x);
    # Initialize variables
    my $fact = $t;
    my $Sum;
    # Sum polynomial
    foreach  my $bi (@b) { 
      $Sum += $bi*$fact;
      $fact *= $t;
    };
    # Calculate probability
    $p = 2*$Sum*exp(-$x*$x/2.0)/(sqrt(2*3.14159265358979323846));
    #
    return $p;
  };
  
  ############ Methods ###################################################################

  sub set_dataset1 {
    my ($self, $dataset1_ref) = @_;

    $dataset1_ref = _check_dataset($dataset1_ref);
    
    my $id = ident $self;
    $dataset1_of{$id} = $dataset1_ref;
    $n1_of{$id} = scalar(@{ $dataset1_ref });
    
    _reset_dependant_datastructures($id);

    return;
  }

  sub set_dataset2 {
    my ($self, $dataset2_ref) = @_;
    
    $dataset2_ref = _check_dataset($dataset2_ref);
    
    my $id = ident $self;
    $dataset1_of{$id} = $dataset2_ref;
    $n2_of{$id} = scalar(@{ $dataset2_ref });
    
    _reset_dependant_datastructures($id);

    return;
  }

  sub load_data {
    my ($self, $dataset1_ref, $dataset2_ref) = @_;

    $dataset1_ref = _check_dataset($dataset1_ref);
    $dataset2_ref = _check_dataset($dataset2_ref);
    
    my $id = ident $self;

    $dataset1_of{$id} = $dataset1_ref;
    $dataset2_of{$id} = $dataset2_ref;
    $n1_of{$id} = scalar(@{ $dataset1_ref });
    $n2_of{$id} = scalar(@{ $dataset2_ref });
    
    _reset_dependant_datastructures($id);

    _compute_N_MaxSum($id);
    
    return;
  }

  sub compute_ranks {
    my ($self) = @_;
    my $id = ident $self;

    croak "Please set/load datasets before computing ranks\n" unless ($dataset1_of{$id} and $dataset2_of{$id});
    
    my @dataset1 = @{ $dataset1_of{$id} };
    my @dataset2 = @{ $dataset2_of{$id} };

    # at this point we are sure we have both data sets, so we may as well compute N and MaxSum - if not already computed
    _compute_N_MaxSum($id);
    
    my %ranks;

    foreach my $el (@dataset1) {
      $ranks{$el}->{in_dataset}->{ds1}++;
    }
    foreach my $el (@dataset2) {
      $ranks{$el}->{in_dataset}->{ds2}++;
    }

    my $rank=0;
    foreach my $value (sort { $a <=> $b } keys %ranks) {

      my $tied_ranks;

      foreach my $ds (keys %{ $ranks{$value}->{in_dataset} }) {
        $tied_ranks += $ranks{$value}->{in_dataset}->{$ds};
      }

      assert $tied_ranks if DEBUG;

      my $rs;
      for my $r ($rank+1 .. $rank+$tied_ranks) {
        $rs += $r;
      }
      $ranks{$value}->{rank} = $rs/$tied_ranks;
      $ranks{$value}->{tied} = $tied_ranks;

      $rank+=$tied_ranks;
    }
    
    $ranks_of{$id} = \%ranks;

    return $ranks_of{$id};
  }

  sub compute_rank_array {
    my ($self) = @_;
    my $id = ident $self;

    my @rank_array;
    if ($rank_array_of{$id} and @{ $rank_array_of{$id} } ) {
      @rank_array = @{ $rank_array_of{$id} };
    } else {

      my %ranks;

      if ($ranks_of{$id} and %{ $ranks_of{$id} } ) {
        %ranks = %{ $ranks_of{$id} };
      } else {
        %ranks = %{ $self->compute_ranks() };
      }

      foreach my $value (sort { $a <=> $b } keys %ranks) {
        foreach my $ds (keys %{ $ranks{$value}->{in_dataset} }) {
          for (1..$ranks{$value}->{in_dataset}->{$ds}) {
            push(@rank_array, [ $ranks{$value}->{rank}, $ds ]);
          }
        }
      }

      $rank_array_of{$id} = \@rank_array;

    }
    
    return (
            SCALAR { scalar @rank_array }      # How many?
            LIST   { @rank_array        }      # What are they?
           );
  }

  sub rank_sum_for {
    my ($self, $for_dataset) = @_;

    my $id = ident $self;

    my $rankSum;
    if ($for_dataset =~ m{1}) {
      if ($rankSum1_of{$id}) {
        return $rankSum1_of{$id};
      } else {
        $rankSum1_of{$id} = $self->_rank_sum_for('ds1');
        return $rankSum1_of{$id};
      }
    } elsif ($for_dataset =~ m{2}) {
      if ($rankSum2_of{$id}) {
        return $rankSum2_of{$id};
      } else {
        $rankSum2_of{$id} = $self->_rank_sum_for('ds2');
        return $rankSum2_of{$id};
      }
    } else {
      croak "Argument must match `1' or `2' (meaning dataset 1 or 2)\n";
    }

    return;

  }

  sub get_smaller_rank_sum {
    my ($self) = @_;

    my $id = ident $self;

    if ($smaller_rank_sum_of{$id} and $smaller_ranks_count_of{$id}) {

      return (
              SCALAR { $smaller_rank_sum_of{$id} }                                  # only the rank sum itselt
              LIST   { ($smaller_rank_sum_of{$id}, $smaller_ranks_count_of{$id} ) } # also the size of the corresponding ds
             );
    }
    
    my $rank_sum_1 = $rankSum1_of{$id};
    my $rank_sum_2 = $rankSum2_of{$id};

    if (not($rank_sum_1) and not($rank_sum_2)) {
      $rank_sum_1 = $self->rank_sum_for('ds1');
    }

    
    if ($rank_sum_1 and $rank_sum_2) {

      _set_smaller_rank_for($id, $rank_sum_1, $rank_sum_2);
      
    } elsif ($rank_sum_1) {
      $rank_sum_2 = $MaxSum_of{$id} - $rank_sum_1;
      $rankSum2_of{$id} = $rank_sum_2;

      _set_smaller_rank_for($id, $rank_sum_1, $rank_sum_2);

    } elsif ($rank_sum_2) {
      $rank_sum_1 = $MaxSum_of{$id} - $rank_sum_2;
      $rankSum1_of{$id} = $rank_sum_1;

      _set_smaller_rank_for($id, $rank_sum_1, $rank_sum_2);

    }

    return (
            SCALAR { $smaller_rank_sum_of{$id} }                                  # only the rank sum itselt
            LIST   { ($smaller_rank_sum_of{$id}, $smaller_ranks_count_of{$id} ) } # also the size of the corresponding ds
           );

    
    return $smaller_rank_sum_of{$id};
  }
  
  sub smaller_rank_sums_count {
    my ($self) = @_;
    my $id = ident $self;

    if ($smaller_rank_sums_count_of{$id}) {
      return $smaller_rank_sums_count_of{$id};
    };
    
    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $N = $N_of{$id};
    my $nB = $N - $nA;
    my $MaxSum = $MaxSum_of{$id};

    my @ranks = map { $_->[0] } $self->compute_rank_array();

    # let's do some checks before starting the big counting
    if ($W > $MaxSum) { croak "Rank sum bound $W is bigger than the maximum possible rank sum $MaxSum\n" };
    if ($N != scalar(@ranks))
      { croak "Sum of $nA and $nB must be equal to number of ranks: ".scalar(@ranks)."\n" }; 

    # compute all possible partitions
    my $s = Set::Partition->new(
                                list => \@ranks,
                                partition => [$nA, $nB],
                               );

    my $count_less_W = 0;
    
    while (my $p = $s->next()) {
      my @pA = @{ $p->[0] };
      my $sumA = sum @pA;
      if ($sumA <= $W) {
        $count_less_W++;
      }
    }
    
    return $count_less_W;
    
  };

  sub rank_sums_other_than_expected_counts {
    my ($self) = @_;
    my $id = ident $self;

    if ($rank_sums_other_than_expected_count_of{$id}) {
      return $rank_sums_other_than_expected_count_of{$id};
    };

    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $W_exp = $self->get_expected_rank_count_for_smaller_ranks_count();

    my $N = $N_of{$id};
    my $nB = $N - $nA;
    my $MaxSum = $MaxSum_of{$id};

    my @ranks = map { $_->[0] } $self->compute_rank_array();

    # let's do some checks before starting the big counting
    if ($W > $MaxSum) { croak "Rank sum bound $W is bigger than the maximum possible rank sum $MaxSum\n" };
    if ($N != scalar(@ranks))
      { croak "Sum of $nA and $nB must be equal to number of ranks: ".scalar(@ranks)."\n" }; 

    # compute all possible partitions
    my $s = Set::Partition->new(
                                list => \@ranks,
                                partition => [$nA, $nB],
                               );

    my $count_other_W = 0;

    if ($W >= $W_exp) { 
    
      while (my $p = $s->next()) {
	my @pA = @{ $p->[0] };
	my $sumA = sum @pA;
	if ($sumA >= $W) {
	  $count_other_W++;
	}
      }

    } else {
    
      while (my $p = $s->next()) {
	my @pA = @{ $p->[0] };
	my $sumA = sum @pA;
	if ($sumA <= $W) {
	  $count_other_W++;
	}
      }

    }
    
    return $count_other_W;
    
  };


  sub probability :NUMERIFY {
    my ($self) = @_;
    my $id = ident $self;

    if ($probability_of{$id}) {
      return $probability_of{$id};
    }

    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $N = $N_of{$id};

    my $p;
    if ($N <= $EXACT_UPTO{$id}) {
      $p = $self->probability_exact();
    } else {
      $p = $self->probability_normal_approx();
    }

    $probability_of{$id} = $p;

    return $probability_of{$id};
  }

  sub probability_exact {
    my ($self) = @_;
    my $id = ident $self;
    
    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $N = $N_of{$id};

    my $partition_count = bcomb($N, $nA);
    my $have_smaller_rank_sums = $self->rank_sums_other_than_expected_counts();
    my $p = Math::BigFloat->new($have_smaller_rank_sums) * 2.0 / Math::BigFloat->new($partition_count);

    if ($p > 1) { $p = 1 };

    return $p;
  }

  sub probability_normal_approx {
    my ($self) = @_;
    my $id = ident $self;
    
    my ($W, $nA) = $self->get_smaller_rank_sum();
    my $N = $N_of{$id};
    my $nB = $N - $nA;
    my $mean = $nA*($N+1)/2;
    my $deviation = sqrt($nA*$nB*($N+1)/12.0);
    my $continuity = (($W - $mean) >= 0) ? -0.5 : +0.5;
    my $z = ($W - $mean + $continuity)/$deviation;
    @{ $probability_normal_approx_of{$id} }{'mean', 'std deviation', 'z'} = ($mean, $deviation, $z);
    
    my $p = 2*Statistics::Distributions::uprob(abs($z));

    return $p;
    
  }
  

  sub probability_status {
    my ($self) = (@_);
    my $id = ident $self;

    my $return_string;
    if ($probability_of{$id}) {
      if ($probability_normal_approx_of{$id}) {
        $return_string = sprintf "Probability: %10f, normal approx w. mean: %10f, std deviation: %10f, z: %10f", $probability_of{$id}, map { $probability_normal_approx_of{$id}->{$_} } ('mean', 'std deviation', 'z');
      } else {
        $return_string = sprintf "Probability: %10f, exact", $probability_of{$id};
      }
    } else {
      $return_string = "Probability not yet computed";
    }

    return (
            STR    { "$return_string" }
            VOID   { print $return_string."\n" }
           );
  }

  sub as_hash :HASHIFY {
    my ($self) = @_;
    my $id = ident $self;

    return {
            dataset_1 => $dataset1_of{$id},
            dataset_2 => $dataset2_of{$id},
            n1 => $n1_of{$id},
            n2 => $n2_of{$id},
            N => $N_of{$id},
            rank_array => $rank_array_of{$id},
            rank_sum_1 => $rankSum1_of{$id},
            rank_sum_2 => $rankSum2_of{$id},
            rank_sum_1_expected => $expected_rank_sum_1_of{$id},
            rank_sum_2_expected => $expected_rank_sum_2_of{$id},
            probability => $probability_of{$id},
            probability_normal_approx => $probability_normal_approx_of{$id},
           };

  }

  
  sub summary :STRINGIFY {
    my ($self) = (@_);
    my $id = ident $self;

    my $hash = $self->as_hash();

    my $return_string;
    if (not($hash->{dataset_1})) {
      $return_string = "Dataset 1 is not yet initialised, no computations could be done\n";
    } elsif (not($hash->{dataset_2})) {
      $return_string = "Dataset 2 is not yet initialised, no computations could be done\n";
    } else {
      my $format = <<END_FORMAT;
----------------------------------------------------------------
dataset |    n      | rank sum: observed / expected 
----------------------------------------------------------------
   1    |%7d    |           %7d      /%7d
----------------------------------------------------------------
   2    |%7d    |           %7d      /%7d
----------------------------------------------------------------
N (size of both datasets): %7d
%s
END_FORMAT
      my $prob = $self->probability_status();
      $return_string = sprintf $format, @{ $hash }{'n1', 'rank_sum_1', 'rank_sum_1_expected', 'n2', 'rank_sum_2', 'rank_sum_2_expected', 'N'}, $prob;
      if ($hash->{probability} >= 0.05) {
        $return_string.="Not significant (at 0.05 level)\n";
      } else {
        $return_string.="Significant (at 0.05 level)\n";
        $return_string.= $hash->{rank_sum_1} > $hash->{rank_sum_1_expected} ?
          "Ranks of dataset 1 are higher than expected\n"
            : "Ranks of dataset 1 are lower than expected\n";
          
      }
      if ($hash->{N} < 5) {
        $return_string.="Warning: sample size ($hash->{N}) too small (<5)!\n";
      }
    }

    return (
            STR    { "$return_string" }
            VOID   { print $return_string }
           );
  }
  
}

1; # Magic true value required at end of module


__END__