Fix regularized incomplete beta function (#440)

Author: ishii-norimi

lib/model/tightest_perceptron.js

@@ -30,10 +30,9 @@ const logGamma = z => {
 	return x
 }
 
-const beta = (p, q) => {
+const logBeta = (p, q) => {
 	// https://www2.math.kyushu-u.ac.jp/~snii/AdvancedCalculus/7-1.pdf
-	// return gamma(p) * gamma(q) / gamma(p + q)
-	return Math.exp(logGamma(p) + logGamma(q) - logGamma(p + q))
+	return logGamma(p) + logGamma(q) - logGamma(p + q)
 }
 
 const hypergeometric = (a, b, c, z) => {
@@ -52,34 +51,37 @@ const hypergeometric = (a, b, c, z) => {
 	return f
 }
 
-const incompleteBeta = (z, a, b) => {
+const logIncompleteBeta = (z, a, b) => {
 	// https://ja.wikipedia.org/wiki/%E4%B8%8D%E5%AE%8C%E5%85%A8%E3%83%99%E3%83%BC%E3%82%BF%E9%96%A2%E6%95%B0
 	// https://math-functions-1.watson.jp/sub1_spec_050.html#section030
 	// https://qiita.com/moriokumura/items/e35025d4ade312b0a017
 	if (b === 1) {
-		return z ** a / a
+		return Math.log(z) * a - Math.log(a)
 	} else if (a === 1) {
-		return (1 - (1 - z) ** b) / b
+		return Math.log(1 - (1 - z) ** b) - Math.log(b)
 	} else if (a === 0.5 && b === 0) {
-		return 2 * Math.atanh(Math.sqrt(z))
+		return Math.log(2 * Math.atanh(Math.sqrt(z)))
 	} else if (a === 0.5 && b === 0.5) {
-		return 2 * Math.asin(Math.sqrt(z))
+		return Math.log(2 * Math.asin(Math.sqrt(z)))
 	} else if (Number.isInteger(b)) {
 		const za = z ** a
 		let ib = za / a
 		for (let i = 1; i < b; i++) {
 			ib = (i * ib + za * (1 - z) ** i) / (a + i)
 		}
-		return ib
+		return Math.log(ib)
 	} else if (Number.isInteger(a)) {
 		const zb = (1 - z) ** b
 		let ib = (1 - zb) / b
 		for (let i = 1; i < a; i++) {
 			ib = (i * ib - z ** i * zb) / (i + b)
 		}
-		return ib
+		return Math.log(ib)
 	}
-	return (z ** a / a) * hypergeometric(a, 1 - b, a + 1, z)
+	if (a < b) {
+		return Math.log(Math.exp(logBeta(a, b)) - Math.exp(logIncompleteBeta(1 - z, b, a)))
+	}
+	return Math.log(z) * a - Math.log(a) + b * Math.log(1 - z) + Math.log(hypergeometric(a + b, 1, a + 1, z))
 }
 
 const regularizedIncompleteBeta = (z, a, b) => {
@@ -93,7 +95,10 @@ const regularizedIncompleteBeta = (z, a, b) => {
 	} else if (a === 1) {
 		return 1 - (1 - z) ** b
 	}
-	return incompleteBeta(z, a, b) / beta(a, b)
+	if (a < b) {
+		return 1 - regularizedIncompleteBeta(1 - z, b, a)
+	}
+	return Math.exp(logIncompleteBeta(z, a, b) - logBeta(a, b))
 }
 
 /**

tests/lib/model/tightest_perceptron.test.js

@@ -9,7 +9,7 @@ import { accuracy } from '../../../lib/evaluate/classification.js'
 describe('classification', () => {
 	describe.each([undefined, 'zero_one', 'hinge'])('accuracyLoss %s', accuracyLoss => {
 		test.each([undefined, 'gaussian'])('kernel %s', kernel => {
-			const model = new TightestPerceptron(10, kernel)
+			const model = new TightestPerceptron(10, kernel, accuracyLoss)
 			const x = Matrix.concat(Matrix.randn(50, 2, 0, 0.2), Matrix.randn(50, 2, 5, 0.2)).toArray()
 			x[50] = [0.1, 0.1]
 			const t = []
@@ -25,26 +25,11 @@ describe('classification', () => {
 			expect(acc).toBeGreaterThan(0.9)
 		})
 
-		test.each(['polynomial'])('kernel %s', kernel => {
-			const model = new TightestPerceptron(10, kernel)
-			const x = Matrix.concat(Matrix.randn(50, 2, 0, 0.2), Matrix.randn(50, 2, 5, 0.2)).toArray()
-			x[50] = [0.1, 0.1]
-			const t = []
-			for (let i = 0; i < x.length; i++) {
-				t[i] = Math.floor(i / 50) * 2 - 1
-			}
-			model.init(x, t)
-			for (let i = 0; i < 10; i++) {
-				model.fit()
-			}
-			const y = model.predict(x)
-			const acc = accuracy(y, t)
-			expect(acc).toBeGreaterThan(0.7)
-		})
-
 		test('custom kernel', () => {
-			const model = new TightestPerceptron(10, (a, b) =>
-				Math.exp(-(a.reduce((s, v, i) => s + (v - b[i]) ** 2, 0) ** 2) / 0.01)
+			const model = new TightestPerceptron(
+				10,
+				(a, b) => Math.exp(-(a.reduce((s, v, i) => s + (v - b[i]) ** 2, 0) ** 2) / 0.01),
+				accuracyLoss
 			)
 			const x = Matrix.concat(Matrix.randn(50, 2, 0, 0.2), Matrix.randn(50, 2, 5, 0.2)).toArray()
 			const t = []
@@ -60,4 +45,39 @@ describe('classification', () => {
 			expect(acc).toBeGreaterThan(0.9)
 		})
 	})
+
+	test.each(['polynomial'])('kernel %s', kernel => {
+		const model = new TightestPerceptron(10, kernel)
+		const x = Matrix.concat(Matrix.randn(50, 2, 0, 0.2), Matrix.randn(50, 2, 5, 0.2)).toArray()
+		x[50] = [0.1, 0.1]
+		const t = []
+		for (let i = 0; i < x.length; i++) {
+			t[i] = Math.floor(i / 50) * 2 - 1
+		}
+		model.init(x, t)
+		for (let i = 0; i < 10; i++) {
+			model.fit()
+		}
+		const y = model.predict(x)
+		const acc = accuracy(y, t)
+		expect(acc).toBeGreaterThan(0.7)
+	})
+
+	test('regularizedIncompleteBeta', () => {
+		const model = new TightestPerceptron(10)
+		model._ap = 1.1
+		const x = Matrix.concat(Matrix.randn(50, 2, 0, 0.2), Matrix.randn(50, 2, 5, 0.2)).toArray()
+		x[50] = [0.1, 0.1]
+		const t = []
+		for (let i = 0; i < x.length; i++) {
+			t[i] = Math.floor(i / 50) * 2 - 1
+		}
+		model.init(x, t)
+		for (let i = 0; i < 10; i++) {
+			model.fit()
+		}
+		const y = model.predict(x)
+		const acc = accuracy(y, t)
+		expect(acc).toBeGreaterThan(0.7)
+	})
 })