Talk:Welch's t-test: Difference between revisions
→cannot run python libraries: move to Hailholyghost's talk page (the problem was unrelated to this task)
(→cannot run python libraries: move to Hailholyghost's talk page (the problem was unrelated to this task)) |
|||
| (43 intermediate revisions by 3 users not shown) | |||
Line 458:
:: Sounds like Math::MPFR needs to get cpan to give better diagnostics for that kind of failure? (But why not show both versions? Are they too big or something?) --[[User:Rdm|Rdm]] ([[User talk:Rdm|talk]]) 02:01, 15 November 2017 (UTC)
I have a more portable version of the current perl code, using the new integration by reduction method, and verifying against the R code as I did with the C code:
<lang perl>
#!/usr/bin/env perl
use strict; use warnings; use Cwd 'getcwd';# use feature 'say';
my $TOP_DIRECTORY = getcwd();
local $SIG{__WARN__} = sub {#kill the program if there are any warnings
my $message = shift;
my $fail_filename = "$TOP_DIRECTORY/$0.FAIL";
open my $fh, '>', $fail_filename or die "Can't write $fail_filename: $!";
printf $fh ("$message @ %s\n", getcwd());
close $fh;
die "$message\n";
};#http:#perlmaven.com/how-to-capture-and-save-warnings-in-perl
sub betain {
my $x = shift;
my $p = shift;
my $q = shift;
my $beta = shift;
#modified from John Burkardt
my $acu = 10**(-15);
my $ai;
my $cx;
my $indx;#int
my $ns;#int
my $pp;
my $psq;
my $qq;
my $rx;
my $temp;
my $term;
my $value;
my $xx;
$value = $x;
# ifault = 0;
#Check the input arguments.
if ( ($p <= 0.0) || ($q <= 0.0 )){
# *ifault = 1;
return $value;
}
if ( $x < 0.0 || 1.0 < $x )
{
# *ifault = 2;
return $value;
}
#/*
# Special cases.
#*/
if ( $x == 0.0 || $x == 1.0 ) {
return $value;
}
$psq = $p + $q;
$cx = 1.0 - $x;
if ( $p < $psq * $x )
{
$xx = $cx;
$cx = $x;
$pp = $q;
$qq = $p;
$indx = 1;
}
else
{
$xx = $x;
$pp = $p;
$qq = $q;
$indx = 0;
}
$term = 1.0;
$ai = 1.0;
$value = 1.0;
$ns = sprintf("%d",( $qq + $cx * $psq ));#this evaluates differently, because $ns is a C integer
#
# Use the Soper reduction formula.
#
$rx = $xx / $cx;
$temp = $qq - $ai;
if ( $ns == 0 )
{
$rx = $xx;
}
for ( ; ; )
{
$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;
if ( $indx )
{
$value = 1.0 - $value;
}
last;
}
$ai = $ai + 1.0;
$ns = $ns - 1;
if ( 0 <= $ns )
{
$temp = $qq - $ai;
if ( $ns == 0 )
{
$rx = $xx;
}
}
else
{
$temp = $psq;
$psq = $psq + 1.0;
}
}
return $value;
}
#lgamma subroutine is less precise than "use POSIX 'lgamma'" but Macs don't have POSIX
sub lgamma {#http://hea-www.harvard.edu/~alexey/calc-src.txt
my $xx = shift;
my ($j, $ser, $tmp, $x, $y);
my @cof = (0.0, 76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155, 0.1208650973866179e-2,
-0.5395239384953e-5);
my $stp = 2.5066282746310005;
$x = $xx; $y = $x;
$tmp = $x + 5.5;
$tmp = ($x+0.5)*log($tmp)-$tmp;
$ser = 1.000000000190015;
foreach $j ( 1 .. 6 ) {
$y+=1.0;
$ser+=$cof[$j]/$y;
}
return $tmp + log($stp*$ser/$x);
}
use List::Util 'sum';
sub calculate_Pvalue {
my $array1 = shift;
my $array2 = shift;
if (scalar @$array1 <= 1) {
return 1.0;
}
if (scalar @$array2 <= 1) {
return 1.0;
}
my $mean1 = sum(@{ $array1 });
$mean1 /= scalar @$array1;
my $mean2 = sum(@{ $array2 });
$mean2 /= scalar @$array2;
if ($mean1 == $mean2) {
return 1.0;
}
# say "mean1 = $mean1 mean2 = $mean2";
my $variance1 = 0.0;
my $variance2 = 0.0;
foreach my $x (@$array1) {
$variance1 += ($x-$mean1)*($x-$mean1);
}
foreach my $x (@$array2) {
$variance2 += ($x-$mean2)*($x-$mean2);
}
if (($variance1 == 0.0) && ($variance2 == 0.0)) {
return 1.0;
}
$variance1 = $variance1/(scalar @$array1-1);
$variance2 = $variance2/(scalar @$array2-1);
my $array1_size = scalar @$array1;
my $array2_size = scalar @$array2;
my $WELCH_T_STATISTIC = ($mean1-$mean2)/sqrt($variance1/$array1_size+$variance2/$array2_size);
my $DEGREES_OF_FREEDOM = (($variance1/$array1_size+$variance2/(scalar @$array2))**2)
/
(
($variance1*$variance1)/($array1_size*$array1_size*($array1_size-1))+
($variance2*$variance2)/($array2_size*$array2_size*($array2_size-1))
);
# say "$WELCH_T_STATISTIC $DEGREES_OF_FREEDOM";
my $A = $DEGREES_OF_FREEDOM/2;
my $x = $DEGREES_OF_FREEDOM/($WELCH_T_STATISTIC*$WELCH_T_STATISTIC+$DEGREES_OF_FREEDOM);
return betain($x, $A, 0.5, lgamma($A)+0.57236494292470009-lgamma($A+0.5));
}
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 = calculate_Pvalue(\@d1,\@d2);
my $error = abs($pvalue - $CORRECT_ANSWERS[0]);
printf("Test sets 1 p-value = %g\n",$pvalue);
$pvalue = calculate_Pvalue(\@d3,\@d4);
$error += abs($pvalue - $CORRECT_ANSWERS[1]);
printf("Test sets 2 p-value = %g\n",$pvalue);
$pvalue = calculate_Pvalue(\@d5,\@d6);
$error += abs($pvalue - $CORRECT_ANSWERS[2]);
printf("Test sets 3 p-value = %g\n",$pvalue);
$pvalue = calculate_Pvalue(\@d7,\@d8);
$error += abs($pvalue - $CORRECT_ANSWERS[3]);
printf("Test sets 4 p-value = %g\n",$pvalue);
$pvalue = calculate_Pvalue(\@x,\@y);
$error += abs($pvalue - $CORRECT_ANSWERS[4]);
printf("Test sets 5 p-value = %g\n",$pvalue);
$pvalue = calculate_Pvalue(\@v1,\@v2);
$error += abs($pvalue - $CORRECT_ANSWERS[5]);
printf("Test sets 6 p-value = %g\n",$pvalue);
$pvalue = calculate_Pvalue(\@s1,\@s2);
$error += abs($pvalue - $CORRECT_ANSWERS[6]);
printf("Test sets 7 p-value = %g\n",$pvalue);
$pvalue = calculate_Pvalue(\@z1,\@z2);
$error += abs($pvalue - $CORRECT_ANSWERS[7]);
printf("Test sets z p-value = %g\n",$pvalue);
printf("the cumulative error is %g\n", $error);</lang>
I think this should replace the current perl code (because frankly, I can't get all the modules surrounding the current one to work). I realize that Macs don't have POSIX, so I replaced the 'lgamma' with something more portable. It should work anywhere with Perl 5.0 and up.
Let me know what you think.--[[User:Hailholyghost|Hailholyghost]] ([[User talk:Hailholyghost|talk]]) 01:49, 5 January 2018 (UTC)
: I like the existing code. Try the following simple changes to the current code:
<pre>
3c3
< use Math::AnyNum qw(gamma pi);
---
> use Math::Cephes qw(gamma $PI);
38c38
< (gamma($sa + 1) * sqrt(pi) / gamma($sa + 1.5)))->numify;
---
> (gamma($sa + 1) * sqrt($PI) / gamma($sa + 1.5)));
</pre>
It gives similar results and the dependency chain should be much easier. There are a lot of machines that don't have MPFR installed, and some that don't even have GMP. The downside is that Math::AnyNum is a much more functional module that builds on others. But Math::Cephes is pretty cool as a standalone native precision library. I've thought about adding gamma to the ntheory module but it's complicated (for me at least). [[User:Danaj|Danaj]] ([[User talk:Danaj|talk]]) 22:31, 5 January 2018 (UTC)
I tried the existing module again, and it still won't even install.
The entire purpose of creating this Rosetta Code page was because I had numerous issues installing libraries in the first place, and I kept getting responses like "but the libraries do work!" despite evidence to the contrary. If a task solution is not portable, it shouldn't be shown. Also, the Perl translation I've done works very differently than the current one, more accurately, as I show. Unlike the current code, I can actually run my code, and it should run anywhere.
What is the Rosetta code policy for using an "alternative solution" on Rosetta Code?
: It's unclear with your message above as to whether you tried installing Math::AnyNum again or are replying to my comment about using Math::Cephes. Without the MPFR library installed in your O/S, Math::MPFR can't install, and therefore Math::AnyNum can't install. The libraries ''do'' work for those who have the requisite software installed. Your argument about portability turns into "don't ever use modules outside of core" which I don't believe is the way many people use Perl.
: This is mostly moot as I've added a translation of the C code which uses only List::Util's sum, based in part on your translation above. Hopefully this will resolve your issues with the page.
: I don't know what sort of official policies there are, but my take for Perl, which is intimately tied with CPAN, is that I like to see a version that uses no modules or only core modules (e.g. List::Util), followed by, if appropriate, versions using modules. The first version is useful in that it can run everywhere with no module installation, and it shows all the work. The versions using modules are one or more of faster, shorter, clearer, more functional, or more idiomatic. Especially if the no-module version is very long, I think it's best to reverse the order -- show the 2-10 line easy solutions, then the long "if you want to write everything by hand:" method.
: It is also not unusual, especially as new tasks get more complex, to see the first code added to a page using a module. This lets us get a Perl solution up quickly and succinctly, and later someone can write a huge everything-by-hand version. I suppose there is another discussion of whether it would make sense to have a "common library" pointer so we don't have to keep pasting in the same helper functions on 10+ tasks, which would let one concentrate more on the task.
thank you for your work, this new Perl version is better than what I had done.--[[User:Hailholyghost|Hailholyghost]] ([[User talk:Hailholyghost|talk]]) 12:47, 10 January 2018 (UTC)
== New title ==
Line 844 ⟶ 1,123:
:With
<lang c>int z = 0;
printf("%g =? %g\n", return_value, betain(x, a, 0.5, lgammal(a)+0.57236494292470009-lgammal(a+0.5), &z));</lang>
:
:There is another problem: you damaged the original function, by replacing <code>*ifault = 0;</code> with <code>ifault = 0;</code>, replacing ifault with a null pointer. This will crash the program with a null pointer dereference whenever the arguments passed to betain are incorrect (hence ifault is assigned a value). I didn't check if you introduced another problem in the code.
:Your statement that ''existing libraries don't calculate correct numbers'' is wrong, at least for this one.
:[[User:Eoraptor|Eoraptor]] ([[User talk:Eoraptor|talk]]) 19:31, 28 December 2017 (UTC)
=== Stata translation of Burkardt's C translation of AS 63 ===
For the record, here is a quick Stata translation of this betain implementation.
The original C code can be simplified in several places using shortcuts (it could be simplified even more in C). Also, ifault is not used here: instead, the function returns a missing value if the arguments are invalid. And log(beta(p,q)) is computed inside betain, rather than passed as an argument.
<lang stata>mata
/*
This code is distributed under the GNU LGPL license.
Modified 28dec2017
Author:
Original Fortran 77 version by KL Majumder, GP Bhattacharjee.
C version by John Burkardt.
Stata version by Jean-Claude Arbaut.
*/
function betain(x,p,q) {
/* Check the input arguments. */
if (p<=0 | q<=0 | x<0 | x>1) return(.)
/* Special cases. */
if (x==0 | x==1) {
return(x)
}
acu = 1e-15
lnbeta = lngamma(p)+lngamma(q)-lngamma(p+q)
/* Change tail if necessary and determine S. */
psq = p+q
if (p<psq*x) {
xx = 1-x
cx = x
pp = q
qq = p
indx = 1
}
else {
xx = x
cx = 1-x
pp = p
qq = q
indx = 0
}
term = ai = value = 1
ns = floor(qq+cx*psq)
/* Use the Soper reduction formula. */
rx = xx/cx
temp = qq-ai
if (ns==0) rx = xx
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)*log(cx)-lnbeta)/pp
return(indx?1-value:value)
}
ai++
if (--ns>=0) {
temp = qq-ai
if (ns==0) rx = xx
}
else temp = psq++
}
}
end</lang>
A note on the license: I used the same license as John Burkardt, but I am not certain it's permitted. The [http://lib.stat.cmu.edu/apstat/ StatLib site] states that ''The Royal Statistical Society holds the copyright to these routines, but has given its permission for their distribution provided that no fee is charged.'' However, LGPL allows redistribution in nonfree software.
[[User:Eoraptor|Eoraptor]] ([[User talk:Eoraptor|talk]]) 21:31, 28 December 2017 (UTC)
| |||