Estimate bias parameter in LTS (#766)

Author: ishii-norimi

lib/model/lts.js

@@ -14,8 +14,13 @@ export default class LeastTrimmedSquaresRegression {
 	}
 
 	_ls(x, y) {
+		const m = x.mean(0)
+		x = Matrix.sub(x, m)
 		const xtx = x.tDot(x)
-		return xtx.solve(x.tDot(y))
+		const w = xtx.solve(x.tDot(y))
+		y = Matrix.sub(y, x.dot(w))
+		const b = y.mean(0)
+		return [w, m, b]
 	}
 
 	/**
@@ -27,10 +32,10 @@ export default class LeastTrimmedSquaresRegression {
 	fit(x, y) {
 		x = Matrix.fromArray(x)
 		y = Matrix.fromArray(y)
-		const xh = Matrix.resize(x, x.rows, x.cols + 1, 1)
 
-		const wls = this._ls(xh, y)
-		const yls = xh.dot(wls)
+		const [wls, mls, bls] = this._ls(x, y)
+		const yls = Matrix.sub(x, mls).dot(wls)
+		yls.add(bls)
 		yls.sub(y)
 		yls.mult(yls)
 
@@ -39,10 +44,10 @@ export default class LeastTrimmedSquaresRegression {
 
 		const h = Math.max(1, Math.floor(r.length * this._h))
 
-		const xlts = xh.row(r.slice(0, h).map(v => v[1]))
+		const xlts = x.row(r.slice(0, h).map(v => v[1]))
 		const ylts = y.row(r.slice(0, h).map(v => v[1]))
 
-		this._w = this._ls(xlts, ylts)
+		;[this._w, this._m, this._b] = this._ls(xlts, ylts)
 	}
 
 	/**
@@ -53,7 +58,9 @@ export default class LeastTrimmedSquaresRegression {
 	 */
 	predict(x) {
 		x = Matrix.fromArray(x)
-		const xh = Matrix.resize(x, x.rows, x.cols + 1, 1)
-		return xh.dot(this._w).toArray()
+		x.sub(this._m)
+		const y = x.dot(this._w)
+		y.add(this._b)
+		return y.toArray()
 	}
 }

tests/lib/model/least_square.test.js

@@ -15,3 +15,28 @@ test('fit', () => {
 	const err = rmse(y, t)[0]
 	expect(err).toBeLessThan(0.5)
 })
+
+test('same bias column', () => {
+	const model = new LeastSquare()
+	const x = Matrix.randn(50, 2, 0, 5).toArray()
+	const t = []
+	for (let i = 0; i < x.length; i++) {
+		t[i] = [x[i][0] + x[i][1] + (Math.random() - 0.5) / 10 + 5]
+	}
+	model.fit(x, t)
+	const y = model.predict(x)
+
+	const x2 = Matrix.resize(Matrix.fromArray(x), x.length, x[0].length + 1, 1)
+	const t2 = Matrix.fromArray(t)
+	const w = x2.tDot(x2).solve(x2.tDot(t2))
+	const y2 = x2.dot(w)
+
+	for (let i = 0; i < y.length; i++) {
+		for (let j = 0; j < y[i].length; j++) {
+			expect(y[i][j]).toBeCloseTo(y2.at(i, j))
+		}
+	}
+	for (let i = 0; i < 2; i++) {
+		expect(model._w.at(i, 0)).toBeCloseTo(w.at(i, 0))
+	}
+})

tests/lib/model/lts.test.js

@@ -15,3 +15,36 @@ test('fit', () => {
 	const err = rmse(y, t)[0]
 	expect(err).toBeLessThan(0.5)
 })
+
+test('same bias column', () => {
+	const model = new LeastTrimmedSquaresRegression()
+	const x = Matrix.randn(50, 2, 0, 5).toArray()
+	const t = []
+	for (let i = 0; i < x.length; i++) {
+		t[i] = [x[i][0] + x[i][1] + (Math.random() - 0.5) / 10]
+	}
+	model.fit(x, t)
+	const y = model.predict(x)
+
+	const x2 = Matrix.resize(Matrix.fromArray(x), x.length, x[0].length + 1, 1)
+	const t2 = Matrix.fromArray(t)
+	const d = Matrix.sub(t2, x2.dot(x2.tDot(x2).solve(x2.tDot(t2))))
+	d.map(v => v ** 2)
+	const r = d.sum(1).value.map((v, i) => [v, i])
+	r.sort((a, b) => a[0] - b[0])
+
+	const h = 45
+	const xlts = x2.row(r.slice(0, h).map(v => v[1]))
+	const ylts = t2.row(r.slice(0, h).map(v => v[1]))
+	const w = xlts.tDot(xlts).solve(xlts.tDot(ylts))
+	const y2 = x2.dot(w)
+
+	for (let i = 0; i < y.length; i++) {
+		for (let j = 0; j < y[i].length; j++) {
+			expect(y[i][j]).toBeCloseTo(y2.at(i, j))
+		}
+	}
+	for (let i = 0; i < 2; i++) {
+		expect(model._w.at(i, 0)).toBeCloseTo(w.at(i, 0))
+	}
+})