#!/usr/athena/bin/perl
# this script stolen from /mit/zeno/Public/stats.pl

while (<>) {
    s/\s+//g;
    next unless /^-?(\.\d+|\d+(\.\d*)?)$/;
    push (@X, $_);
    $n++;
    $sumX  += $_;
    $sumX2 += $_ * $_;
}

$n > 1 || die "Must have > 1 datum!\n";

sub bynum { $a <=> $b; }
@X = sort bynum @X;
push (@X, $X[$#X]);     # just so the highest percentiles return the top val

sub percentile {
    local($j, $intj, $p) = (0, 0, @_);
    $j = ($n + 1) * $p /100;
    $intj = int(--$j);
    return ( $X[$intj] + ($j - $intj) * ($X[$intj + 1] - $X[$intj]) );
}

  #print @x;
  #print "\n======\n";

$mean = $sumX / $n;
$median = ( $X[int(($n - 1)/2)] + $X[int($n/2)] ) / 2;
$range = $X[$#X] - $X[0];
$var = ( $sumX2 - ($sumX * $sumX / $n) ) / ($n - 1);
$popvar = ($n - 1) * $var / $n;
$std = sqrt($var);
$popstd = sqrt($popvar);
printf "n: %d\nmean: %.3f\nmedian: %.3f\nrange: %.3f\nvariance: %.3f  (as population: %.3f)\nstandard deviation: %.3f  (as population: %.3f)\n",
    $n, $mean, $median, $range, $var, $popvar, $std, $popstd;

print "PERCENTILES:\n";
foreach $p (5, 10, 15, 20, 25) {
    foreach $col (0, 25, 50, 75) {
	printf "%3d%%: %12.4f |", ($p + $col), &percentile($p + $col);
    }
    print "\n";
}


exit 0;