Skip to content
Snippets Groups Projects
Commit 9f7574a3 authored by Eddie McWhirter's avatar Eddie McWhirter
Browse files

Remove commented code and unneccessary comments from SQDistAlgorithm.

parent 701f775b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
geomagio/algorithm/SQDistAlgorithm.py
+ 17
− 64
View file @ 9f7574a3
......@@ -136,7 +136,7 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
if b0 is None:
b = (np.nanmean(yobs[m:2 * m]) - np.nanmean(yobs[0:m])) / m
b = 0 if np.isnan(b) else b # replace NaN with 0
b = 0 if np.isnan(b) else b # replace NaN with 0
else:
b = b0
if not np.isscalar(b0):
......@@ -151,14 +151,13 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
if s0 is None:
s = [yobs[i] - l for i in range(m)]
s = [i if ~np.isnan(i) else 0 for i in s] # replace NaNs with 0s
s = [i if ~np.isnan(i) else 0 for i in s] # replace NaNs with 0s
else:
s = list(s0)
if len(s) != m:
raise AlgorithmException("s0 must have length %d " % m)
if sigma0 is None:
# NOTE: maybe default should be vector of zeros???
sigma = [np.sqrt(np.nanvar(yobs))] * (hstep + 1)
else:
sigma = list(sigma0)
......@@ -186,7 +185,7 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
initial_values = [alpha]
else:
boundaries = [(0, 1)]
initial_values = [0.3] # FIXME: should add alpha0 option
initial_values = [0.3] # FIXME: should add alpha0 option
if beta != None:
# allows us to fix beta
......@@ -194,7 +193,7 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
initial_values.append(beta)
else:
boundaries.append((0, 1))
initial_values.append(0.1) # FIXME: should add beta0 option
initial_values.append(0.1) # FIXME: should add beta0 option
if gamma != None:
# allows us to fix gamma
......@@ -202,7 +201,7 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
initial_values.append(gamma)
else:
boundaries.append((0, 1))
initial_values.append(0.1) # FIXME: should add gamma0 option
initial_values.append(0.1) # FIXME: should add gamma0 option
if phi != None:
# allows us to fix phi
......@@ -210,7 +209,7 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
initial_values.append(phi)
else:
boundaries.append((0, 1))
initial_values.append(0.9) # FIXME: should add phi0 option
initial_values.append(0.9) # FIXME: should add phi0 option
initial_values = np.array(initial_values)
method = 'additive'
......@@ -231,14 +230,9 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
# r equal to mean(s)
r = [np.nanmean(s)]
# determine h-step vector of phis for damped trends
# NOTE: Did away with phiVec altogether, and just use phiHminus1 now;
#
#phiVec = np.array([phi**i for i in range(1,hstep)])
# determine sum(c^2) and phi_(j-1) for hstep "prediction interval" outside
# of
# loop; initialize variables for jstep (beyond hstep) prediction intervals
# of loop; initialize variables for jstep (beyond hstep) prediction
# intervals
sumc2_H = 1
phiHminus1 = 0
for h in range(1, hstep):
......@@ -250,7 +244,6 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
jstep = hstep
# convert to, and pre-allocate numpy arrays
# FIXME: this should just be done when checking inputs above
yobs = np.array(yobs)
sigma = np.concatenate((sigma, np.zeros(yobs.size + fc)))
yhat = np.concatenate((yhat, np.zeros(yobs.size + fc)))
......@@ -259,10 +252,9 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
# smooth/simulate/forecast yobs
for i in range(len(yobs) + fc):
# update/append sigma for h steps ahead of i following
# Hyndman-et-al-2005
# NOTE: this will be over-written if valid observations exist at step i
# Update/append sigma for h steps ahead of i following
# Hyndman-et-al-2005. This will be over-written if valid observations
# exist at step i
if jstep == hstep:
sigma2 = sigma[i] * sigma[i]
sigma[i + hstep + 1] = np.sqrt(sigma2 * sumc2)
......@@ -270,28 +262,14 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
# predict h steps ahead
yhat[i + hstep] = l + phiHminus1 * b + s[i + hstep % m]
## NOTE: this was a misguided attempt to smooth s that led to
## oscillatory
## behavior; this makes perfect sense in hindsight, but I'm
## leaving
## comments here as a reminder to NOT try this again. -EJR 6/2015
# yhat[i+hstep] = l + (phiVec*b).sum() + np.nanmean(s[i+ssIdx])
# determine discrepancy between observation and prediction at step i
# FIXME: this if-block becomes unneccessary if we remove the fc option,
# and force the user to place NaNs at the end of yobs if/when
# they want forecasts beyond the last available observation
if i < len(yobs):
et = yobs[i] - yhat[i]
else:
et = np.nan
# this if/else block is not strictly necessary, but it makes the logic
# somewhat easier to follow (for me at least -EJR 5/2015)
if (np.isnan(et) or np.abs(et) > zthresh * sigma[i]):
#
# forecast (i.e., update l, b, and s assuming et==0)
#
# no change in seasonal adjustments
r[i + 1] = 0
......@@ -303,27 +281,23 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
if np.isnan(et):
# when forecasting, grow sigma=sqrt(var) like a prediction
# interval;
# sumc2 and jstep will be reset with the next valid observation
# interval; sumc2 and jstep will be reset with the next
# valid observation
phiJminus1 = phiJminus1 + phi**jstep
sumc2 = sumc2 + (alpha * (1 + phiJminus1 * beta) +
gamma * (1 if (jstep % m == 0) else 0))**2
jstep = jstep + 1
else:
# still update sigma using et when et > zthresh*sigma
# still update sigma using et when et > zthresh * sigma
# (and is not NaN)
# NOTE: Bodenham-et-Adams-2013 may have a more robust method
sigma[i + 1] = alpha * np.abs(et) + (1 - alpha) * sigma[i]
else:
#
# smooth (i.e., update l, b, and s by filtering et)
#
# renormalization could occur inside loop, but we choose to
# integrate
# r, and adjust a and s outside the loop to improve performance.
# integrate r, and adjust a and s outside the loop to improve
# performance.
r[i + 1] = gamma * (1 - alpha) * et / m
# update and append to s using equation-error formulation
......@@ -334,29 +308,14 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
b = phi * b + alpha * beta * et
# update sigma with et, then reset prediction interval variables
# NOTE: Bodenham-et-Adams-2013 may have a more robust method
sigma[i + 1] = alpha * np.abs(et) + (1 - alpha) * sigma[i]
sumc2 = sumc2_H
phiJminus1 = phiHminus1
jstep = hstep
# endif (np.isnan(et) or np.abs(et) > zthresh * sigma[i])
# endfor i in range(len(yobs) + fc - hstep)
"""
NOTE: Seasonal adjustments s[i+1:i+m] should be normalized so their mean
is zero, at least until the next observation, or else the notion of a
"seasonal" adjustment loses all meaning. In order to ensure that the
predictions yhat[:] remain unchanged, the baseline a is shifted too.
Archibald-et-Koehler-2003 recommend doing all this inside the loop,
but this slows the algorithm significantly. A&K-2003 note, however,
that r can be integrated, and used to adjust s[:] *outside* the loop,
and Gardner-2006 recommends this approach. A&K-2003 provide valid
reasons for their recommendation (online optimization of alpha will
be impacted), but since ours is not currently such an estimator, we
choose the more computationally efficient approach.
"""
r = np.cumsum(r)
l = l + r[-1]
s = list(np.array(s) - np.hstack((r, np.tile(r[-1], m - 1))))
......@@ -387,10 +346,4 @@ def additive(yobs, m, alpha=None, beta=None, gamma=None, phi=1,
if __name__ == '__main__':
"""
This might be expanded to call HoltWinters.py as a script. More likely,
HoltWinters.py will be incorporated into another module or class, which
will have it's own command-line functionality.
"""
print 'HELLO'
pass
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment