Update: 2026-07-06 18:23:46
This commit is contained in:
Hamza-Ayed
@@ -1,5 +1,9 @@
|
||||
import { median, mean } from 'simple-statistics';
|
||||
|
||||
// ─────────────────────────────────────────────
|
||||
// Internal helpers
|
||||
// ─────────────────────────────────────────────
|
||||
|
||||
function pearsonCorr(x: number[], y: number[]): number {
|
||||
const n = Math.min(x.length, y.length);
|
||||
if (n < 3) return 0;
|
||||
@@ -17,90 +21,6 @@ function pearsonCorr(x: number[], y: number[]): number {
|
||||
return denom === 0 ? 0 : num / denom;
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
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;
|
||||
|
||||
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, 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;
|
||||
}
|
||||
|
||||
const A = [
|
||||
[n, sumX1, sumX2],
|
||||
[sumX1, sumX1Sq + lambda, sumX1X2],
|
||||
[sumX2, sumX1X2, sumX2Sq + lambda],
|
||||
];
|
||||
const B = [sumY, sumX1Y, sumX2Y];
|
||||
|
||||
try {
|
||||
const beta = gaussianElimination(A, B);
|
||||
let baseFare = Math.max(0, beta[0]);
|
||||
let kmRate = Math.max(0, beta[1]);
|
||||
let minRate = Math.max(0, beta[2]);
|
||||
|
||||
const predFull = samples.map(s => baseFare + kmRate * s.distance_km + minRate * s.duration_min);
|
||||
const rmseFull = calcRMSE(prices, predFull);
|
||||
|
||||
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 };
|
||||
} catch {
|
||||
return null;
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
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((_, 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;
|
||||
}
|
||||
|
||||
function gaussianElimination(A: number[][], B: number[]): number[] {
|
||||
const n = A.length;
|
||||
const a = A.map(row => [...row]);
|
||||
@@ -130,10 +50,208 @@ function gaussianElimination(A: number[][], B: number[]): number[] {
|
||||
return x;
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────
|
||||
// Core regression — full 3-parameter model
|
||||
// price = baseFare + kmRate × dist + minRate × dur
|
||||
// ─────────────────────────────────────────────
|
||||
|
||||
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;
|
||||
|
||||
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);
|
||||
// Stronger ridge when predictors are collinear (typical in taxi data)
|
||||
const lambda = Math.abs(corr) > 0.85 ? 1.5 : 0.05;
|
||||
|
||||
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, 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;
|
||||
}
|
||||
|
||||
const A = [
|
||||
[n, sumX1, sumX2 ],
|
||||
[sumX1, sumX1Sq + lambda, sumX1X2 ],
|
||||
[sumX2, sumX1X2, sumX2Sq + lambda],
|
||||
];
|
||||
const B = [sumY, sumX1Y, sumX2Y];
|
||||
|
||||
try {
|
||||
const beta = gaussianElimination(A, B);
|
||||
return {
|
||||
baseFare: Math.max(0, beta[0]),
|
||||
kmRate: Math.max(0, beta[1]),
|
||||
minRate: Math.max(0, beta[2]),
|
||||
};
|
||||
} catch {
|
||||
return null;
|
||||
}
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────
|
||||
// Two-Stage Regression
|
||||
// Stage 1: estimate flag fall from shortest rides
|
||||
// Stage 2: regress residuals on (dist, dur) with no intercept
|
||||
// ─────────────────────────────────────────────
|
||||
|
||||
/**
|
||||
* Stage 1 — Estimate the flag fall (فتحة العداد / meter opening charge).
|
||||
*
|
||||
* Takes the shortest 20% of rides by distance (min 5 samples) and fits
|
||||
* a simple linear model: price ~ intercept + slope × dist.
|
||||
* The intercept is the flag fall estimate.
|
||||
*
|
||||
* Clamped to [0, 65% of median price] to avoid unreasonable values.
|
||||
*/
|
||||
export function estimateFlagFall(
|
||||
samples: Array<{ distance_km: number; duration_min: number; price: number }>
|
||||
): number {
|
||||
if (samples.length < 5) return 0;
|
||||
|
||||
const sorted = [...samples].sort((a, b) => a.distance_km - b.distance_km);
|
||||
const shortCount = Math.max(5, Math.floor(sorted.length * 0.20));
|
||||
const shortRides = sorted.slice(0, shortCount);
|
||||
|
||||
// Simple OLS: price ~ a + b × dist on short rides only
|
||||
const n = shortRides.length;
|
||||
const sumX = shortRides.reduce((s, r) => s + r.distance_km, 0);
|
||||
const sumY = shortRides.reduce((s, r) => s + r.price, 0);
|
||||
const sumXX = shortRides.reduce((s, r) => s + r.distance_km ** 2, 0);
|
||||
const sumXY = shortRides.reduce((s, r) => s + r.distance_km * r.price, 0);
|
||||
|
||||
const denom = n * sumXX - sumX * sumX;
|
||||
if (Math.abs(denom) < 1e-10) {
|
||||
// Degenerate case — return a safe lower-bound estimate
|
||||
return Math.min(...shortRides.map(r => r.price)) * 0.4;
|
||||
}
|
||||
|
||||
const slope = (n * sumXY - sumX * sumY) / denom;
|
||||
const intercept = (sumY - slope * sumX) / n;
|
||||
|
||||
const medianPrice = median(samples.map(s => s.price));
|
||||
return Math.max(0, Math.min(intercept, medianPrice * 0.65));
|
||||
}
|
||||
|
||||
/**
|
||||
* Stage 2 — Two-variable regression with no intercept.
|
||||
* Fits: (price − fixedBase) ~ kmRate × dist + minRate × dur
|
||||
*
|
||||
* Uses ridge regularization (lambda = 2.0 when dist/dur are collinear)
|
||||
* to distribute the effect between km and min rather than collapsing to one.
|
||||
*/
|
||||
function twoVarNoIntercept(
|
||||
samples: Array<{ distance_km: number; duration_min: number; price: number }>,
|
||||
fixedBase: number
|
||||
): { kmRate: number; minRate: number } | null {
|
||||
const n = samples.length;
|
||||
if (n < 3) return null;
|
||||
|
||||
const dists = samples.map(s => s.distance_km);
|
||||
const durs = samples.map(s => s.duration_min);
|
||||
const corr = pearsonCorr(dists, durs);
|
||||
|
||||
// Higher ridge when predictors are correlated — forces balance between km and min
|
||||
const lambda = Math.abs(corr) > 0.85 ? 2.0 : 0.5;
|
||||
|
||||
let s11 = 0, s22 = 0, s12 = 0, s1y = 0, s2y = 0;
|
||||
for (let i = 0; i < n; i++) {
|
||||
const x1 = dists[i], x2 = durs[i];
|
||||
const y = samples[i].price - fixedBase;
|
||||
s11 += x1 * x1; s22 += x2 * x2; s12 += x1 * x2;
|
||||
s1y += x1 * y; s2y += x2 * y;
|
||||
}
|
||||
|
||||
// Solve 2×2 ridge system:
|
||||
// [ s11+λ s12 ] [ kmRate ] [ s1y ]
|
||||
// [ s12 s22+λ ] [ minRate ] = [ s2y ]
|
||||
const a = s11 + lambda, b = s12, d = s22 + lambda;
|
||||
const det = a * d - b * b;
|
||||
if (Math.abs(det) < 1e-12) return null;
|
||||
|
||||
return {
|
||||
kmRate: Math.max(0, (s1y * d - s2y * b) / det),
|
||||
minRate: Math.max(0, (a * s2y - b * s1y) / det),
|
||||
};
|
||||
}
|
||||
|
||||
/**
|
||||
* Two-Stage Regression (main entry point for the engine).
|
||||
*
|
||||
* Properly decomposes taxi pricing into three components:
|
||||
* price = baseFare (flag fall) + kmRate × dist + minRate × dur
|
||||
*
|
||||
* Stage 1 fixes baseFare from shortest rides.
|
||||
* Stage 2 fits kmRate and minRate on residuals.
|
||||
*
|
||||
* This avoids the distance/duration collinearity problem by removing
|
||||
* the constant component first.
|
||||
*/
|
||||
export function twoStageRegression(
|
||||
samples: Array<{ distance_km: number; duration_min: number; price: number }>
|
||||
): { baseFare: number; kmRate: number; minRate: number } | null {
|
||||
if (samples.length < 5) return null;
|
||||
|
||||
const baseFare = estimateFlagFall(samples);
|
||||
const rates = twoVarNoIntercept(samples, baseFare);
|
||||
if (!rates) return null;
|
||||
|
||||
return { baseFare, kmRate: rates.kmRate, minRate: rates.minRate };
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────
|
||||
// Robust regression (iterative outlier removal)
|
||||
// ─────────────────────────────────────────────
|
||||
|
||||
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;
|
||||
|
||||
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((_, 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;
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────
|
||||
// Statistics utilities
|
||||
// ─────────────────────────────────────────────
|
||||
|
||||
export function calcRMSE(actual: number[], predicted: number[]): number {
|
||||
const n = Math.min(actual.length, predicted.length);
|
||||
if (n === 0) return Infinity;
|
||||
return Math.sqrt(actual.reduce((sum, a, i) => i < predicted.length ? sum + (a - predicted[i]) ** 2 : sum, 0) / n);
|
||||
return Math.sqrt(
|
||||
actual.reduce((sum, a, i) => i < predicted.length ? sum + (a - predicted[i]) ** 2 : sum, 0) / n
|
||||
);
|
||||
}
|
||||
|
||||
export function calcRSquared(actual: number[], predicted: number[]): number {
|
||||
@@ -149,27 +267,26 @@ export function findInliersMAD(values: number[], threshold: number = 3.5): numbe
|
||||
const med = median(values);
|
||||
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, z: 0.6745 * Math.abs(v - med) / mad }))
|
||||
.filter(x => x.z < 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);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute K-Means inertia (sum of squared distances from each point to its centroid).
|
||||
* Lower inertia = better clustering.
|
||||
*/
|
||||
// ─────────────────────────────────────────────
|
||||
// K-Means clustering (multi-run for stability)
|
||||
// ─────────────────────────────────────────────
|
||||
|
||||
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
|
||||
// K-Means++ seeding
|
||||
const centroids: number[] = [];
|
||||
centroids.push(values[Math.floor(Math.random() * values.length)]);
|
||||
for (let c = 1; c < k; c++) {
|
||||
@@ -180,7 +297,6 @@ function kMeansOnce(
|
||||
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]);
|
||||
}
|
||||
|
||||
@@ -202,14 +318,12 @@ function kMeansOnce(
|
||||
}
|
||||
}
|
||||
|
||||
const inertia = calcInertia(values, assignments, centroids);
|
||||
return { assignments, centroids, inertia };
|
||||
return { assignments, centroids, inertia: calcInertia(values, assignments, centroids) };
|
||||
}
|
||||
|
||||
/**
|
||||
* K-Means clustering with multiple restarts.
|
||||
* Runs `runs` times and returns the assignment with the lowest inertia,
|
||||
* eliminating randomness instability across different executions.
|
||||
* K-Means with multiple restarts — picks the run with lowest inertia
|
||||
* to eliminate randomness instability across executions.
|
||||
*/
|
||||
export function kMeans(
|
||||
values: number[],
|
||||
@@ -219,30 +333,30 @@ export function kMeans(
|
||||
): number[] {
|
||||
if (values.length < k) return values.map(() => 0);
|
||||
|
||||
let bestResult: { assignments: number[]; centroids: number[]; inertia: number } | null = null;
|
||||
|
||||
let best: ReturnType<typeof kMeansOnce> | null = null;
|
||||
for (let run = 0; run < runs; run++) {
|
||||
const result = kMeansOnce(values, k, maxIterations);
|
||||
if (bestResult === null || result.inertia < bestResult.inertia) {
|
||||
bestResult = result;
|
||||
}
|
||||
if (!best || result.inertia < best.inertia) best = result;
|
||||
}
|
||||
|
||||
// Re-order cluster indices so label 0 = lowest centroid (economy), etc.
|
||||
const { assignments, centroids } = bestResult!;
|
||||
const { assignments, centroids } = best!;
|
||||
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]));
|
||||
const map = new Map(order.map((item, idx) => [item.i, idx]));
|
||||
return assignments.map(a => map.get(a)!);
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────
|
||||
// Minimum fare detection
|
||||
// ─────────────────────────────────────────────
|
||||
|
||||
export function detectMinimumFare(distances: number[], prices: number[], kmRate: number): number | null {
|
||||
if (distances.length < 5) return null;
|
||||
const pairs = distances.map((d, i) => ({ d, p: prices[i] })).sort((a, b) => a.d - b.d);
|
||||
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;
|
||||
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