Fix Markov switching

js/view/markov_switching.js

@@ -6,8 +6,8 @@ var dispMSM = function (elm, platform) {
 		platform.fit((tx, ty, pred_cb, thup) => {
 			const regime = +elm.select('[name=regime]').property('value')
 			const trial = +elm.select('[name=trial]').property('value')
-			const model = new MarkovSwitching(regime, tx[0].length)
-			model.fit(tx, 0.1, trial)
+			const model = new MarkovSwitching(regime)
+			model.fit(tx, 1, trial)
 			const threshold = +elm.select('[name=threshold]').property('value')
 			const pred = model.predict(tx)
 			thupdater = thup

lib/model/markov_switching.js

@@ -7,18 +7,11 @@ export default class MarkovSwitching {
 	// https://qiita.com/9uant/items/0abf942fac26aee1fc3f
 	/**
 	 * @param {number} regime
-	 * @param {number} d
 	 */
-	constructor(regime, d) {
+	constructor(regime) {
 		this._regime = regime
-		this._d = d
 		this._mu = []
 		this._sigma = []
-		for (let i = 0; i < this._regime; i++) {
-			this._mu[i] = Matrix.randn(1, d)
-			this._sigma[i] = Matrix.random(d, d, 0, 0.1)
-			this._sigma[i] = this._sigma[i].tDot(this._sigma[i])
-		}
 	}
 
 	_stationary_prob(p) {
@@ -33,13 +26,14 @@ export default class MarkovSwitching {
 		const n = x.rows
 		const lh = new Matrix(n, this._regime)
 		for (let i = 0; i < this._regime; i++) {
-			const dv = Math.pow(2 * Math.PI, this._d / 2) * Math.sqrt(sigma[i].det())
+			const dv = Math.pow(2 * Math.PI, mu[i].cols / 2) * Math.sqrt(sigma[i].det())
 			const sinv = sigma[i].inv()
-			const sx = Matrix.sub(x, mu[i])
-			for (let k = 0; k < n; k++) {
-				const s = sx.row(k)
-				lh.set(k, i, Math.exp(-s.dot(sinv).dot(s.t).toScaler() / 2) / dv)
-			}
+			const cx = Matrix.sub(x, mu[i])
+			const xsx = cx.dot(sinv)
+			xsx.mult(cx)
+			const s = xsx.sum(1)
+			s.map(v => Math.exp(-v / 2) / dv)
+			lh.set(0, i, s)
 		}
 		return lh
 	}
@@ -86,33 +80,30 @@ export default class MarkovSwitching {
 		const ngp = genProb.copy()
 
 		for (let i = 0; i < this._regime; i++) {
-			nmu[i].add(Math.random(1, this._d, -eps, eps))
-			for (let j = 0; j < this._d; j++) {
-				for (let k = 0; k <= j; k++) {
-					let s = nsi[i].at(j, k)
-					s = Math.exp(Math.log(s) + (2 * Math.random() - 1) * eps)
-					nsi[i].set(j, k, s)
-					nsi[i].set(k, j, s)
-				}
+			nmu[i].add(Matrix.random(1, nmu[i].cols, -eps, eps))
+			for (let j = 0; j < nmu[i].cols; j++) {
+				let s = nsi[i].at(j, j)
+				s = Math.exp(Math.log(s) + (2 * Math.random() - 1) * eps)
+				nsi[i].set(j, j, s)
 			}
 		}
 
 		ngp.add(Matrix.map(Matrix.random(this._regime, this._regime), v => (2 * v - 1) * eps * 0.1))
 		for (let i = 0; i < this._regime; i++) {
 			ngp.set(i, i, 0)
 		}
-		const np = Matrix.map(genProb, v => Math.exp(v) / (1 + v))
+		const np = Matrix.map(genProb, v => Math.exp(v) / (1 + Math.exp(v)))
 		np.add(Matrix.diag(Matrix.sub(1, np.sum(0).t).value))
 
 		return [nmu, nsi, ngp, np]
 	}
 
 	_mcmc(x, eps, trial) {
-		let genProb = new Matrix(this._regime, this._regime, -3)
+		let genProb = new Matrix(this._regime, this._regime, -this._regime)
 		for (let i = 0; i < this._regime; i++) {
 			genProb.set(i, i, 0)
 		}
-		let prob = Matrix.map(genProb, v => Math.exp(v) / (1 + v))
+		let prob = Matrix.map(genProb, v => Math.exp(v) / (1 + Math.exp(v)))
 		prob.add(Matrix.diag(Matrix.sub(1, prob.sum(0).t).value))
 
 		const mus = [this._mu.concat()]
@@ -153,9 +144,30 @@ export default class MarkovSwitching {
 	 */
 	fit(datas, eps, trial) {
 		const x = Matrix.fromArray(datas)
+		if (this._mu.length === 0) {
+			for (let i = 0; i < this._regime; i++) {
+				this._mu[i] = Matrix.randn(1, x.cols, 0, 0.1)
+				this._sigma[i] = Matrix.zeros(x.cols, x.cols)
+				for (let j = 0; j < x.cols; j++) {
+					this._sigma[i].set(j, j, Math.random())
+				}
+			}
+		}
 		const [ms, ss, ps, ls] = this._mcmc(x, eps, trial)
 	}
 
+	/**
+	 * Returns probabilities.
+	 *
+	 * @param {Array<Array<number>>} datas
+	 * @returns {Array<Array<number>>}
+	 */
+	probability(datas) {
+		const x = Matrix.fromArray(datas)
+		const probs = this._prob(x)
+		return probs.map(p => p.value)
+	}
+
 	/**
 	 * Returns anomaly degrees.
 	 *
@@ -165,10 +177,11 @@ export default class MarkovSwitching {
 	predict(datas) {
 		const x = Matrix.fromArray(datas)
 		const probs = this._prob(x)
+		const norms = probs.map(p => p.norm())
 
 		const pred = []
 		for (let i = 0; i < probs.length - 1; i++) {
-			pred.push(1 - Matrix.mult(probs[i], probs[i + 1]).sum() / (probs[i].norm() * probs[i + 1].norm()))
+			pred.push(1 - Matrix.mult(probs[i], probs[i + 1]).sum() / (norms[i] * norms[i + 1]))
 		}
 		return pred
 	}

tests/lib/model/markov_switching.test.js

@@ -0,0 +1,28 @@
+import { jest } from '@jest/globals'
+jest.retryTimes(3)
+
+import Matrix from '../../../lib/util/matrix.js'
+import MarkovSwitching from '../../../lib/model/markov_switching.js'
+
+test('anomaly detection', () => {
+	const model = new MarkovSwitching(2)
+	const n = 50
+	const x = Matrix.concat(Matrix.random(n, 2, 0, 1), Matrix.random(n, 2, 2, 3)).toArray()
+
+	model.fit(x, 1, 2000)
+	const p = model.predict(x)
+
+	const threshold = 0.1
+	const range = 5
+	let c = 0
+	for (let i = 0; i < p.length; i++) {
+		if (i < n - range || n + range < i) {
+			expect(p[i]).toBeLessThan(threshold)
+		} else {
+			if (p[i] >= threshold) {
+				c++
+			}
+		}
+	}
+	expect(c).toBeGreaterThan(0)
+})