Welch's t-test: Difference between revisions
Content added Content deleted
(Added Julia language) |
(→{{header|Perl}}: I cannot install the Math::AnyNum package. I added my own Perl translation of C which should work on any perl installation.) |
||
Line 535: | Line 535: | ||
=={{header|Perl}}== |
=={{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}} |
{{trans|Sidef}} |
||
<lang perl>use utf8; |
<lang perl>use utf8; |