Rosetta Code

Welch's t-test: Difference between revisions

Page Discussion

Welch's t-test (view source)

Revision as of 17:04, 11 January 2020

151 bytes removed ,  4 years ago
→‎Using Burkhardt's 'incomplete beta': general trimming and clean-up
m (→‎Using Burkardt's betain: fix so it actually compiles)
(→‎Using Burkhardt's 'incomplete beta': general trimming and clean-up)
Line 956:
</pre>
 
=== Using BurkardtBurkhardt's betain'incomplete beta' ===
We use a slightly more accurate lgamma than the C code. Note that Perl can be compiled with different underlying floating point representations -- double, long double, or quad double.
{{trans|C}}
<lang perl>use strict;
use warnings;
use strictList::Util 'sum';
 
sub lgamma {
Line 978:
}
 
sub calculate_P_value {
use List::Util 'sum';
my ($array1,$array2) = (shift, shift);
sub calculate_Pvalue {
my return 1 if @$array1 <= 1 or @$array2 <= shift1;
 
my $array2 = shift;
my $mean1 = sum(@$array1);
if (scalar @$array1 <= 1) {
my $mean2 = sum(@$array2);
return 1.0;
$mean1 /= scalar @$array1;
}
$mean2 /= scalar @$array2;
if (scalar @$array2 <= 1) {
return 1.0 if $mean1 == $mean2;
}
my ($variance1,$variance2);
my $mean1 = sum(@{ $array1 });
$mean1variance1 /+= scalar($mean1-$_)**2 for @$array1;
my $mean2variance2 += sum(@{$mean2-$_)**2 for @$array2 });
return 1 if $variance1 == 0 and $variance2 == 0;
$mean2 /= scalar @$array2;
if ($mean1 == $mean2) {
$variance1 = $variance1/(@$array1-1);
return 1.0;
$variance2 = $variance2/(@$array2-1);
}
my $Welch_t_statistic = ($mean1-$mean2)/sqrt($variance1/$array1_size+$variance2/$array2_size);
my $variance1 = 0.0;
my $DoF = (($variance1/@$array1+$variance2/@$array2)**2) /
my $variance2 = 0.0;
(
foreach my $x (@$array1) {
($variance1*$variance1)/(@$array1*@$array1*(@$array1-1)) +
$variance1 += ($x-$mean1)*($x-$mean1);
($variance2*$variance2)/(@$array2*@$array2*(@$array2-1))
}
);
foreach my $x (@$array2) {
my $A = $DoF / 2;
$variance2 += ($x-$mean2)*($x-$mean2);
my $value = $DoF / ($Welch_t_statistic**2 + $DoF);
}
return $value if $A <= 0 or $value <= 0 or 1 <= $value;
if (($variance1 == 0.0) && ($variance2 == 0.0)) {
 
return 1.0;
# from here, translation of John Burkhardt's C code
}
my $beta = lgamma($A) + 0.57236494292470009 - lgamma($A+0.5); constant is logΓ(.5), more precise than 'lgamma' routine
$variance1 = $variance1/(scalar @$array1-1);
my $eps = 10**-15;
$variance2 = $variance2/(scalar @$array2-1);
my($ai,$cx,$indx,$ns,$pp,$psq,$qq,$qq_ai,$rx,$term,$xx);
my $array1_size = scalar @$array1;
my $array2_size = scalar @$array2;
$psq = $A + 0.5;
my $WELCH_T_STATISTIC = ($mean1-$mean2)/sqrt($variance1/$array1_size+$variance2/$array2_size);
$cx = 1 - $value;
my $DEGREES_OF_FREEDOM = (($variance1/$array1_size+$variance2/(scalar @$array2))**2)
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 = $ai = $value = 1;
($variance1*$variance1)/($array1_size*$array1_size*($array1_size-1))+
$ns = int $qq + $cx * $psq;
($variance2*$variance2)/($array2_size*$array2_size*($array2_size-1))
 
);
# Soper reduction formula
my $A = $DEGREES_OF_FREEDOM/2;
$qq_ai = $qq - $ai;
my $value = $DEGREES_OF_FREEDOM/($WELCH_T_STATISTIC*$WELCH_T_STATISTIC+$DEGREES_OF_FREEDOM);
$rx = $ns == 0 ? $xx : $xx / $cx;
#from here, translation of John Burkhardt's C
while (1) {
my $beta = lgamma($A)+0.57236494292470009-lgamma($A+0.5);
$term = $term * $qq_ai * $rx / ( $pp + $ai );
my $acu = 10**(-15);
$value = $value + $term;
my($ai,$cx,$indx,$ns,$pp,$psq,$qq,$rx,$temp,$term,$xx);
$qq_ai = abs($term);
# Check the input arguments.
return $value if ($Aqq_ai <= 0.0;#$eps ||&& $qqq_ai <= 0.0;$eps * $value) {
return $value if= $value <* 0.0exp ||($pp * log($xx) + ($qq - 1.0) <* log($valuecx) - $beta) / $pp;
$value = 1 - $value if $indx;
# Special cases
last;
return $value if $value == 0.0 || $value == 1.0;
$psq = $A + 0.5; }
$cx = 1.0 - $valueai++;
$ns--;
if ($A < $psq * $value) {
if ($ns >= 0) {
($xx, $cx, $pp, $qq, $indx) = ($cx, $value, 0.5, $A, 1);
$qq_ai = $qq - $ai;
} else {
($xx, $pp, $qq, $indx)rx = ($value,xx if $A,ns 0.5,== 0);
} else {
}
$qq_ai = $psq;
$term = 1.0;
$psq = $psq + 1;
$ai = 1.0;
}
$value = 1.0;
}
$ns = int($qq + $cx * $psq);
$value
#Soper reduction formula.
$rx = $xx / $cx;
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
while (1) {
$term = $term * $temp * $rx / ( $pp + $ai );
$value = $value + $term;
$temp = abs ($term);
if ($temp <= $acu && $temp <= $acu * $value) {
$value = $value * exp ($pp * log($xx)
+ ($qq - 1.0) * log($cx) - $beta) / $pp;
$value = 1.0 - $value if $indx;
last;
}
$ai = $ai + 1.0;
$ns = $ns - 1;
if (0 <= $ns) {
$temp = $qq - $ai;
$rx = $xx if $ns == 0;
} else {
$temp = $psq;
$psq = $psq + 1.0;
}
}
return $value;
}
 
Line 1,106 ⟶ 1,082:
while (@tests) {
my ($left, $right) = splice(@tests, 0, 2);
my $pvalue = calculate_Pvaluecalculate_P_value($left,$right);
$error += abs($pvalue - shift @answers);
printf("p-value = %.14g\n",$pvalue);
SqrtNegInf
2,392

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/258319"