Welch's t-test: Difference between revisions
Content added Content deleted
Thundergnat (talk | contribs) m (Fix erroneous header which caused page to be not indexed) |
Thundergnat (talk | contribs) (→{{header|Perl 6}}: Add a Perl 6 example) |
||
Line 362: | Line 362: | ||
d9 p2_tail da |
d9 p2_tail da |
||
0.0107377</lang> |
0.0107377</lang> |
||
=={{header|Perl 6}}== |
|||
{{works with|Rakudo|2017.08}} |
|||
{{trans|C}} |
|||
Perhaps "inspired by C example" may be more accurate. Gamma subroutine from [[Gamma_function#Perl_6 |Gamma function task]]. |
|||
<lang perl6>sub Γ(\z) { |
|||
constant g = 9; |
|||
z < .5 ?? π / sin(π * z) / Γ(1 - z) !! |
|||
τ.sqrt * (z + g - 1/2)**(z - 1/2) * |
|||
exp(-(z + g - 1/2)) * |
|||
[+] < |
|||
1.000000000000000174663 |
|||
5716.400188274341379136 |
|||
-14815.30426768413909044 |
|||
14291.49277657478554025 |
|||
-6348.160217641458813289 |
|||
1301.608286058321874105 |
|||
-108.1767053514369634679 |
|||
2.605696505611755827729 |
|||
-0.7423452510201416151527e-2 |
|||
0.5384136432509564062961e-7 |
|||
-0.4023533141268236372067e-8 |
|||
> Z* 1, |map 1/(z + *), 0..* |
|||
} |
|||
sub p-value (@A, @B) { |
|||
return 1 if @A <= 1 or @B <= 1; |
|||
my $a-mean = @A.sum/@A.elems; |
|||
my $b-mean = @B.sum/@B.elems; |
|||
my $a-variance = @A.map( { ($a-mean - $_)² } ).sum / (@A - 1); |
|||
my $b-variance = @B.map( { ($b-mean - $_)² } ).sum / (@B - 1); |
|||
return 1 unless $a-variance && $b-variance; |
|||
my \Welsh-𝒕-statistic = ($a-mean - $b-mean)/($a-variance/+@A + $b-variance/+@B).sqrt; |
|||
my $DoF = ($a-variance/+@A + $b-variance/+@B)² / |
|||
(($a-variance² / (+@A² * (+@A - 1))) + |
|||
($b-variance² / (+@B² * (+@B - 1)))); |
|||
my $sa = $DoF / 2 - 1; |
|||
my $x = $DoF / (Welsh-𝒕-statistic² + $DoF); |
|||
my $N = 65355; |
|||
my $h = $x / $N; |
|||
my ( $sum1, $sum2 ); |
|||
for ^$N »*» $h -> $i { |
|||
$sum1 += (($i + $h / 2) ** $sa ) / ((1 - ($i + $h / 2)).sqrt); |
|||
$sum2 += $i ** $sa / (1 - $i).sqrt; |
|||
} |
|||
(($h / 6) * ( $x ** $sa / (1 - $x).sqrt + 4 * $sum1 + 2 * $sum2)) / |
|||
( Γ($sa + 1) * 1.77245385090551610 / Γ($sa + 1.5) ); |
|||
} |
|||
# Testing |
|||
for ( |
|||
[<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>] |
|||
) -> @left, @right { say p-value @left, @right }</lang> |
|||
{{out}} |
|||
<pre>0.0213780014628669 |
|||
0.148841696605328 |
|||
0.0359722710297969 |
|||
0.0907733242856673 |
|||
0.010751534033393 |
|||
</pre> |
|||
=={{header|R}}== |
=={{header|R}}== |