On the other hand, it's not that hard to do these by hand. Particularly in perl.

  use strict;
  use List::Util qw(sum);
  my @d    = ( 1 .. 10_000 );
  my $s    = sum @d;
  my $mean = $s/@d;
  my $var  = (sum map { ($_-$mean)**2 } @d);
  my $std  = sqrt($var/@d);

# etc...
[download]

... The more I think about it though, if you're pulling these from a database, you don't really even need to do the stddev by hand. I imagine your database of choice has a stddev() built in. The co-varience would probably have to be calculated by hand though. Maybe "select sum( (cola - avg(cola))*(colb - avg(colb))/count(cola) ) from tablename" ... or something like that.

Your method requires that all data to be in memory at once, but it's simple to refactor it so that's not the case.

my ($cnt, $sum, $squ);
while (my ($d) = $iter->()) {
   $cnt++;
   $sum += $d;
   $squ += $d * $d;
}

my $mean = $sum / $cnt;
my $var = $squ + -2*$mean*$sum + $mean*$mean*$cnt;
my $std  = sqrt($var/$cnt);
[download]

while (my ($d) = $iter->()) can be replaced with any loop, including a file reading loop or a database fetching loop.

Update: If you don't need $var anywhere else, the last two lines can be simplified to

my $std  = sqrt($squ/$cnt - $mean*$mean);
[download]

Sure. But I started by imagining piddles, and I think you'd have to have them all in memory to use that well also. I've used PDL for a grand total of 30 minutes though, so I could be wrong.

-Paul

If you prefer a simple specialised module over full-blown solutions like PDL or Statistics::R, you could use Statistics::LineFit (which is 27k pure Perl, without any external dependencies).

Here's a simple example, which computes intercept, slope and R² for one data set (100 x/y data points):

use Statistics::LineFit;

my $x = [ 1..100 ];
my $y = [ map $_*2, 11..110 ];

my $lfit = Statistics::LineFit->new();
$lfit->setData($x, $y);

printf "a=%.4f, b=%.4f,  R2=%.4f\n", $lfit->coefficients(), $lfit->rSq
+uared();


# which outputs:  a=20.0000, b=2.0000,  R2=1.0000
[download]

( Looking at the input data, it's not too surprising that the R² is exactly 1.0 )

The module also allows to compute a number of related values, like residuals, standard error, etc. - in case you should need them.

You are correct with all your assumptions. Further, you may be interested to know that I do not want to pull out unique data coordinates for each plot. I am planning to take random samples from a base table and perform the correlation analysis for each of these random sample sets.