Update: 2026-07-06 17:00:43
This commit is contained in:
Hamza-Ayed
@@ -0,0 +1,326 @@
|
||||
/**
|
||||
* Matrix and statistical utilities for pricing analysis.
|
||||
* Pure math — no external dependencies except 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;
|
||||
const mx = x.reduce((a, b) => a + b, 0) / n;
|
||||
const my = y.reduce((a, b) => a + b, 0) / n;
|
||||
let num = 0, dx2 = 0, dy2 = 0;
|
||||
for (let i = 0; i < n; i++) {
|
||||
const dx = x[i] - mx;
|
||||
const dy = y[i] - my;
|
||||
num += dx * dy;
|
||||
dx2 += dx * dx;
|
||||
dy2 += dy * dy;
|
||||
}
|
||||
const denom = Math.sqrt(dx2 * dy2);
|
||||
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 corr = pearsonCorr(dists, durs);
|
||||
|
||||
const lambda = Math.abs(corr) > 0.85 ? 0.5 : 0.01; // ridge penalty
|
||||
|
||||
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;
|
||||
|
||||
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],
|
||||
];
|
||||
|
||||
const B = [sumY, sumX1Y, sumX2Y];
|
||||
|
||||
try {
|
||||
const beta = gaussianElimination(A, B);
|
||||
const 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;
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
|
||||
return { baseFare, kmRate, minRate };
|
||||
} catch {
|
||||
return null;
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* 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
|
||||
): { baseFare: number; kmRate: number; minRate: number } | null {
|
||||
let currentSamples = [...samples];
|
||||
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);
|
||||
|
||||
for (let i = 0; i < maxIterations; i++) {
|
||||
const predicted = currentSamples.map(
|
||||
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;
|
||||
});
|
||||
|
||||
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;
|
||||
for (let k = i + 1; k < n; 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];
|
||||
}
|
||||
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];
|
||||
}
|
||||
}
|
||||
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);
|
||||
}
|
||||
|
||||
/**
|
||||
* 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;
|
||||
}
|
||||
|
||||
/**
|
||||
* Median Absolute Deviation outlier detection.
|
||||
* Returns indices of inlier samples.
|
||||
*/
|
||||
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);
|
||||
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);
|
||||
}
|
||||
|
||||
/**
|
||||
* K-Means clustering (for PPK-based tier detection).
|
||||
* Returns cluster assignments (0..k-1) for each sample.
|
||||
*/
|
||||
export function kMeans(
|
||||
values: number[],
|
||||
k: number,
|
||||
maxIterations: number = 100
|
||||
): 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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)!);
|
||||
}
|
||||
|
||||
/**
|
||||
* 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 {
|
||||
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;
|
||||
}
|
||||
Reference in New Issue
Block a user