Changeset 17602

Timestamp:

May 9, 2008, 11:02:21 AM (18 years ago)

Author:

eugene

Message:

added clipmean and clipstdev options

File:

• trunk/PS-IPP-Config/lib/PS/IPP/Metadata/Stats.pm (modified) (4 diffs)

trunk/PS-IPP-Config/lib/PS/IPP/Metadata/Stats.pm

@@ -1 +1 @@
 # Copyright (c) 2006  Paul Price, Joshua Hoblitt
 #
-# $Id: Stats.pm,v 1.27 2008-04-11 07:37:14 eugene Exp $
+# $Id: Stats.pm,v 1.28 2008-05-09 21:02:21 eugene Exp $
 
 package PS::IPP::Metadata::Stats;
@@ -109 +109 @@
         }
 
+        # calculate a 3-sigma clipped mean (if Nvalues > 8, reject > 3sigma outliers)
+        if ($type eq "clipmean") {
+            if (scalar @$data > 0) {
+                # Get statistics on the value
+                my $stats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
+                $stats->add_data(@$data);
+                my $rawMean = $stats->mean();
+                my $rawStdev = $stats->standard_deviation();
+                my $rawLimit = 3.0*$rawStdev;
+
+                # clip, if we can
+                if (scalar @$data > 8) {
+                    my @clipdata = ();
+                    foreach my $value (@$data){ 
+                        my $delta = abs($value - $rawMean);
+                        if ($delta > $rawLimit) { next; }
+                        push @clipdata, $value;
+                    }
+                    my $clipstats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
+                    $clipstats->add_data(@clipdata);
+                    $entry->{value} = $clipstats->mean();
+                } else {
+                    $entry->{value} = $stats->mean();
+                }
+            } else {
+                $self->_null_for_type ($entry);
+            }
+            next;
+        }
+
         if ($type eq "stdev") {
             if (scalar @$data > 1) {
@@ -122 +152 @@
             } 
             $self->_null_for_type ($entry);
+            next;
+        }
+
+        # calculate a 3-sigma clipped stdev (if Nvalues > 8, reject > 3sigma outliers)
+        if ($type eq "clipstdev") {
+            if (scalar @$data > 0) {
+                # Get statistics on the value
+                my $stats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
+                $stats->add_data(@$data);
+                my $rawMean = $stats->mean();
+                my $rawStdev = $stats->standard_deviation();
+                my $rawLimit = 3.0*$rawStdev;
+
+                # clip, if we can
+                if (scalar @$data > 8) {
+                    my @clipdata = ();
+                    foreach my $value (@$data){ 
+                        my $delta = abs($value - $rawMean);
+                        if ($delta > $rawLimit) { next; }
+                        push @clipdata, $value;
+                    }
+                    my $clipstats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
+                    $clipstats->add_data(@clipdata);
+                    $entry->{value} = $clipstats->standard_deviation();
+                } else {
+                    $entry->{value} = $stats->standard_deviation();
+                }
+            } else {
+                $self->_null_for_type ($entry);
+            }
             next;
         }
@@ -270 +330 @@
             }
 
+            if ($type eq 'clipmean') {
+                push @$data, $value;
+                next;
+            }
+
             if ($type eq 'stdev') {
+                push @$data, $value;
+                next;
+            }
+
+            if ($type eq 'clipstdev') {
                 push @$data, $value;
                 next;