sign_test.pl

# Reference:
#   http://aoki2.si.gunma-u.ac.jp/lecture/Average/sign-test.html
#   http://kusuri-jouhou.com/statistics/fugou.html

use strict;
use warnings;

use List::Util;

my @sign = qw(1 0 1 1 -1 1 1 -1 1 1);
#my @sign = qw(1 -1 1 1 -1 1 1 0 -1 1 1 0);
#my @sign = qw(1 -1 1 1 -1 1 1 0 -1 1 1 -1 1 1 -1 1 1 1 -1 1 -1 1 1 -1 -1 -1 1);

printf("N = %d, P = %f", sign_test(@sign));

sub sign_test {
	my @list = @_;

	my @sign;
	foreach my $val (@_) {
		if ($val != 0) {
			push @sign, $val;
		}
	}

	if (scalar(@sign) <= 5) {
		return scalar(@sign), -1;
	} elsif (scalar(@sign) <= 25) {
		return scalar(@sign), sign_test_small(@sign);
	} else {
		return scalar(@sign), sign_test_large(@sign);
	}
}

sub sign_test_small {
	my @list = @_;

	my $is_plus  = 0;
	my $is_minus = 0;
	foreach my $val (@list) {
		if ($val > 0) {
			$is_plus++;
		} elsif ($val < 0) {
			$is_minus++;
		}
	}

	my $p = 0;
	foreach my $i (0 ... List::Util::min($is_plus, $is_minus)) {
		$p += combination($is_plus + $is_minus, $i) * (1 / 2) ** ($is_plus + $is_minus);
	}

	return $p * 2;
}

sub sign_test_large {
	my @list = @_;

	my $is_plus  = 0;
	my $is_minus = 0;
	foreach my $val (@list) {
		if ($val > 0) {
			$is_plus++;
		} elsif ($val < 0) {
			$is_minus++;
		}
	}

	my $z = (abs($is_plus - $is_minus) - 1) / sqrt($is_plus + $is_minus);

	return uprob($z) * 2;
}

# from Statistics::Distributions
#   http://search.cpan.org/perldoc?Statistics%3A%3ADistributions
use Math::Trig;
sub uprob {
	my ($x) = @_;
	my $p = 0; # if ($absx > 100)
	my $absx = abs($x);

	if ($absx < 1.9) {
		$p = (1 +
			$absx * (.049867347
			  + $absx * (.0211410061
			  	+ $absx * (.0032776263
				  + $absx * (.0000380036
					+ $absx * (.0000488906
					  + $absx * .000005383)))))) ** -16/2;
	} elsif ($absx <= 100) {
		for (my $i = 18; $i >= 1; $i--) {
			$p = $i / ($absx + $p);
		}
		$p = exp(-.5 * $absx * $absx) 
			/ sqrt(2 * pi) / ($absx + $p);
	}

	$p = 1 - $p if ($x<0);
	return $p;
}

sub combination {
	my ($n, $m) = @_;

	my $numerator   = 1;
	my $denominator = 1;
	foreach my $k (1 ... $m) {
		$numerator   *= ($n - $k + 1);
		$denominator *= $k;
	}

	return $numerator / $denominator;
}