Anonymous user
Welch's t-test: Difference between revisions
→{{header|Perl 6}}: added Soper reduction version
(→{{header|Perl 6}}: added Soper reduction version) |
|||
Line 1,211:
0.0907733242856673
0.010751534033393
</pre>
=== Using Burkardt's betain ===
{{works with|Rakudo|2018.10}}
{{trans|Perl}}
This uses the Soper reduction formula to evaluate the integral, which converges much more quickly than Simpson's formula.
<lang perl6>
sub lgamma ( Num(Real) \n --> Num ){
use NativeCall;
sub lgamma (num64 --> num64) is native {}
lgamma( n )
}
sub pvalue (@a, @b) {
if @a.elems <= 1 {
return 1.0;
}
if @b.elems <= 1 {
return 1.0;
}
my Rat $mean1 = @a.sum / @a.elems;
my Rat $mean2 = @b.sum / @b.elems;
if $mean1 == $mean2 {
return 1.0;
}
my Rat $variance1 = 0.0;
my Rat $variance2 = 0.0;
for @a -> $i {
$variance1 += ($mean1 - $i)**2#";" unnecessary for last statement in block
}
for @b -> $i {
$variance2 += ($mean2 - $i)**2
}
if ($variance1 == 0 && $variance2 == 0) {
return 1.0;
}
$variance1 /= (@a.elems - 1);
$variance2 /= (@b.elems - 1);
my $WELCH_T_STATISTIC = ($mean1-$mean2)/sqrt($variance1/@a.elems+$variance2/@b.elems);
my $DEGREES_OF_FREEDOM = (($variance1/@a.elems+$variance2/@b.elems)**2)
/
(
($variance1*$variance1)/(@a.elems*@a.elems*(@a.elems-1))+
($variance2*$variance2)/(@b.elems*@b.elems*(@b.elems-1))
);
my $A = $DEGREES_OF_FREEDOM/2;
my $value = $DEGREES_OF_FREEDOM/($WELCH_T_STATISTIC*$WELCH_T_STATISTIC+$DEGREES_OF_FREEDOM);
my $beta = lgamma($A)+0.57236494292470009-lgamma($A+0.5);
my Rat $acu = 10**(-15);
my ($ai,$cx,$indx,$ns,$pp,$psq,$qq,$rx,$temp,$term,$xx);
# Check the input arguments.
return $value if $A <= 0.0;# || $q <= 0.0;
return $value if $value < 0.0 || 1.0 < $value;
# Special cases
return $value if $value == 0.0 || $value == 1.0;
$psq = $A + 0.5;
$cx = 1.0 - $value;
if $A < $psq * $value {
($xx, $cx, $pp, $qq, $indx) = ($cx, $value, 0.5, $A, 1);
} else {
($xx, $pp, $qq, $indx) = ($value, $A, 0.5, 0);
}
$term = 1.0;
$ai = 1.0;
$value = 1.0;
$ns = $qq + $cx * $psq;
$ns = $ns.Int;
#Soper reduction formula.
$rx = $xx / $cx;
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
while (True) {
$term = $term * $temp * $rx / ( $pp + $ai );
$value = $value + $term;
$temp = $term.abs;
if $temp <= $acu && $temp <= $acu * $value {
$value = $value * ($pp * $xx.log + ($qq - 1.0) * $cx.log - $beta).exp / $pp;
$value = 1.0 - $value if $indx;
last;
}
$ai = $ai + 1.0;
$ns--;
if 0 <= $ns {
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
} else {
$temp = $psq;
$psq = $psq + 1.0;
}
}
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 @s1 = 1.0/15,10.0/62.0;
my @s2 = 1.0/10,2/50.0;
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;
my @CORRECT_ANSWERS = (0.021378001462867,
0.148841696605327,
0.0359722710297968,
0.090773324285671,
0.0107515611497845,
0.00339907162713746,
0.52726574965384,
0.545266866977794);
my $pvalue = pvalue(@d1, @d2);
my $error = abs($pvalue - @CORRECT_ANSWERS[0]);
printf("Test sets 1 p-value = %.14g\n",$pvalue);
$pvalue = pvalue(@d3, @d4);
$error += abs($pvalue - @CORRECT_ANSWERS[1]);
printf("Test sets 2 p-value = %.14g\n",$pvalue);
$pvalue = pvalue(@d5, @d6);
$error += abs($pvalue - @CORRECT_ANSWERS[2]);
printf("Test sets 3 p-value = %.14g\n",$pvalue);
$pvalue = pvalue(@d7, @d8);
$error += abs($pvalue - @CORRECT_ANSWERS[3]);
printf("Test sets 4 p-value = %.14g\n",$pvalue);
$pvalue = pvalue(@x, @y);
$error += abs($pvalue - @CORRECT_ANSWERS[4]);
printf("Test sets 5 p-value = %.14g\n",$pvalue);
$pvalue = pvalue(@v1, @v2);
$error += abs($pvalue - @CORRECT_ANSWERS[5]);
printf("Test sets 6 p-value = %.14g\n",$pvalue);
$pvalue = pvalue(@s1, @s2);
$error += abs($pvalue - @CORRECT_ANSWERS[6]);
printf("Test sets 7 p-value = %.14g\n",$pvalue);
$pvalue = pvalue(@z1, @z2);
$error += abs($pvalue - @CORRECT_ANSWERS[7]);
printf("Test sets 8 p-value = %.14g\n",$pvalue);
printf("the cumulative error is %g\n", $error);
</lang>
{{out}}
<pre>
Test sets 1 p-value = 0.021378001462867
Test sets 2 p-value = 0.14884169660533
Test sets 3 p-value = 0.035972271029797
Test sets 4 p-value = 0.090773324285667
Test sets 5 p-value = 0.010751561149784
Test sets 6 p-value = 0.0033990716271375
Test sets 7 p-value = 0.52726574965384
Test sets 8 p-value = 0.54526686697779
the cumulative error is 5.50254e-15
</pre>
| |||