|
| 1 | +const logGamma = z => { |
| 2 | + // https://en.wikipedia.org/wiki/Lanczos_approximation |
| 3 | + // https://slpr.sakura.ne.jp/qp/gamma-function/ |
| 4 | + let x = 0 |
| 5 | + if (Number.isInteger(z)) { |
| 6 | + for (let i = 2; i < z; i++) { |
| 7 | + x += Math.log(i) |
| 8 | + } |
| 9 | + } else if (Number.isInteger(z - 0.5)) { |
| 10 | + const n = z - 0.5 |
| 11 | + x = Math.log(Math.sqrt(Math.PI)) - Math.log(2) * n |
| 12 | + for (let i = 2 * n - 1; i > 0; i -= 2) { |
| 13 | + x += Math.log(i) |
| 14 | + } |
| 15 | + } else if (z < 0.5) { |
| 16 | + x = Math.log(Math.PI) - Math.log(Math.sin(Math.PI * z)) - logGamma(1 - z) |
| 17 | + } else { |
| 18 | + const p = [ |
| 19 | + 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, |
| 20 | + -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7, |
| 21 | + ] |
| 22 | + z -= 1 |
| 23 | + x = 0.99999999999980993 |
| 24 | + for (let i = 0; i < p.length; i++) { |
| 25 | + x += p[i] / (z + i + 1) |
| 26 | + } |
| 27 | + const t = z + p.length - 0.5 |
| 28 | + x = Math.log(Math.sqrt(2 * Math.PI)) + Math.log(t) * (z + 0.5) - t + Math.log(x) |
| 29 | + } |
| 30 | + return x |
| 31 | +} |
| 32 | + |
| 33 | +const beta = (p, q) => { |
| 34 | + // https://www2.math.kyushu-u.ac.jp/~snii/AdvancedCalculus/7-1.pdf |
| 35 | + // return gamma(p) * gamma(q) / gamma(p + q) |
| 36 | + return Math.exp(logGamma(p) + logGamma(q) - logGamma(p + q)) |
| 37 | +} |
| 38 | + |
| 39 | +const hypergeometric = (a, b, c, z) => { |
| 40 | + // https://qiita.com/moriokumura/items/e35025d4ade312b0a017 |
| 41 | + let f = 1 |
| 42 | + let p = 0 |
| 43 | + const lnz = Math.log(z) |
| 44 | + for (let n = 0; n < 1000; n++) { |
| 45 | + p = p + lnz + Math.log(a + n) + Math.log(b + n) - Math.log(c + n) - Math.log(1 + n) |
| 46 | + const ep = Math.exp(p) |
| 47 | + f += ep |
| 48 | + if (Math.abs(ep / f) < 1.0e-14) { |
| 49 | + break |
| 50 | + } |
| 51 | + } |
| 52 | + return f |
| 53 | +} |
| 54 | + |
| 55 | +const incompleteBeta = (z, a, b) => { |
| 56 | + // 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 |
| 57 | + // https://math-functions-1.watson.jp/sub1_spec_050.html#section030 |
| 58 | + // https://qiita.com/moriokumura/items/e35025d4ade312b0a017 |
| 59 | + if (b === 1) { |
| 60 | + return z ** a / a |
| 61 | + } else if (a === 1) { |
| 62 | + return (1 - (1 - z) ** b) / b |
| 63 | + } else if (a === 0.5 && b === 0) { |
| 64 | + return 2 * Math.atanh(Math.sqrt(z)) |
| 65 | + } else if (a === 0.5 && b === 0.5) { |
| 66 | + return 2 * Math.asin(Math.sqrt(z)) |
| 67 | + } else if (Number.isInteger(b)) { |
| 68 | + const za = z ** a |
| 69 | + let ib = za / a |
| 70 | + for (let i = 1; i < b; i++) { |
| 71 | + ib = (i * ib + za * (1 - z) ** i) / (a + i) |
| 72 | + } |
| 73 | + return ib |
| 74 | + } else if (Number.isInteger(a)) { |
| 75 | + const zb = (1 - z) ** b |
| 76 | + let ib = (1 - zb) / b |
| 77 | + for (let i = 1; i < a; i++) { |
| 78 | + ib = (i * ib - z ** i * zb) / (i + b) |
| 79 | + } |
| 80 | + return ib |
| 81 | + } |
| 82 | + return (z ** a / a) * hypergeometric(a, 1 - b, a + 1, z) |
| 83 | +} |
| 84 | + |
| 85 | +const regularizedIncompleteBeta = (z, a, b) => { |
| 86 | + // beta distribution of the first kind |
| 87 | + if (z === 0) { |
| 88 | + return 0 |
| 89 | + } else if (z === 1) { |
| 90 | + return 1 |
| 91 | + } else if (b === 1) { |
| 92 | + return z ** a |
| 93 | + } else if (a === 1) { |
| 94 | + return 1 - (1 - z) ** b |
| 95 | + } |
| 96 | + return incompleteBeta(z, a, b) / beta(a, b) |
| 97 | +} |
| 98 | + |
| 99 | +/** |
| 100 | + * Tightest Perceptron |
| 101 | + */ |
| 102 | +export default class TightestPerceptron { |
| 103 | + // Online Learning: A Comprehensive Survey |
| 104 | + // https://arxiv.org/abs/1802.02871 |
| 105 | + // Tighter perceptron with improved dual use of cached data for model representation and validation. |
| 106 | + // https://www.dabi.temple.edu/external/vucetic/documents/wang09ijcnn.pdf |
| 107 | + /** |
| 108 | + * @param {number} [b=10] Budget size |
| 109 | + * @param {'gaussian' | 'polynomial' | function (number[], number[]): number} [kernel=gaussian] Kernel name |
| 110 | + * @param {'zero_one' | 'hinge'} [accuracyLoss=hinge] Accuracy loss type name |
| 111 | + */ |
| 112 | + constructor(b = 10, kernel = 'gaussian', accuracyLoss = 'hinge') { |
| 113 | + this._b = b |
| 114 | + |
| 115 | + if (typeof kernel === 'function') { |
| 116 | + this._kernel = kernel |
| 117 | + } else { |
| 118 | + switch (kernel) { |
| 119 | + case 'gaussian': |
| 120 | + this._s = 1 |
| 121 | + this._kernel = (a, b) => |
| 122 | + Math.exp(-(a.reduce((s, v, i) => s + (v - b[i]) ** 2, 0) ** 2) / this._s ** 2) |
| 123 | + break |
| 124 | + case 'polynomial': |
| 125 | + this._d = 2 |
| 126 | + this._kernel = (a, b) => (1 + a.reduce((s, v, i) => s + v * b[i])) ** this._d |
| 127 | + break |
| 128 | + } |
| 129 | + } |
| 130 | + |
| 131 | + if (accuracyLoss === 'hinge') { |
| 132 | + this._accuracyLossP = y => { |
| 133 | + return Math.max(0, 1 - y) |
| 134 | + } |
| 135 | + this._accuracyLossN = y => { |
| 136 | + return Math.max(0, 1 + y) |
| 137 | + } |
| 138 | + } else { |
| 139 | + this._accuracyLossP = y => { |
| 140 | + return y < 0 ? 1 : 0 |
| 141 | + } |
| 142 | + this._accuracyLossN = y => { |
| 143 | + return y < 0 ? 0 : 1 |
| 144 | + } |
| 145 | + } |
| 146 | + this._ap = 1 |
| 147 | + this._an = 1 |
| 148 | + } |
| 149 | + |
| 150 | + /** |
| 151 | + * Initialize this model. |
| 152 | + * |
| 153 | + * @param {Array<Array<number>>} train_x Training data |
| 154 | + * @param {Array<1 | -1>} train_y Target values |
| 155 | + */ |
| 156 | + init(train_x, train_y) { |
| 157 | + this._x = train_x |
| 158 | + this._y = train_y |
| 159 | + |
| 160 | + this._sv = [] |
| 161 | + } |
| 162 | + |
| 163 | + /** |
| 164 | + * Fit model parameters. |
| 165 | + */ |
| 166 | + fit() { |
| 167 | + for (let i = 0; i < this._x.length; i++) { |
| 168 | + let s = 0 |
| 169 | + for (let k = 0; k < this._sv.length; k++) { |
| 170 | + const sk = this._sv[k] |
| 171 | + s += sk.y * this._kernel(this._x[i], sk.x) |
| 172 | + } |
| 173 | + if (s * this._y[i] <= 0) { |
| 174 | + if (this._y[i] === 1) { |
| 175 | + this._sv.push({ x: this._x[i], y: this._y[i], cp: 1, cn: 0 }) |
| 176 | + } else { |
| 177 | + this._sv.push({ x: this._x[i], y: this._y[i], cp: 0, cn: 1 }) |
| 178 | + } |
| 179 | + if (this._sv.length > this._b) { |
| 180 | + let min_l = Infinity |
| 181 | + let min_r = -1 |
| 182 | + for (let k = 0; k < this._sv.length; k++) { |
| 183 | + let loss = 0 |
| 184 | + for (let j = 0; j < this._sv.length; j++) { |
| 185 | + let f = 0 |
| 186 | + for (let m = 0; m < this._sv.length; m++) { |
| 187 | + if (m === k) { |
| 188 | + continue |
| 189 | + } |
| 190 | + const sk = this._sv[m] |
| 191 | + f += sk.y * this._kernel(this._sv[j].x, sk.x) |
| 192 | + } |
| 193 | + const lp = this._accuracyLossP(f) |
| 194 | + const ln = this._accuracyLossN(f) |
| 195 | + const wp = |
| 196 | + 1 - regularizedIncompleteBeta(0.5, this._sv[j].cp + this._ap, this._sv[j].cn + this._an) |
| 197 | + loss += wp * lp + (1 - wp) * ln |
| 198 | + } |
| 199 | + loss /= this._sv.length |
| 200 | + if (loss < min_l) { |
| 201 | + min_l = loss |
| 202 | + min_r = k |
| 203 | + } |
| 204 | + } |
| 205 | + const sv = this._sv.splice(min_r, 1)[0] |
| 206 | + this._updateSummary(sv.x, sv.cp, sv.cn) |
| 207 | + } |
| 208 | + } else { |
| 209 | + if (this._y[i] === 1) { |
| 210 | + this._updateSummary(this._x[i], 1, 0) |
| 211 | + } else { |
| 212 | + this._updateSummary(this._x[i], 0, 1) |
| 213 | + } |
| 214 | + } |
| 215 | + } |
| 216 | + } |
| 217 | + |
| 218 | + _updateSummary(x, cp, cn) { |
| 219 | + if (this._sv.length === 0) { |
| 220 | + return |
| 221 | + } |
| 222 | + let min_d = Infinity |
| 223 | + let min_k = -1 |
| 224 | + for (let i = 0; i < this._sv.length; i++) { |
| 225 | + const d = this._sv[i].x.reduce((s, v, d) => s + (v - x[d]) ** 2, 0) |
| 226 | + if (d < min_d) { |
| 227 | + min_d = d |
| 228 | + min_k = i |
| 229 | + } |
| 230 | + } |
| 231 | + const kn = this._kernel(x, this._sv[min_k].x) |
| 232 | + this._sv[min_k].cp += cp * kn |
| 233 | + this._sv[min_k].cn += cn * kn |
| 234 | + } |
| 235 | + |
| 236 | + /** |
| 237 | + * Returns predicted values. |
| 238 | + * |
| 239 | + * @param {Array<Array<number>>} data Sample data |
| 240 | + * @returns {(1 | -1)[]} Predicted values |
| 241 | + */ |
| 242 | + predict(data) { |
| 243 | + const p = [] |
| 244 | + for (let i = 0; i < data.length; i++) { |
| 245 | + let s = 0 |
| 246 | + for (let k = 0; k < this._sv.length; k++) { |
| 247 | + const sk = this._sv[k] |
| 248 | + s += sk.y * this._kernel(data[i], sk.x) |
| 249 | + } |
| 250 | + p[i] = s < 0 ? -1 : 1 |
| 251 | + } |
| 252 | + return p |
| 253 | + } |
| 254 | +} |
0 commit comments