|
</pre>
=== Using BurkardtBurkhardt's betain'incomplete beta' ===
{{works with|Rakudo|20182019.1011}}
{{trans|Perl}}
}
sub pvaluep-value (@a, @b) {
return 1 if @a.elems <=| 1@b.elems {≤ 1;
my Rat my $mean1 = @a.sum / @a.elems; ▼
return 1.0;
my Rat my $mean2 = @b.sum / @b.elems; ▼
return 1 if $mean1 == $mean2 {;▼
if @b.elems <= 1 {
return 1.0;
}
▲ my Rat $mean1 = @a.sum / @a.elems;
▲ my Rat $mean2 = @b.sum / @b.elems;
return 1.0;
}
my Rat $variance1 = 0.0; ▼ my Rat $variance2 = 0.0;
my $variance1 /= sum (@a .elems «- » 1$mean1) X**2; ▼
for @a -> $i {
my $variance2 += sum ( $mean2@b «- » $ imean2) X**2 ;▼
$variance1 += ($mean1 - $i)**2#";" unnecessary for last statement in block
return 1 if ($variance1 == 0 &&| $variance2 == 0 ) {;▼
}
for @b -> $i {
▲ $variance2 += ($mean2 - $i)**2
}
▲ if ($variance1 == 0 && $variance2 == 0) {
return 1.0;
}
▲ $variance1 /= (@a.elems - 1);
$variance2 /= (@b.elems - 1);
▲ my Rat $variance1 /= 0@a. 0elems - 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) - 1;
▲ my $ WELCH_T_STATISTICWelchs-𝒕-statistic = ($mean1-$mean2)/sqrt($variance1/@a.elems+$variance2/@b.elems);
/
my $DoF = ($variance1/@a.elems + $variance2/@b.elems)² /
(
(($variance1* × $variance1)/(@a.elems* × @a.elems* × (@a.elems-1)) +
($variance2* × $variance2)/(@b.elems* × @b.elems* × (@b.elems-1))
);
);
my $A = $DEGREES_OF_FREEDOMDoF / 2;
my $value = $DEGREES_OF_FREEDOMDoF / ($WELCH_T_STATISTIC*$WELCH_T_STATISTICWelchs-𝒕-statistic² + $DEGREES_OF_FREEDOMDoF);
return $value if $A | $value ==≤ 0.0 ||or $value ==≥ 1 .0; ▼ my $beta = lgamma($A)+0.57236494292470009-lgamma($A+0.5); ▼
my Rat $acu = 10**(-15);
# from here, translation of John Burkhardt's C
my ($ai,$cx,$indx,$ns,$pp,$psq,$qq,$rx,$temp,$term,$xx);
▲ my $beta = lgamma($A) + 0.57236494292470009 - lgamma($A+0.5); # constant is logΓ(.5), more precise than 'lgamma' routine
# Check the input arguments.
return $value if $A <=my 0.0;#$eps || $q <= 0.010**-15;
return my $valuepsq if= $valueA <+ 0.0 || 1.0 < $value5;
# Special cases
my ($xx,$pp,$qq,$indx);
▲ return $value if $value == 0.0 || $value == 1.0;
if $A < $psq × $value { ($xx , $cx, $pp, $qq, $indx) = ($ valuecx, $ Avalue, 0.5, 0); $A, 1 }▼ $psq = $A + 0.5;
else { ($xx, $pp, $qq, $indx) = $value, $A, 0.5, 0 }
if my $Aterm <= my $psqai *= $value {= 1;
($xx, $cx, $pp, $qq, my $indx)ns = (floor $cx,qq + $value,cx 0.5,× $A, 1)psq;
} else {
# Soper reduction formula .▼
▲ ($xx, $pp, $qq, $indx) = ($value, $A, 0.5, 0);
my $ tempqq-ai = $qq - $ai; ▼ }
my $rx = $ns == 0 ?? $xx !! $xx / $cx;
$term = 1.0;
loop {
$ai = 1.0;
$term ×= $qq-ai × $rx / ($pp + $ai);
$value = 1.0;
$ns = $qq + $cxvalue *+= $psqterm;
$nsqq-ai = $nsterm.Intabs;
if $qq-ai ≤ $eps & $eps×$value {
▲#Soper reduction formula.
$value = $value × ($pp × $xx.log + ($qq - 1) × $cx.log - $beta).exp / $pp;
$rx = $xx / $cx;
$tempvalue = $qq1 - $aivalue if $indx;
last
}
while (True) {
$term = $term * $temp * $rx / ( $pp ai++ $ai );
$value = $value + $termns--;
▲ $rx = $xx if $ns ==≥ 0 ; {
$temp = $term.abs;
if $temp <= $acu && $tempqq-ai <= $acuqq *- $value {ai;
$value = $value * ($pp * $xx.log + ( $qqrx - 1.0) * = $cx.logxx -if $beta).expns /== $pp0;
$value = 1.0 - $value if} $indx;else {
$qq-ai = $psq;
last;
$psq += 1;
}
$ai = $ai + 1.0; }
$ns--;
if 0 <= $ns {
$rx = $xx if $ns == 0;
} else {
$temp = $psq;
$psq = $psq + 1.0;
}
}
}
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>],
[< 0.010268 0.000167 0.000167>],
[< 0.159258 0.136278 0.122389>],
[< 1.0/15 10.0/62.0>],
[< 1.0/10 2/50.0>],
[< 9/23.0 21/45.0 0/38.0>],
[< 0/44.0 42/94.0 0/22.0>],
) -> @left, @right {
my $pvaluep-value = pvaluep-value @left, @right;
printf("p-value = %.14g\n",$pvalue);
$error += abs($pvaluep-value - shift @answers);
}
printf("cumulative error is %g\n", $error);</lang>
|