Welch's t-test: Difference between revisions
Added Perl
(Added Sidef) |
(Added Perl) |
||
Line 479:
0.090773
0.010751
</pre>
=={{header|Perl}}==
{{trans|Sidef}}
<lang perl>use utf8;
use List::Util qw(sum);
use Math::AnyNum qw(gamma pi);
sub p_value ($$) {
my ($A, $B) = @_;
(@$A > 1 && @$B > 1) || return 1;
my $x̄_a = sum(@$A) / @$A;
my $x̄_b = sum(@$B) / @$B;
my $a_var = sum(map { ($x̄_a - $_)**2 } @$A) / (@$A - 1);
my $b_var = sum(map { ($x̄_b - $_)**2 } @$B) / (@$B - 1);
($a_var && $b_var) || return 1;
my $Welsh_𝒕_statistic = ($x̄_a - $x̄_b) / sqrt($a_var/@$A + $b_var/@$B);
my $DoF = ($a_var/@$A + $b_var/@$B)**2 / (
$a_var**2 / (@$A**3 - @$A**2) +
$b_var**2 / (@$B**3 - @$B**2));
my $sa = $DoF / 2 - 1;
my $x = $DoF / ($Welsh_𝒕_statistic**2 + $DoF);
my $N = 65355;
my $h = $x / $N;
my ($sum1, $sum2) = (0, 0);
foreach my $k (0 .. $N - 1) {
my $i = $h * $k;
$sum1 += ($i + $h/2)**$sa / sqrt(1 - ($i + $h/2));
$sum2 += $i**$sa / sqrt(1-$i);
}
($h/6 * ($x**$sa / sqrt(1-$x) + 4*$sum1 + 2*$sum2) /
(gamma($sa + 1) * sqrt(pi) / gamma($sa + 1.5)))->numify;
}
my @tests = (
[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],
[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],
[17.2, 20.9, 22.6, 18.1, 21.7, 21.4, 23.5, 24.2, 14.7, 21.8],
[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],
[19.8, 20.4, 19.6, 17.8, 18.5, 18.9, 18.3, 18.9, 19.5, 22.0],
[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],
[30.02, 29.99, 30.11, 29.97, 30.01, 29.99],
[29.89, 29.93, 29.72, 29.98, 30.02, 29.98],
[3.0, 4.0, 1.0, 2.1],
[490.2, 340.0, 433.9],
);
while (@tests) {
my ($left, $right) = splice(@tests, 0, 2);
print p_value($left, $right), "\n";
}</lang>
{{out}}
<pre>
0.0213780014628667
0.148841696605327
0.0359722710297968
0.0907733242856612
0.0107515340333929
</pre>
| |||