Revert "simplify some calculations" as causes numerical problems - fix #34

This reverts commit 826ec94.

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,7 @@ pp_def('stdv',
       a2 += $a() * $a();
     %}
     if (N < 1) { $SETBAD(b()); continue; }
-    $b() = sqrt( (a2 - sa*sa/N) / N );
+    $b() = sqrt( a2 / N - (sa/N)*(sa/N) );
   ',
   Doc      => 'Sample standard deviation.',
 );
@@ -77,7 +77,7 @@ pp_def('stdv_unbiased',
       a2 += $a() * $a();
     %}
     if (N < 2) { $SETBAD(b()); continue; }
-    $b() = pow( (a2 - sa*sa/N) / (N-1), .5 );
+    $b() = pow( a2/(N-1) - sa*sa/(N*(N-1)), .5 );
   ',
   Doc      => 'Unbiased estimate of population standard deviation.',
 );
@@ -95,7 +95,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.',
 );
@@ -275,7 +275,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/glm.t

@@ -11,6 +11,11 @@ is_pdl pdl('BAD 1 2 3 4')->fill_m, pdl('2.5 1 2 3 4'), "fill_m replaces bad valu
 {
 my $stdv = pdl('BAD 1 2 3 4')->fill_rand->stdv;
 ok PDL::Core::approx( $stdv, 1.01980390271856 ) || PDL::Core::approx( $stdv, 1.16619037896906 ), "fill_rand replaces bad values with random sample of good values from same variable";
+my $x = pdl [(0.001) x 6];
+my $var = $x->var;
+ok $var >= 0, 'var >= 0' or diag "var = $var";
+$stdv = $x->stdv;
+ok $stdv >= 0, 'stdv >= 0' or diag "stdv = $stdv";
 }
 
 my $a = sequence 5;