Update: 2026-07-06 17:51:16
This commit is contained in:
@@ -1,13 +1,5 @@
|
||||
/**
|
||||
* Matrix and statistical utilities for pricing analysis.
|
||||
* Pure math — no external dependencies except simple-statistics.
|
||||
*/
|
||||
import { median, mean } from 'simple-statistics';
|
||||
|
||||
import { median, mean, standardDeviation } from 'simple-statistics';
|
||||
|
||||
/**
|
||||
* Compute Pearson correlation between two arrays.
|
||||
*/
|
||||
function pearsonCorr(x: number[], y: number[]): number {
|
||||
const n = Math.min(x.length, y.length);
|
||||
if (n < 3) return 0;
|
||||
@@ -25,64 +17,53 @@ function pearsonCorr(x: number[], y: number[]): number {
|
||||
return denom === 0 ? 0 : num / denom;
|
||||
}
|
||||
|
||||
/**
|
||||
* Multiple linear regression via Gaussian elimination with ridge regularization.
|
||||
* Solves: price = baseFare + kmRate*distance + minRate*duration
|
||||
*
|
||||
* Uses L2 ridge (lambda=0.1) when distance≈duration are collinear.
|
||||
* Falls back to distance-only model if necessary.
|
||||
*/
|
||||
export function multipleLinearRegression(
|
||||
samples: Array<{ distance_km: number; duration_min: number; price: number }>
|
||||
): { baseFare: number; kmRate: number; minRate: number } | null {
|
||||
const n = samples.length;
|
||||
if (n < 3) return null;
|
||||
|
||||
// Check collinearity: if distance and duration are highly correlated
|
||||
const dists = samples.map(s => s.distance_km);
|
||||
const durs = samples.map(s => s.duration_min);
|
||||
const prices = samples.map(s => s.price);
|
||||
const corr = pearsonCorr(dists, durs);
|
||||
|
||||
const lambda = Math.abs(corr) > 0.85 ? 0.5 : 0.01; // ridge penalty
|
||||
const lambda = Math.abs(corr) > 0.85 ? 0.5 : 0.01;
|
||||
|
||||
let sumX1 = 0, sumX2 = 0, sumY = 0;
|
||||
let sumX1Sq = 0, sumX2Sq = 0, sumX1X2 = 0;
|
||||
let sumX1Y = 0, sumX2Y = 0;
|
||||
|
||||
for (const s of samples) {
|
||||
const x1 = s.distance_km;
|
||||
const x2 = s.duration_min;
|
||||
const y = s.price;
|
||||
|
||||
const x1 = s.distance_km, x2 = s.duration_min, y = s.price;
|
||||
sumX1 += x1; sumX2 += x2; sumY += y;
|
||||
sumX1Sq += x1 * x1; sumX2Sq += x2 * x2; sumX1X2 += x1 * x2;
|
||||
sumX1Y += x1 * y; sumX2Y += x2 * y;
|
||||
}
|
||||
|
||||
// Ridge: add lambda to diagonal of X^T X (except intercept)
|
||||
const A = [
|
||||
[n, sumX1, sumX2],
|
||||
[sumX1, sumX1Sq + lambda, sumX1X2],
|
||||
[sumX2, sumX1X2, sumX2Sq + lambda],
|
||||
[n, sumX1, sumX2],
|
||||
[sumX1, sumX1Sq + lambda, sumX1X2],
|
||||
[sumX2, sumX1X2, sumX2Sq + lambda],
|
||||
];
|
||||
|
||||
const B = [sumY, sumX1Y, sumX2Y];
|
||||
|
||||
try {
|
||||
const beta = gaussianElimination(A, B);
|
||||
const baseFare = Math.max(0, beta[0]);
|
||||
let baseFare = Math.max(0, beta[0]);
|
||||
let kmRate = Math.max(0, beta[1]);
|
||||
let minRate = Math.max(0, beta[2]);
|
||||
|
||||
// If minRate is essentially zero after ridge, keep it minimal
|
||||
if (minRate < 0.001) minRate = 0;
|
||||
const predFull = samples.map(s => baseFare + kmRate * s.distance_km + minRate * s.duration_min);
|
||||
const rmseFull = calcRMSE(prices, predFull);
|
||||
|
||||
// If both non-intercept terms are zero, try distance-only model
|
||||
if (kmRate === 0 && minRate === 0) {
|
||||
const k = sumX1Y / (sumX1Sq + lambda);
|
||||
if (k > 0) {
|
||||
kmRate = k;
|
||||
}
|
||||
const ratioK = dists.reduce((a, d, i) => d > 0 ? a + prices[i] / d : a, 0) / dists.filter(d => d > 0).length;
|
||||
if (!isFinite(ratioK)) return { baseFare, kmRate, minRate };
|
||||
|
||||
const predDistOnly = dists.map(d => ratioK * d);
|
||||
const rmseDist = calcRMSE(prices, predDistOnly);
|
||||
|
||||
if (rmseDist <= rmseFull * 1.10) {
|
||||
return { baseFare: 0, kmRate: Math.round(ratioK * 1000) / 1000, minRate: 0 };
|
||||
}
|
||||
|
||||
return { baseFare, kmRate, minRate };
|
||||
@@ -91,10 +72,6 @@ export function multipleLinearRegression(
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Robust iterative regression to find floor pricing (exclude surge outliers).
|
||||
* Uses a fixed JOD/SYP threshold per iteration instead of tightening RMSE.
|
||||
*/
|
||||
export function robustMultipleLinearRegression(
|
||||
samples: Array<{ distance_km: number; duration_min: number; price: number }>,
|
||||
maxIterations: number = 4
|
||||
@@ -103,7 +80,6 @@ export function robustMultipleLinearRegression(
|
||||
let bestModel = multipleLinearRegression(currentSamples);
|
||||
if (!bestModel) return null;
|
||||
|
||||
// Determine threshold from data scale (median price × 0.3)
|
||||
const prices = samples.map(s => s.price).sort((a, b) => a - b);
|
||||
const medianPrice = prices[Math.floor(prices.length / 2)];
|
||||
const fixedThreshold = Math.max(medianPrice * 0.3, 0.1);
|
||||
@@ -113,214 +89,160 @@ export function robustMultipleLinearRegression(
|
||||
s => bestModel!.baseFare + bestModel!.kmRate * s.distance_km + bestModel!.minRate * s.duration_min
|
||||
);
|
||||
const actual = currentSamples.map(s => s.price);
|
||||
|
||||
const inliers = currentSamples.filter((s, idx) => {
|
||||
const residual = actual[idx] - predicted[idx];
|
||||
return residual < fixedThreshold;
|
||||
});
|
||||
const inliers = currentSamples.filter((_, idx) => actual[idx] - predicted[idx] < fixedThreshold);
|
||||
|
||||
if (inliers.length < Math.max(5, samples.length * 0.3)) break;
|
||||
if (inliers.length === currentSamples.length) break;
|
||||
|
||||
currentSamples = inliers;
|
||||
const newModel = multipleLinearRegression(currentSamples);
|
||||
if (!newModel) break;
|
||||
bestModel = newModel;
|
||||
}
|
||||
|
||||
return bestModel;
|
||||
}
|
||||
|
||||
/**
|
||||
* Gaussian elimination for solving Ax = B (3x3 system).
|
||||
*/
|
||||
function gaussianElimination(A: number[][], B: number[]): number[] {
|
||||
const n = A.length;
|
||||
const a = A.map(row => [...row]);
|
||||
const b = [...B];
|
||||
|
||||
for (let i = 0; i < n; i++) {
|
||||
let maxEl = Math.abs(a[i][i]);
|
||||
let maxRow = i;
|
||||
let maxEl = Math.abs(a[i][i]), maxRow = i;
|
||||
for (let k = i + 1; k < n; k++) {
|
||||
if (Math.abs(a[k][i]) > maxEl) {
|
||||
maxEl = Math.abs(a[k][i]);
|
||||
maxRow = k;
|
||||
}
|
||||
if (Math.abs(a[k][i]) > maxEl) { maxEl = Math.abs(a[k][i]); maxRow = k; }
|
||||
}
|
||||
|
||||
[a[maxRow], a[i]] = [a[i], a[maxRow]];
|
||||
[b[maxRow], b[i]] = [b[i], b[maxRow]];
|
||||
|
||||
if (Math.abs(a[i][i]) < 1e-12) continue;
|
||||
|
||||
for (let k = i + 1; k < n; k++) {
|
||||
const c = -a[k][i] / a[i][i];
|
||||
for (let j = i; j < n; j++) {
|
||||
if (i === j) a[k][j] = 0;
|
||||
else a[k][j] += c * a[i][j];
|
||||
if (i === j) a[k][j] = 0; else a[k][j] += c * a[i][j];
|
||||
}
|
||||
b[k] += c * b[i];
|
||||
}
|
||||
}
|
||||
|
||||
const x = new Array(n).fill(0);
|
||||
for (let i = n - 1; i >= 0; i--) {
|
||||
if (Math.abs(a[i][i]) < 1e-12) continue;
|
||||
x[i] = b[i] / a[i][i];
|
||||
for (let k = i - 1; k >= 0; k--) {
|
||||
b[k] -= a[k][i] * x[i];
|
||||
}
|
||||
for (let k = i - 1; k >= 0; k--) b[k] -= a[k][i] * x[i];
|
||||
}
|
||||
return x;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate RMSE between predicted and actual values.
|
||||
*/
|
||||
export function calcRMSE(actual: number[], predicted: number[]): number {
|
||||
const n = Math.min(actual.length, predicted.length);
|
||||
if (n === 0) return Infinity;
|
||||
const sumSq = actual.reduce((sum, a, i) => {
|
||||
if (i >= predicted.length) return sum;
|
||||
return sum + (a - predicted[i]) ** 2;
|
||||
}, 0);
|
||||
return Math.sqrt(sumSq / n);
|
||||
return Math.sqrt(actual.reduce((sum, a, i) => i < predicted.length ? sum + (a - predicted[i]) ** 2 : sum, 0) / n);
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate R² coefficient of determination.
|
||||
*/
|
||||
export function calcRSquared(actual: number[], predicted: number[]): number {
|
||||
const n = Math.min(actual.length, predicted.length);
|
||||
if (n < 2) return 0;
|
||||
const meanActual = mean(actual);
|
||||
const ssTot = actual.reduce((sum, y) => sum + (y - meanActual) ** 2, 0);
|
||||
if (ssTot === 0) return 1;
|
||||
const ssRes = actual.reduce((sum, y, i) => {
|
||||
if (i >= predicted.length) return sum;
|
||||
return sum + (y - predicted[i]) ** 2;
|
||||
}, 0);
|
||||
return 1 - ssRes / ssTot;
|
||||
return 1 - actual.reduce((sum, y, i) => i < predicted.length ? sum + (y - predicted[i]) ** 2 : sum, 0) / ssTot;
|
||||
}
|
||||
|
||||
/**
|
||||
* Median Absolute Deviation outlier detection.
|
||||
* Returns indices of inlier samples.
|
||||
*/
|
||||
export function findInliersMAD(
|
||||
values: number[],
|
||||
threshold: number = 3.5
|
||||
): number[] {
|
||||
export function findInliersMAD(values: number[], threshold: number = 3.5): number[] {
|
||||
const med = median(values);
|
||||
const absDevs = values.map(v => Math.abs(v - med));
|
||||
const mad = median(absDevs);
|
||||
const mad = median(values.map(v => Math.abs(v - med)));
|
||||
if (mad === 0) return values.map((_, i) => i);
|
||||
|
||||
return values
|
||||
.map((v, i) => ({ v, i, modifiedZ: 0.6745 * Math.abs(v - med) / mad }))
|
||||
.filter(x => x.modifiedZ < threshold)
|
||||
.map(x => x.i);
|
||||
return values.map((v, i) => ({ v, i, z: 0.6745 * Math.abs(v - med) / mad }))
|
||||
.filter(x => x.z < threshold).map(x => x.i);
|
||||
}
|
||||
|
||||
/**
|
||||
* K-Means clustering (for PPK-based tier detection).
|
||||
* Returns cluster assignments (0..k-1) for each sample.
|
||||
* Compute K-Means inertia (sum of squared distances from each point to its centroid).
|
||||
* Lower inertia = better clustering.
|
||||
*/
|
||||
function calcInertia(values: number[], assignments: number[], centroids: number[]): number {
|
||||
return values.reduce((sum, v, i) => sum + (v - centroids[assignments[i]]) ** 2, 0);
|
||||
}
|
||||
|
||||
/**
|
||||
* Single K-Means run. Returns assignments, centroids, and inertia.
|
||||
*/
|
||||
function kMeansOnce(
|
||||
values: number[],
|
||||
k: number,
|
||||
maxIterations: number
|
||||
): { assignments: number[]; centroids: number[]; inertia: number } {
|
||||
// K-Means++ seeding for better initialisation
|
||||
const centroids: number[] = [];
|
||||
centroids.push(values[Math.floor(Math.random() * values.length)]);
|
||||
for (let c = 1; c < k; c++) {
|
||||
const dists = values.map(v => Math.min(...centroids.map(cent => (v - cent) ** 2)));
|
||||
const total = dists.reduce((a, b) => a + b, 0);
|
||||
let r = Math.random() * total;
|
||||
for (let i = 0; i < dists.length; i++) {
|
||||
r -= dists[i];
|
||||
if (r <= 0) { centroids.push(values[i]); break; }
|
||||
}
|
||||
// Fallback: if loop exits without pushing (floating point edge case)
|
||||
if (centroids.length < c + 1) centroids.push(values[values.length - 1]);
|
||||
}
|
||||
|
||||
const assignments = new Array(values.length).fill(0);
|
||||
for (let iter = 0; iter < maxIterations; iter++) {
|
||||
let changed = false;
|
||||
for (let i = 0; i < values.length; i++) {
|
||||
let minDist = Infinity, best = 0;
|
||||
for (let c = 0; c < k; c++) {
|
||||
const dist = Math.abs(values[i] - centroids[c]);
|
||||
if (dist < minDist) { minDist = dist; best = c; }
|
||||
}
|
||||
if (assignments[i] !== best) { assignments[i] = best; changed = true; }
|
||||
}
|
||||
if (!changed) break;
|
||||
for (let c = 0; c < k; c++) {
|
||||
const clusterVals = values.filter((_, i) => assignments[i] === c);
|
||||
if (clusterVals.length > 0) centroids[c] = mean(clusterVals);
|
||||
}
|
||||
}
|
||||
|
||||
const inertia = calcInertia(values, assignments, centroids);
|
||||
return { assignments, centroids, inertia };
|
||||
}
|
||||
|
||||
/**
|
||||
* K-Means clustering with multiple restarts.
|
||||
* Runs `runs` times and returns the assignment with the lowest inertia,
|
||||
* eliminating randomness instability across different executions.
|
||||
*/
|
||||
export function kMeans(
|
||||
values: number[],
|
||||
k: number,
|
||||
maxIterations: number = 100
|
||||
maxIterations: number = 100,
|
||||
runs: number = 8
|
||||
): number[] {
|
||||
if (values.length < k) return values.map(() => 0);
|
||||
|
||||
// Initialize centroids using k-means++
|
||||
let centroids: number[] = [];
|
||||
centroids.push(values[Math.floor(Math.random() * values.length)]);
|
||||
for (let c = 1; c < k; c++) {
|
||||
const dists = values.map(v => Math.min(
|
||||
...centroids.map(cent => Math.abs(v - cent))
|
||||
));
|
||||
const totalDist = dists.reduce((a, b) => a + b, 0);
|
||||
let r = Math.random() * totalDist;
|
||||
for (let i = 0; i < dists.length; i++) {
|
||||
r -= dists[i];
|
||||
if (r <= 0) {
|
||||
centroids.push(values[i]);
|
||||
break;
|
||||
}
|
||||
let bestResult: { assignments: number[]; centroids: number[]; inertia: number } | null = null;
|
||||
|
||||
for (let run = 0; run < runs; run++) {
|
||||
const result = kMeansOnce(values, k, maxIterations);
|
||||
if (bestResult === null || result.inertia < bestResult.inertia) {
|
||||
bestResult = result;
|
||||
}
|
||||
}
|
||||
|
||||
const assignments = new Array(values.length).fill(0);
|
||||
|
||||
for (let iter = 0; iter < maxIterations; iter++) {
|
||||
// Assign
|
||||
let changed = false;
|
||||
for (let i = 0; i < values.length; i++) {
|
||||
let minDist = Infinity;
|
||||
let bestCluster = 0;
|
||||
for (let c = 0; c < k; c++) {
|
||||
const dist = Math.abs(values[i] - centroids[c]);
|
||||
if (dist < minDist) {
|
||||
minDist = dist;
|
||||
bestCluster = c;
|
||||
}
|
||||
}
|
||||
if (assignments[i] !== bestCluster) {
|
||||
assignments[i] = bestCluster;
|
||||
changed = true;
|
||||
}
|
||||
}
|
||||
|
||||
if (!changed) break;
|
||||
|
||||
// Update centroids
|
||||
for (let c = 0; c < k; c++) {
|
||||
const clusterVals = values.filter((_, i) => assignments[i] === c);
|
||||
if (clusterVals.length > 0) {
|
||||
centroids[c] = mean(clusterVals);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Sort clusters by centroid value (ascending: economy < standard < premium)
|
||||
const centroidOrder = centroids
|
||||
.map((c, i) => ({ centroid: c, index: i }))
|
||||
.sort((a, b) => a.centroid - b.centroid);
|
||||
|
||||
const labelMap = new Map<number, number>();
|
||||
centroidOrder.forEach((item, newIdx) => labelMap.set(item.index, newIdx));
|
||||
|
||||
return assignments.map(a => labelMap.get(a)!);
|
||||
// Re-order cluster indices so label 0 = lowest centroid (economy), etc.
|
||||
const { assignments, centroids } = bestResult!;
|
||||
const order = centroids.map((c, i) => ({ c, i })).sort((a, b) => a.c - b.c);
|
||||
const map = new Map(order.map((item, idx) => [item.i, idx]));
|
||||
return assignments.map(a => map.get(a)!);
|
||||
}
|
||||
|
||||
/**
|
||||
* Find "knee point" in price-vs-distance curve for minimum fare detection.
|
||||
* Uses simple piecewise linear fit.
|
||||
*/
|
||||
export function detectMinimumFare(
|
||||
distances: number[],
|
||||
prices: number[],
|
||||
kmRate: number
|
||||
): number | null {
|
||||
export function detectMinimumFare(distances: number[], prices: number[], kmRate: number): number | null {
|
||||
if (distances.length < 5) return null;
|
||||
|
||||
// Sort by distance
|
||||
const pairs = distances.map((d, i) => ({ d, p: prices[i] }))
|
||||
.sort((a, b) => a.d - b.d);
|
||||
|
||||
// Compute expected price without min fare
|
||||
const residuals = pairs.map(({ d, p }) => p - kmRate * d);
|
||||
|
||||
// Find where actual price consistently exceeds predicted
|
||||
// The minimum fare is the max of (price - kmRate*dist) for short rides
|
||||
const shortRides = pairs.filter(({ d }) => d < 10);
|
||||
if (shortRides.length < 3) return null;
|
||||
|
||||
const minFareEstimate = Math.max(
|
||||
...shortRides.map(({ d, p }) => p - kmRate * d)
|
||||
);
|
||||
|
||||
return minFareEstimate > 0 ? Math.round(minFareEstimate * 100) / 100 : null;
|
||||
const pairs = distances.map((d, i) => ({ d, p: prices[i] })).sort((a, b) => a.d - b.d);
|
||||
const shortCount = Math.max(5, Math.floor(pairs.length * 0.3));
|
||||
const shortRides = pairs.slice(0, shortCount);
|
||||
const residuals = shortRides.map(({ d, p }) => p - kmRate * d).sort((a, b) => a - b);
|
||||
const trimIdx = Math.max(0, Math.floor(residuals.length * 0.1));
|
||||
const trimmed = residuals.slice(trimIdx, residuals.length - trimIdx);
|
||||
const estimate = trimmed.length > 0 ? Math.max(...trimmed) : 0;
|
||||
return estimate > 0 ? Math.round(estimate * 100) / 100 : null;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user