Revert "simplify some calculations" as causes numerical problems, fix variance to >= 0 - fix #34

Author: mohawk2

Changes

@@ -1,3 +1,5 @@
+- fix recent change that allowed negative variance and standard deviation (#34) - thanks @shawnlaffan for report
+
 0.854 2025-02-25
 - fix anova_rptd problem depending on grouping of IVs (part of RT#97925)
 - do all plotting with PDL::Graphics::Simple

lib/PDL/Stats/Basic.pd

@@ -59,7 +59,9 @@ pp_def('stdv',
       a2 += $a() * $a();
     %}
     if (N < 1) { $SETBAD(b()); continue; }
-    $b() = sqrt( (a2 - sa*sa/N) / N );
+    $GENERIC() var = a2 / N - (sa/N)*(sa/N);
+    if (var < 0) var = 0;
+    $b() = sqrt(var);
   ',
   Doc      => 'Sample standard deviation.',
 );
@@ -77,7 +79,9 @@ pp_def('stdv_unbiased',
       a2 += $a() * $a();
     %}
     if (N < 2) { $SETBAD(b()); continue; }
-    $b() = pow( (a2 - sa*sa/N) / (N-1), .5 );
+    $GENERIC() var = a2/(N-1) - sa*sa/(N*(N-1));
+    if (var < 0) var = 0;
+    $b() = sqrt(var);
   ',
   Doc      => 'Unbiased estimate of population standard deviation.',
 );
@@ -95,7 +99,7 @@ pp_def('var',
       a2 += $a() * $a();
     %}
     if (N < 1) { $SETBAD(b()); continue; }
-    $b() = (a2 - sa*sa/N) / N;
+    $b() = a2 / N - sa*sa/(N*N);
   ',
   Doc      => 'Sample variance.',
 );
@@ -130,8 +134,9 @@ pp_def('se',
       sa += $a();
       a2 += $a() * $a();
     %}
-    if (N < 2) { $SETBAD(b()); continue; }
-    $b() = sqrt( (a2 - sa*sa/N) / (N*(N-1)) );
+    $GENERIC() se2 = (a2 - sa*sa/N) / (N*(N-1));
+    if (se2 < 0) se2 = 0;
+    $b() = sqrt(se2);
   ',
   Doc      => '
 =for ref
@@ -275,7 +280,7 @@ pp_def('cov',
       sb += $b();
     %}
     if (N < 1) { $SETBAD(c()); continue; }
-    $c() = (ab - sa*sb/N) / N;
+    $c() = ab / N - (sa/N) * (sb/N);
   ',
   Doc      => 'Sample covariance. see B<corr> for ways to call',
 );

t/basic.t

@@ -18,6 +18,14 @@ is_pdl $a->skew_unbiased, pdl( 0 ), "unbiased sample skewness of $a";
 is_pdl $a->kurt, pdl( -1.3 ), "sample kurtosis of $a";
 is_pdl $a->kurt_unbiased, pdl( -1.2 ), "unbiased sample kurtosis of $a";
 
+{
+my $x = pdl [(0.001) x 6];
+my $var = $x->var;
+ok $var >= 0, 'var >= 0' or diag "var = $var";
+my $stdv = $x->stdv;
+ok $stdv >= 0, 'stdv >= 0' or diag "stdv = $stdv";
+}
+
 is_pdl $_->ss, (($_ - $_->avg)**2)->sumover, "ss for $_" for
   pdl('[1 1 1 1 2 3 4 4 4 4 4 4]'),
   pdl('[1 2 2 2 3 3 3 3 4 4 5 5]'),