From 5634a0c4cb77bbbb8475cf4020133c9f91e17604 Mon Sep 17 00:00:00 2001
From: "E. Joshua Rigler" <erigler@usgs.gov>
Date: Thu, 20 Oct 2016 10:12:27 -0600
Subject: [PATCH] Added option to smooth SQ locally
By default, each element of the SQ vector (i.e., shat returned from additive())
is treated as if it is independent of its nearest neighbors in time. This is
not generally true, but to treat this in a completely statistically robust
manner would be computationally prohibitive. Instead, we apply a Gaussian
smoother/filter to a limited set of SQ elements in the neighborhood of the
current time step. The configuration parameter 'smooth' specifies the period
(in samples) at which one requires the amplitude to be attenuated by half.
Longer periods will be less attenuated, while shorter periods will be more
attenuated.
Also, added the running estimate of the standard error (i.e., sigmahat) to
the stream returned by process_one().
---
geomagio/algorithm/SqDistAlgorithm.py | 49 +++++++++++++++++++++------
1 file changed, 38 insertions(+), 11 deletions(-)
diff --git a/geomagio/algorithm/SqDistAlgorithm.py b/geomagio/algorithm/SqDistAlgorithm.py
index 1141a1fdc..5ebca90db 100644
--- a/geomagio/algorithm/SqDistAlgorithm.py
+++ b/geomagio/algorithm/SqDistAlgorithm.py
@@ -230,12 +230,13 @@ class SqDistAlgorithm(Algorithm):
out += self.create_trace(channel + '_Dist', trace.stats, dist)
out += self.create_trace(channel + '_SQ', trace.stats, sq)
out += self.create_trace(channel + '_SV', trace.stats, sv)
+ out += self.create_trace(channel + '_Sigma', trace.stats, sigmahat)
return out
@classmethod
def additive(cls, yobs, m, alpha, beta, gamma, phi=1,
yhat0=None, s0=None, l0=None, b0=None, sigma0=None,
- zthresh=6, fc=0, hstep=0):
+ zthresh=6, fc=0, hstep=0, smooth=1):
"""Primary function for Holt-Winters smoothing/forecasting with
damped linear trend and additive seasonal component.
@@ -291,6 +292,9 @@ class SqDistAlgorithm(Algorithm):
hstep : int
the number of steps ahead to predict yhat[i]
which forces an hstep prediction at each time step
+ smooth: int
+ period (in samples) at which Gaussian smoother will attenuate
+ signal power by half
Returns
-------
@@ -358,6 +362,27 @@ class SqDistAlgorithm(Algorithm):
raise AlgorithmException(
"sigma0 must have length %d" % (hstep + 1))
+
+ # generate a vector of weights that will "smooth" seasonal
+ # variations locally...the quotes are because what we really
+ # do is distribute the error correction across a range of
+ # seasonal corrections; we do NOT convovle this filter with
+ # a signal to dampen noise, as is typical. -EJR 10/2016
+
+ # smooth parameter should specify the required cut-off period in terms
+ # of discrete samples; for now, we generate a Gaussian filter according
+ # to White et al. (USGS SIR 2014-5045).
+ fom = 10**(-3/20.) # halve power at corner frequency
+ omg = np.pi / np.float64(smooth) # corner angular frequency
+ sig = np.sqrt(-2*np.log(fom)/omg**2) + np.finfo(float).eps # prevent sig==0
+ ts = np.linspace(np.max((-m, -3*np.round(sig))),
+ np.min((m, 3*np.round(sig))),
+ np.min((2*m, 6*np.round(sig))) + 1) # this is alwyays symmetric
+ nts = ts.size
+ weights = np.exp(-0.5 * (ts / sig)**2)
+ weights = weights / np.sum(weights)
+
+
#
# Now begin the actual Holt-Winters algorithm
#
@@ -436,17 +461,15 @@ class SqDistAlgorithm(Algorithm):
# performance.
r[i + 1] = gamma * (1 - alpha) * et / m
- ## this properly smooths SQ locally; figure out if we need an
- ## input parameter to control this, or if we just fix it to a
- ## reasonable, but static value; also, figure out how/if r
- ## needs to be changed -EJR 8/2016
- #s[i + m] = s[i] + gamma * (1 - alpha) * et / 61.
- #s[i+m-30:i+m] = s[i+m-30:i+m] + gamma * (1-alpha) * et / 61.
- #s[i+1:i+30+1] = s[i+1:i+30+1] + gamma * (1-alpha) * et / 61.
-
- # update and append to s using equation-error formulation
- s[i + m] = s[i] + gamma * (1 - alpha) * et
+ ## update and append to s using equation-error formulation
+ #s[i + m] = s[i] + gamma * (1 - alpha) * et
+ # distribute error correction across range of seasonal
+ # corrections according to weights calculated above
+ s[i + m] = s[i] + gamma * (1-alpha) * et * weights[nts/2]
+ s[i+m-nts/2:i+m] = s[i+m-nts/2:i+m] + gamma * (1-alpha) * et * weights[:nts/2]
+ s[i+1:i+nts/2+1] = s[i+1:i+nts/2+1] + gamma * (1-alpha) * et * weights[nts/2+1:]
+
# update l and b using equation-error formulation
l = l + phi * b + alpha * et
b = phi * b + alpha * beta * et
@@ -617,6 +640,9 @@ class SqDistAlgorithm(Algorithm):
parser.add_argument('--sqdist-zthresh',
default=6,
help='Set Z-score threshold')
+ parser.add_argument('--sqdist-smooth',
+ default=1,
+ help='Local SQ smoothing parameter')
def configure(self, arguments):
"""Configure algorithm using comand line arguments.
@@ -634,4 +660,5 @@ class SqDistAlgorithm(Algorithm):
self.mag = arguments.sqdist_mag
self.statefile = arguments.sqdist_statefile
self.zthresh = arguments.sqdist_zthresh
+ self.smooth = arguments.sqdist_smooth
self.load_state()
--
GitLab