Skip to content

Instantly share code, notes, and snippets.

@ceekz
Last active Dec 15, 2015
Embed
What would you like to do?
# 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment