Anonymous user
Welch's t-test: Difference between revisions
Undo revision 255357 by Thundergnat (talk) This was removed for no reason. The other perl translation doesn't work on my computer.
m (→{{header|Stata}}: should be more accurate) |
(Undo revision 255357 by Thundergnat (talk) This was removed for no reason. The other perl translation doesn't work on my computer.) |
||
Line 535:
=={{header|Perl}}==
{{trans|C}}
This is a portable translation of C.
<lang perl>
#!/usr/bin/env perl
use strict; use warnings; use Cwd;# use feature 'say';
my $TOP_DIRECTORY = getcwd();
local $SIG{__WARN__} = sub {#kill the program if there are any warnings
my $message = shift;
my $fail_filename = "$TOP_DIRECTORY/$0.FAIL";
open my $fh, '>', $fail_filename or die "Can't write $fail_filename: $!";
printf $fh ("$message @ %s\n", getcwd());
close $fh;
die "$message\n";
};#http://perlmaven.com/how-to-capture-and-save-warnings-in-perl
use POSIX 'lgamma';
use List::Util 'sum';
sub calculate_Pvalue {
my $array1 = shift;
my $array2 = shift;
if (scalar @$array1 <= 1) {
return 1.0;
}
if (scalar @$array2 <= 1) {
return 1.0;
}
my $mean1 = sum(@{ $array1 });
$mean1 /= scalar @$array1;
my $mean2 = sum(@{ $array2 });
$mean2 /= scalar @$array2;
if ($mean1 == $mean2) {
return 1.0;
}
# say "mean1 = $mean1 mean2 = $mean2";
my $variance1 = 0.0;
my $variance2 = 0.0;
foreach my $x (@$array1) {
$variance1 += ($x-$mean1)*($x-$mean1);
}
foreach my $x (@$array2) {
$variance2 += ($x-$mean2)*($x-$mean2);
}
if (($variance1 == 0.0) && ($variance2 == 0.0)) {
return 1.0;
}
$variance1 = $variance1/(scalar @$array1-1);
$variance2 = $variance2/(scalar @$array2-1);
my $array1_size = scalar @$array1;
my $array2_size = scalar @$array2;
my $WELCH_T_STATISTIC = ($mean1-$mean2)/sqrt($variance1/$array1_size+$variance2/$array2_size);
my $DEGREES_OF_FREEDOM = (($variance1/$array1_size+$variance2/(scalar @$array2))**2)
/
(
($variance1*$variance1)/($array1_size*$array1_size*($array1_size-1))+
($variance2*$variance2)/($array2_size*$array2_size*($array2_size-1))
);
# say "$WELCH_T_STATISTIC $DEGREES_OF_FREEDOM";
my $A = $DEGREES_OF_FREEDOM/2;
my $x = $DEGREES_OF_FREEDOM/($WELCH_T_STATISTIC*$WELCH_T_STATISTIC+$DEGREES_OF_FREEDOM);
my $N = 65536;
my $h = $x/$N;
my $sum1 = 0.0;
my $sum2 = 0.0;
for(my $i = 0;$i < $N; $i++) {
$sum1 += (($h * $i + $h / 2.0)**($A-1))/(sqrt(1-($h * $i + $h / 2.0)));
$sum2 += (($h * $i)**($A-1))/(sqrt(1-$h * $i));
}
# say "sum1 = $sum1 sum2 = $sum2";
my $return_value = (($h / 6.0) * ((($x)**($A-1))/(sqrt(1-$x)) + 4.0*$sum1 + 2.0*$sum2))/(exp(lgamma($A)+0.57236494292470009-lgamma($A+0.5)));
if ($return_value > 1.0) {
return 1.0;
} else {
return $return_value;
}
}
my @d1 = (27.5,21.0,19.0,23.6,17.0,17.9,16.9,20.1,21.9,22.6,23.1,19.6,19.0,21.7,21.4);
my @d2 = (27.1,22.0,20.8,23.4,23.4,23.5,25.8,22.0,24.8,20.2,21.9,22.1,22.9,20.5,24.4);
my @d3 = (17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8);
my @d4 = (21.5,22.8,21.0,23.0,21.6,23.6,22.5,20.7,23.4,21.8,20.7,21.7,21.5,22.5,23.6,21.5,22.5,23.5,21.5,21.8);
my @d5 = (19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0);
my @d6 = (28.2,26.6,20.1,23.3,25.2,22.1,17.7,27.6,20.6,13.7,23.2,17.5,20.6,18.0,23.9,21.6,24.3,20.4,24.0,13.2);
my @d7 = (30.02,29.99,30.11,29.97,30.01,29.99);
my @d8 = (29.89,29.93,29.72,29.98,30.02,29.98);
my @x = (3.0,4.0,1.0,2.1);
my @y = (490.2,340.0,433.9);
my @v1 = (0.010268,0.000167,0.000167);
my @v2 = (0.159258,0.136278,0.122389);
my @z1 = (9/23.0,21/45.0,0/38.0);
my @z2 = (0/44.0,42/94.0,0/22.0);
printf("Test sets 1 p-value = %g\n",&calculate_Pvalue(\@d1,\@d2));
printf("Test sets 2 p-value = %g\n",&calculate_Pvalue(\@d3,\@d4));
printf("Test sets 3 p-value = %g\n",&calculate_Pvalue(\@d5,\@d6));
printf("Test sets 4 p-value = %g\n",&calculate_Pvalue(\@d7,\@d8));
printf("Test sets 5 p-value = %g\n",&calculate_Pvalue(\@x,\@y));
printf("Test sets 6 p-value = %g\n",&calculate_Pvalue(\@v1,\@v2));
printf("Test sets 7 p-value = %g\n",&calculate_Pvalue(\@z1,\@z2));
</lang>
<pre>
Test sets 1 p-value = 0.021378
Test sets 2 p-value = 0.148842
Test sets 3 p-value = 0.0359723
Test sets 4 p-value = 0.0907733
Test sets 5 p-value = 0.0107515
Test sets 6 p-value = 0.00339907
Test sets 7 p-value = 0.545267
</pre>
{{trans|Sidef}}
<lang perl>use utf8;
| |||