Github user myui commented on a diff in the pull request:
https://github.com/apache/incubator-hivemall/pull/92#discussion_r125170752
--- Diff: core/src/main/java/hivemall/anomaly/SDAR2D.java ---
@@ -33,129 +33,123 @@
public final class SDAR2D {
- // --------------------------------
- // hyperparameters
-
- /**
- * Forgetfulness parameter
- */
- private final double _r;
+ // --------------------------------
+ // hyperparameters
+
+ /**
+ * Forgetfulness parameter
+ */
+ private final double _r;
+
+ // --------------------------------
+ // parameters
+
+ private final RealMatrix[] _C;
+ private RealVector _mu;
+ private RealMatrix _sigma;
+ private RealVector _muOld;
+ private RealMatrix _sigmaOld;
+
+ private boolean _initialized;
+
+ public SDAR2D(final double r, final int k) {
+ Preconditions.checkArgument(0.d < r && r < 1.d, "Invalid
forgetfullness parameter r: " + r);
+ Preconditions.checkArgument(k >= 1, "Invalid smoothing parameter k: "
+ k);
+ this._r = r;
+ this._C = new RealMatrix[k + 1];
+ this._initialized = false;
+ }
+
+ /**
+ * @param x series of input in LIFO order
+ * @param k AR window size
+ * @return x_hat predicted x
+ * @link
https://en.wikipedia.org/wiki/Matrix_multiplication#Outer_product
+ */
+ @Nonnull
+ public RealVector update(@Nonnull final ArrayRealVector[] x, final int
k) {
+ Preconditions.checkArgument(x.length >= 1, "x.length MUST be greather
than 1: " + x.length);
+ Preconditions.checkArgument(k >= 0, "k MUST be greather than or equals
to 0: ", k);
+ Preconditions.checkArgument(k < _C.length, "k MUST be less than |C|
but " + "k=" + k + ", |C|="
+ + _C.length);
+
+ final ArrayRealVector x_t = x[0];
+ final int dims = x_t.getDimension();
+
+ if (_initialized == false) {
+ this._mu = x_t.copy();
+ this._sigma = new BlockRealMatrix(dims, dims);
+ assert (_sigma.isSquare());
+ this._initialized = true;
+ return new ArrayRealVector(dims);
+ }
+ Preconditions.checkArgument(k >= 1, "k MUST be greater than 0: ", k);
- // --------------------------------
- // parameters
+ // old parameters are accessible to compute the Hellinger distance
+ this._muOld = _mu.copy();
+ this._sigmaOld = _sigma.copy();
- private final RealMatrix[] _C;
- private RealVector _mu;
- private RealMatrix _sigma;
- private RealVector _muOld;
- private RealMatrix _sigmaOld;
+ // update mean vector
+ // \hat{mu} := (1-r) \hat{µ} + r x_t
+ this._mu = _mu.mapMultiply(1.d - _r).add(x_t.mapMultiply(_r));
- private boolean _initialized;
+ // compute residuals (x - \hat{µ})
+ final RealVector[] xResidual = new RealVector[k + 1];
+ for (int j = 0; j <= k; j++) {
+ xResidual[j] = x[j].subtract(_mu);
+ }
- public SDAR2D(final double r, final int k) {
- Preconditions.checkArgument(0.d < r && r < 1.d, "Invalid
forgetfullness parameter r: " + r);
- Preconditions.checkArgument(k >= 1, "Invalid smoothing parameter
k: " + k);
- this._r = r;
- this._C = new RealMatrix[k + 1];
- this._initialized = false;
+ // update covariance matrices
+ // C_j := (1-r) C_j + r (x_t - \hat{µ}) (x_{t-j} - \hat{µ})'
+ final RealMatrix[] C = this._C;
+ final RealVector rxResidual0 = xResidual[0].mapMultiply(_r); // r (x_t
- \hat{µ})
+ for (int j = 0; j <= k; j++) {
+ RealMatrix Cj = C[j];
+ if (Cj == null) {
+ C[j] = rxResidual0.outerProduct(x[j].subtract(_mu));
+ } else {
+ C[j] = Cj.scalarMultiply(1.d -
_r).add(rxResidual0.outerProduct(x[j].subtract(_mu)));
+ }
}
- /**
- * @param x series of input in LIFO order
- * @param k AR window size
- * @return x_hat predicted x
- * @link
https://en.wikipedia.org/wiki/Matrix_multiplication#Outer_product
+ // solve A in the following Yule-Walker equation
+ // C_j = â_{i=1}^{k} A_i C_{j-i} where j = 1..k, C_{-i} = C_i'
+ /*
+ * /C_1\ /A_1\ /C_0 |C_1' |C_2' | . . . |C_{k-1}' \ |---| |---|
|--------+--------+--------+
--- End diff --
This comment should not be formatted.
---
If your project is set up for it, you can reply to this email and have your
reply appear on GitHub as well. If your project does not have this feature
enabled and wishes so, or if the feature is enabled but not working, please
contact infrastructure at [email protected] or file a JIRA ticket
with INFRA.
---
