From d1b24f46135bd7fd0ce63b0f1835785541e9264c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Offringa?= <offringa@astron.nl>
Date: Sat, 25 Feb 2012 13:54:30 +0000
Subject: [PATCH] Task #1892: error values on slope
---
.../include/AOFlagger/quality/loghistogram.h | 28 +++++++++++++++++++
.../src/gui/quality/histogrampage.cpp | 12 ++++----
2 files changed, 35 insertions(+), 5 deletions(-)
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h b/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h
index 4f335323d7a..6b10279e8e3 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h
@@ -263,6 +263,34 @@ class LogHistogram : public Serializable
return (double) (sumOffset/(long double) n);
}
+ double NormalizedSlopeStdDev(double startAmplitude, double endAmplitude, double slope, double offset) const
+ {
+ long double sqErrorSum = 0.0, xsqErrorSum = 0.0, xSum = 0.0;
+ unsigned long n = 0;
+ // determine the 'average' x
+ for(const_iterator i=begin();i!=end();++i)
+ {
+ if(i.value() >= startAmplitude && i.value() < endAmplitude)
+ {
+ xSum += log10(i.value());
+ ++n;
+ }
+ }
+ const long double avgX = xSum / n;
+ for(const_iterator i=begin();i!=end();++i)
+ {
+ if(i.value() >= startAmplitude && i.value() < endAmplitude)
+ {
+ long double y = log10(i.normalizedCount());
+ long double x = log10(i.value());
+ long double ySlope = x*slope + offset;
+ sqErrorSum += (y-ySlope)*(y-ySlope);
+ xsqErrorSum += (x - avgX)*(x - avgX);
+ }
+ }
+ return (double) sqrtl(sqErrorSum/(xsqErrorSum * (long double) (n-2)));
+ }
+
double PowerLawUpperLimit(double constrainingAmplitude, double exponent, double factor) const
{
const double count = NormalizedCountAbove(constrainingAmplitude);
diff --git a/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp b/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp
index 3053bae58f9..0f7115f5f75 100644
--- a/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp
+++ b/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp
@@ -453,11 +453,13 @@ void HistogramPage::addSlopeText(std::stringstream &str, const LogHistogram &his
}
double rfiRatio = atof(_slopeRFIRatio.get_text().c_str());
- double slope = histogram.NormalizedSlope(minRange, maxRange);
- double offset = histogram.NormalizedSlopeOffset(minRange, maxRange, slope);
- double upperLimit = histogram.PowerLawUpperLimit(minRange, slope, pow10(offset));
- double lowerLimit = histogram.PowerLawLowerLimit(minRange, slope, pow10(offset), rfiRatio);
- str << '\n' << slope << '[' << log10(lowerLimit) << ';' << log10(upperLimit) << ']';
+ const double
+ slope = histogram.NormalizedSlope(minRange, maxRange),
+ offset = histogram.NormalizedSlopeOffset(minRange, maxRange, slope),
+ error = histogram.NormalizedSlopeStdDev(minRange, maxRange, slope, offset),
+ upperLimit = histogram.PowerLawUpperLimit(minRange, slope, pow10(offset)),
+ lowerLimit = histogram.PowerLawLowerLimit(minRange, slope, pow10(offset), rfiRatio);
+ str << '\n' << slope << "±" << error << '[' << log10(lowerLimit) << ';' << log10(upperLimit) << ']';
}
void HistogramPage::updateDataWindow()
--
GitLab