diff --git a/CHANGELOG.md b/CHANGELOG.md index 4bdba941dd..da026b8a4e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ### Added +- Added `LINEST` for simple and multiple linear regression, with optional intercept and regression statistics. The `stats` argument must be constant because result dimensions are determined before evaluation. [#1769](https://github.com/handsontable/hyperformula/pull/1769) - Added support for the new license key format. A proprietary key can now grant a subset of the library: a function your key does not include evaluates to a `#LIC!` error, and the matching parts of the API throw a `LicenseCapabilityMissingError`. `getAvailableFunctions()` and `getFunctionDetails()` describe only the functions your key includes. [#1728](https://github.com/handsontable/hyperformula/pull/1728) ### Changed diff --git a/docs/guide/known-limitations.md b/docs/guide/known-limitations.md index d78cac018d..8fd335900a 100644 --- a/docs/guide/known-limitations.md +++ b/docs/guide/known-limitations.md @@ -61,6 +61,16 @@ a circular reference. * Ordering (including mixed types, empty cells, and text collation) follows HyperFormula's own comparison rules, which honor the `caseSensitive` and `accentSensitive` configuration options. Numbers sort before text, and text before logical values. +### LINEST function + +* The `stats` argument controls whether the result has one or five rows. It must be omitted or supplied as a constant: `TRUE()`, `FALSE()`, a number, or the text `"TRUE"` or `"FALSE"`. Parentheses and numeric unary signs are supported. Cell references, named expressions, and calculated expressions such as `1=1` return `#VALUE!`, including when `LINEST` is nested inside another function. + +* Result width is determined from the predictor dimensions before evaluation. Changes to values within fixed input ranges and to a referenced `const` argument recalculate normally. Changes to dependency values do not resize the result automatically; this is the existing array-sizing limitation. + +* If an input expression returns a different predictor count than predicted, `LINEST` returns `#VALUE!`. For example, filtering predictor columns can change the required result width. + +* Very large or very small input magnitudes can reduce numerical accuracy. Intermediate calculations, such as residual sums of squares, can still overflow or underflow, producing `#NUM!` or inaccurate statistics. Rescale inputs to more moderate units where possible. See [LINEST numerical differences](list-of-differences.md#linest) for compatibility considerations. + ### OFFSET function HyperFormula resolves the OFFSET function at parse time rather than during evaluation. The parser inspects the arguments and rewrites the expression into a plain cell reference or range. This keeps the dependency graph accurate but imposes several restrictions. diff --git a/docs/guide/list-of-differences.md b/docs/guide/list-of-differences.md index a109f82616..e4a3b28fdc 100644 --- a/docs/guide/list-of-differences.md +++ b/docs/guide/list-of-differences.md @@ -34,6 +34,14 @@ See a full list of differences between HyperFormula, Microsoft Excel, and Google **Contents:** [[toc]] +## LINEST + +Unlike Excel, HyperFormula requires a constant `stats` argument because the output dimensions are determined before evaluation. For example, `=LINEST(A1:A10,B1:B10,TRUE(),D1)` returns `#VALUE!`; use `TRUE()` or `FALSE()` directly for `stats`. See [LINEST limitations](known-limitations.md#linest-function) for supported constants and input-sizing requirements. + +Numerical results can differ for nearly dependent predictors and nearly perfect fits. Coefficient standard errors are sensitive to conditioning, and the F statistic is sensitive to residuals close to machine precision. For an effectively perfect multiple regression, Excel and HyperFormula can return different large finite F values even when the coefficients agree. + +Very large or very small input scales can also cause substantial differences in coefficients and statistics, including for a single predictor with a non-perfect fit. HyperFormula does not reproduce Excel's loss of predictors or zero standard errors observed at extreme scales. Rescaling inputs to more moderate units can reduce numerical errors in both engines. + ## General functionalities | Functionality | Examples | HyperFormula | Google Sheets | Microsoft Excel | diff --git a/src/error-message.ts b/src/error-message.ts index d42040901c..9f0ef381db 100644 --- a/src/error-message.ts +++ b/src/error-message.ts @@ -7,6 +7,8 @@ * This is a class for detailed error messages across HyperFormula. */ export class ErrorMessage { + public static LinestStaticSize = 'LINEST requires input dimensions that determine a fixed result size.' + public static LinestStaticStats = 'LINEST requires a constant stats argument to determine its result size.' public static DistinctSigns = 'Distinct signs.' public static WrongArgNumber = 'Wrong number of arguments.' public static EmptyArg = 'Empty function argument.' diff --git a/src/i18n/languages/csCZ.ts b/src/i18n/languages/csCZ.ts index fea286083c..d3e0575403 100644 --- a/src/i18n/languages/csCZ.ts +++ b/src/i18n/languages/csCZ.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'JE.TEXT', LEFT: 'ZLEVA', LEN: 'DÉLKA', + LINEST: 'LINREGRESE', LN: 'LN', LOG10: 'LOG', LOG: 'LOGZ', diff --git a/src/i18n/languages/daDK.ts b/src/i18n/languages/daDK.ts index 9264551249..f3737cefae 100644 --- a/src/i18n/languages/daDK.ts +++ b/src/i18n/languages/daDK.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ER.TEKST', LEFT: 'VENSTRE', LEN: 'LÆNGDE', + LINEST: 'LINREGR', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/deDE.ts b/src/i18n/languages/deDE.ts index 4d3d39319b..0b1a180f0b 100644 --- a/src/i18n/languages/deDE.ts +++ b/src/i18n/languages/deDE.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ISTTEXT', LEFT: 'LINKS', LEN: 'LÄNGE', + LINEST: 'RGP', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/enGB.ts b/src/i18n/languages/enGB.ts index d271b732cd..3a7185fdad 100644 --- a/src/i18n/languages/enGB.ts +++ b/src/i18n/languages/enGB.ts @@ -144,6 +144,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ISTEXT', LEFT: 'LEFT', LEN: 'LEN', + LINEST: 'LINEST', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/esES.ts b/src/i18n/languages/esES.ts index 9c6b817b98..df2038ff36 100644 --- a/src/i18n/languages/esES.ts +++ b/src/i18n/languages/esES.ts @@ -143,6 +143,7 @@ export const dictionary: RawTranslationPackage = { ISTEXT: 'ESTEXTO', LEFT: 'IZQUIERDA', LEN: 'LARGO', + LINEST: 'ESTIMACION.LINEAL', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/fiFI.ts b/src/i18n/languages/fiFI.ts index 99ff05708f..87ab6be3cd 100644 --- a/src/i18n/languages/fiFI.ts +++ b/src/i18n/languages/fiFI.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ONTEKSTI', LEFT: 'VASEN', LEN: 'PITUUS', + LINEST: 'LINREGR', LN: 'LUONNLOG', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/frFR.ts b/src/i18n/languages/frFR.ts index f7b9be45f7..bc962747b3 100644 --- a/src/i18n/languages/frFR.ts +++ b/src/i18n/languages/frFR.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ESTTEXTE', LEFT: 'GAUCHE', LEN: 'NBCAR', + LINEST: 'DROITEREG', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/huHU.ts b/src/i18n/languages/huHU.ts index e40219553a..245709fa3b 100644 --- a/src/i18n/languages/huHU.ts +++ b/src/i18n/languages/huHU.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'SZÖVEG.E', LEFT: 'BAL', LEN: 'HOSSZ', + LINEST: 'LIN.ILL', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/idID.ts b/src/i18n/languages/idID.ts index 719cb11db6..532ee25e31 100644 --- a/src/i18n/languages/idID.ts +++ b/src/i18n/languages/idID.ts @@ -144,6 +144,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ADALAH.TEKS', LEFT: 'KIRI', LEN: 'PANJANG', + LINEST: 'LINEST', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/itIT.ts b/src/i18n/languages/itIT.ts index 1a494628f8..1a863ea46a 100644 --- a/src/i18n/languages/itIT.ts +++ b/src/i18n/languages/itIT.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'VAL.TESTO', LEFT: 'SINISTRA', LEN: 'LUNGHEZZA', + LINEST: 'REGR.LIN', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/nbNO.ts b/src/i18n/languages/nbNO.ts index ff15415a98..de4a1cf667 100644 --- a/src/i18n/languages/nbNO.ts +++ b/src/i18n/languages/nbNO.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ERTEKST', LEFT: 'VENSTRE', LEN: 'LENGDE', + LINEST: 'RETTLINJE', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/nlNL.ts b/src/i18n/languages/nlNL.ts index c491e2f44b..a7999b8913 100644 --- a/src/i18n/languages/nlNL.ts +++ b/src/i18n/languages/nlNL.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ISTEKST', LEFT: 'LINKS', LEN: 'PITUUS', + LINEST: 'LIJNSCH', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/plPL.ts b/src/i18n/languages/plPL.ts index 6fa4f016b4..7548158e88 100644 --- a/src/i18n/languages/plPL.ts +++ b/src/i18n/languages/plPL.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'CZY.TEKST', LEFT: 'LEWY', LEN: 'DŁ', + LINEST: 'REGLINP', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/ptPT.ts b/src/i18n/languages/ptPT.ts index 6c1e68c328..a74fba28b5 100644 --- a/src/i18n/languages/ptPT.ts +++ b/src/i18n/languages/ptPT.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ÉTEXTO', LEFT: 'ESQUERDA', LEN: 'NÚM.CARACT', + LINEST: 'PROJ.LIN', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/ruRU.ts b/src/i18n/languages/ruRU.ts index 6b35c1ce2f..e808973691 100644 --- a/src/i18n/languages/ruRU.ts +++ b/src/i18n/languages/ruRU.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ЕТЕКСТ', LEFT: 'ЛЕВСИМВ', LEN: 'ДЛСТР', + LINEST: 'ЛИНЕЙН', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/svSE.ts b/src/i18n/languages/svSE.ts index 5f87580087..7d869a89f1 100644 --- a/src/i18n/languages/svSE.ts +++ b/src/i18n/languages/svSE.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ÄRTEXT', LEFT: 'VÄNSTER', LEN: 'LÄNGD', + LINEST: 'REGR', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/trTR.ts b/src/i18n/languages/trTR.ts index 93c2176d37..991b50a31f 100644 --- a/src/i18n/languages/trTR.ts +++ b/src/i18n/languages/trTR.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'EMETİNSE', LEFT: 'SOL', LEN: 'UZUNLUK', + LINEST: 'DOT', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/interpreter/functionMetadata/categories/statistical.ts b/src/interpreter/functionMetadata/categories/statistical.ts index 1752c3eaf6..6a404977f5 100644 --- a/src/interpreter/functionMetadata/categories/statistical.ts +++ b/src/interpreter/functionMetadata/categories/statistical.ts @@ -304,6 +304,18 @@ export const STATISTICAL_DOCS: Record = { documentationUrl: 'https://hyperformula.handsontable.com/docs/guide/built-in-functions.html', examples: ['=LARGE(A1:A10, 1)', '=LARGE(A1:A10, 3)'], }, + LINEST: { + category: 'Statistical', + shortDescription: 'Returns linear regression coefficients and optional statistics.', + parameters: [ + {name: 'known_y', description: 'A numeric range of observed dependent values.'}, + {name: 'known_x', description: 'Optional numeric predictors. If omitted, uses sequential values starting at 1 with the shape of known_y.'}, + {name: 'const', description: 'Whether to fit an intercept. Defaults to TRUE; FALSE fits through zero.'}, + {name: 'stats', description: 'Whether to return five rows including regression statistics. Defaults to FALSE. Must be a constant; cell references and computed expressions are unsupported.'}, + ], + documentationUrl: 'https://hyperformula.handsontable.com/docs/guide/built-in-functions.html', + examples: ['=LINEST(A1:A10, B1:C10)', '=LINEST(A1:A10, B1:C10, TRUE(), TRUE())'], + }, 'LOGNORM.DIST': { category: 'Statistical', shortDescription: 'Returns density of lognormal distribution.', diff --git a/src/interpreter/plugin/RegressionPlugin.ts b/src/interpreter/plugin/RegressionPlugin.ts new file mode 100644 index 0000000000..9df3160c1e --- /dev/null +++ b/src/interpreter/plugin/RegressionPlugin.ts @@ -0,0 +1,216 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +import {ArraySize} from '../../ArraySize' +import {CellError, ErrorType} from '../../Cell' +import {ErrorMessage} from '../../error-message' +import {Ast, AstNodeType, ProcedureAst} from '../../parser' +import {SimpleRangeValue} from '../../SimpleRangeValue' +import {coerceScalarToBoolean, coerceToRange} from '../ArithmeticHelper' +import {InterpreterState} from '../InterpreterState' +import {EmptyValue, getRawValue, InternalScalarValue, InterpreterValue, isExtendedNumber} from '../InterpreterValue' +import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' +import {fitLinearRegression, LinearRegressionResult} from './regression/LinearRegression' + +/** The orientation and predictor count shared by size prediction and runtime validation. */ +interface RegressionShape { + predictorCount: number, + orientation: 'paired' | 'columns' | 'rows', +} + +/** Classifies the original dimensions before any values are flattened. */ +function regressionShape(y: ArraySize, x?: ArraySize): RegressionShape | undefined { + if (x === undefined || (y.width === x.width && y.height === x.height)) { + return {predictorCount: 1, orientation: 'paired'} + } + if (y.width === 1 && y.height === x.height) { + return {predictorCount: x.width, orientation: 'columns'} + } + if (y.height === 1 && y.width === x.width) { + return {predictorCount: x.height, orientation: 'rows'} + } + return undefined +} + +/** + * Converts validated, flattened predictor values into one row per observation. + * In columns orientation, each source row is an observation; in rows orientation, + * each source row is a predictor. Paired ranges supply one predictor per observation. + */ +function buildPredictorRows(predictorValues: number[], observationCount: number, shape: RegressionShape): number[][] { + return Array.from({length: observationCount}, (_, observationIndex) => { + if (shape.orientation === 'paired') { + return [predictorValues[observationIndex]] + } + return Array.from({length: shape.predictorCount}, (_, predictorIndex) => { + let flatIndex: number + if (shape.orientation === 'columns') { + flatIndex = observationIndex * shape.predictorCount + predictorIndex + } else { + flatIndex = predictorIndex * observationCount + observationIndex + } + return predictorValues[flatIndex] + }) + }) +} + +/** LINEST rejects empty strings as Boolean options, including formula-generated strings. */ +function regressionBoolean(value: InternalScalarValue): boolean | CellError | undefined { + return value === '' ? undefined : coerceScalarToBoolean(value) +} + +/** Recognizes a numeric literal with optional parentheses and unary signs. */ +function isNumericConstant(ast: Ast): boolean { + if (ast.type === AstNodeType.PARENTHESIS) { + return isNumericConstant(ast.expression) + } + if (ast.type === AstNodeType.PLUS_UNARY_OP || ast.type === AstNodeType.MINUS_UNARY_OP) { + return isNumericConstant(ast.value) + } + return ast.type === AstNodeType.NUMBER +} + +/** Resolves scalar constants without evaluating dependencies during array-size prediction. */ +function staticBoolean(ast: Ast | undefined): boolean | undefined { + if (ast === undefined || ast.type === AstNodeType.EMPTY) { + return false + } + if (ast.type === AstNodeType.PARENTHESIS) { + return staticBoolean(ast.expression) + } + if (ast.type === AstNodeType.NUMBER || ast.type === AstNodeType.STRING) { + const value = regressionBoolean(ast.value) + return typeof value === 'boolean' ? value : undefined + } + if (ast.type === AstNodeType.PLUS_UNARY_OP || ast.type === AstNodeType.MINUS_UNARY_OP) { + return isNumericConstant(ast.value) ? staticBoolean(ast.value) : undefined + } + if (ast.type === AstNodeType.FUNCTION_CALL && ast.args.length === 0) { + if (ast.procedureName === 'TRUE') { + return true + } + if (ast.procedureName === 'FALSE') { + return false + } + } + return undefined +} + +/** Converts numerical failures to spreadsheet errors, including inside an array result. */ +function finiteValue(value: number): number | CellError { + return Number.isFinite(value) ? value : new CellError(ErrorType.NUM, ErrorMessage.NaN) +} + +/** Formats the five-row result while preserving statistics errors independently of coefficients. */ +function regressionOutput(fit: LinearRegressionResult, statistics: boolean, fitIntercept: boolean): SimpleRangeValue { + const result: InternalScalarValue[][] = [[...fit.coefficients].reverse().concat(fit.intercept).map(finiteValue)] + if (!statistics) { + return SimpleRangeValue.onlyValues(result) + } + const regressionSumSquares = fit.totalSumSquares - fit.residualSumSquares + const variance = fit.degreesOfFreedom === 0 ? 0 : fit.residualSumSquares / fit.degreesOfFreedom + const rSquared = fit.totalSumSquares === 0 ? 1 : regressionSumSquares / fit.totalSumSquares + const f = finiteValue((regressionSumSquares / fit.retainedPredictorCount) / variance) + const interceptError = fitIntercept ? finiteValue(fit.interceptError) : new CellError(ErrorType.NA) + result.push([...fit.standardErrors].reverse().map(finiteValue).concat(interceptError)) + result.push([finiteValue(rSquared), finiteValue(Math.sqrt(variance))]) + result.push([f, fit.degreesOfFreedom]) + result.push([finiteValue(regressionSumSquares), finiteValue(fit.residualSumSquares)]) + for (const row of result) { + while (row.length < fit.coefficients.length + 1) { + row.push(new CellError(ErrorType.NA)) + } + } + return SimpleRangeValue.onlyValues(result) +} + +/** Implements linear regression with statically predictable result dimensions. */ +export class RegressionPlugin extends FunctionPlugin implements FunctionPluginTypecheck { + public static implementedFunctions: ImplementedFunctions = { + 'LINEST': { + method: 'linest', + sizeOfResultArrayMethod: 'linestArraySize', + vectorizationForbidden: true, + enableArrayArithmeticForArguments: true, + parameters: [ + {argumentType: FunctionArgumentType.RANGE}, + {argumentType: FunctionArgumentType.ANY, defaultValue: EmptyValue}, + {argumentType: FunctionArgumentType.SCALAR, defaultValue: true, emptyAsDefault: true}, + {argumentType: FunctionArgumentType.SCALAR, defaultValue: false, emptyAsDefault: true}, + ], + }, + } + + /** + * Evaluates LINEST(known_y, [known_x], [const], [stats]). + * Returns coefficients in reverse predictor order, followed by the intercept. + * With stats enabled, adds standard errors; R-squared and residual standard error; + * F and residual degrees of freedom; regression and residual sums of squares. + * Unused statistics cells and the intercept standard error for const=FALSE are #N/A. + * + * A dynamic stats option is rejected here as well as during size prediction so that + * nested consumers cannot accidentally bypass the fixed-size contract. + */ + public linest(ast: ProcedureAst, state: InterpreterState): InterpreterValue { + const generatedX = ast.args.length < 2 || ast.args[1].type === AstNodeType.EMPTY + return this.runFunction(ast.args, state, this.metadata('LINEST'), + (knownY: SimpleRangeValue, knownX: InterpreterValue, constArg: InternalScalarValue, statsArg: InternalScalarValue) => { + const fitIntercept = regressionBoolean(constArg) + const statistics = regressionBoolean(statsArg) + if (fitIntercept instanceof CellError) { + return fitIntercept + } + if (statistics instanceof CellError) { + return statistics + } + if (fitIntercept === undefined || statistics === undefined) { + return new CellError(ErrorType.VALUE, ErrorMessage.WrongType) + } + if (staticBoolean(ast.args[3]) === undefined) { + return new CellError(ErrorType.VALUE, ErrorMessage.LinestStaticStats) + } + const x = generatedX ? undefined : coerceToRange(knownX) + const shape = regressionShape(knownY.size, x?.size) + if (shape === undefined) { + return new CellError(ErrorType.REF, ErrorMessage.ArrayDimensions) + } + if (shape.predictorCount + 1 > this.config.maxColumns) { + return new CellError(ErrorType.VALUE, ErrorMessage.ValueLarge) + } + if (this.linestArraySize(ast, state).width !== shape.predictorCount + 1) { + return new CellError(ErrorType.VALUE, ErrorMessage.LinestStaticSize) + } + const yValues = Array.from(knownY.valuesFromTopLeftCorner()) + const xValues = x === undefined ? yValues.map((_, i) => i + 1) : Array.from(x.valuesFromTopLeftCorner()) + if (!yValues.every(isExtendedNumber) || !xValues.every(isExtendedNumber)) { + return new CellError(ErrorType.VALUE, ErrorMessage.NumberRange) + } + const observations = yValues.map(getRawValue) + const numericX = xValues.map(getRawValue) + const predictors = buildPredictorRows(numericX, observations.length, shape) + const fit = fitLinearRegression(predictors, observations, fitIntercept, statistics) + return regressionOutput(fit, statistics, fitIntercept) + }) + } + + /** Predicts output width from predictor geometry and height from a constant stats option. */ + public linestArraySize(ast: ProcedureAst, state: InterpreterState): ArraySize { + if (ast.args.length < 1 || ast.args.length > 4) { + return ArraySize.error() + } + const statistics = staticBoolean(ast.args[3]) + if (statistics === undefined) { + return ArraySize.error() + } + const arrayState = new InterpreterState(state.formulaAddress, true) + const y = this.arraySizeForAst(ast.args[0], arrayState) + const x = ast.args.length < 2 || ast.args[1].type === AstNodeType.EMPTY ? undefined : this.arraySizeForAst(ast.args[1], arrayState) + const shape = regressionShape(y, x) + if (shape === undefined || shape.predictorCount + 1 > this.config.maxColumns) { + return ArraySize.error() + } + return new ArraySize(shape.predictorCount + 1, statistics ? 5 : 1) + } +} diff --git a/src/interpreter/plugin/index.ts b/src/interpreter/plugin/index.ts index 2a3005375b..505776ea2b 100644 --- a/src/interpreter/plugin/index.ts +++ b/src/interpreter/plugin/index.ts @@ -51,3 +51,4 @@ export {StatisticalPlugin} from './StatisticalPlugin' export {MathPlugin} from './MathPlugin' export {ComplexPlugin} from './ComplexPlugin' export {StatisticalAggregationPlugin} from './StatisticalAggregationPlugin' +export {RegressionPlugin} from './RegressionPlugin' diff --git a/src/interpreter/plugin/regression/LinearRegression.ts b/src/interpreter/plugin/regression/LinearRegression.ts new file mode 100644 index 0000000000..86b7b1400a --- /dev/null +++ b/src/interpreter/plugin/regression/LinearRegression.ts @@ -0,0 +1,386 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +/** Numerical fit and uncertainties, before spreadsheet output formatting. */ +export interface LinearRegressionResult { + coefficients: number[], + intercept: number, + standardErrors: number[], + interceptError: number, + residualSumSquares: number, + totalSumSquares: number, + degreesOfFreedom: number, + retainedPredictorCount: number, +} + +/** + * Column-major storage: columns[columnIndex].values[rowIndex] is a matrix entry. + * Factorization overwrites the values with the triangular system. predictorIndex + * preserves the original predictor position through swaps; -1 identifies the intercept. + */ +interface RegressionColumn { + values: number[], + predictorIndex: number, + /** Norm after optional centering, before any reflections. */ + originalNorm: number, +} + +/** Working arrays owned by one fit, with centering offsets kept in the original predictor order. */ +interface PreparedRegression { + columns: RegressionColumn[], + predictorMeans: number[], + transformedResponse: number[], +} + +/** Residual statistics shared by coefficient uncertainty calculation and spreadsheet output. */ +interface ResidualStatistics { + totalSumSquares: number, + residualSumSquares: number, + degreesOfFreedom: number, + residualVariance: number, +} + +/** A reflection and the leading value it produces in the selected column tail. */ +interface HouseholderReflection { + vector: number[], + factor: number, + diagonalValue: number, +} + +/** Coefficient uncertainties restored to the original predictor order. */ +interface CoefficientStandardErrors { + standardErrors: number[], + interceptError: number, +} + +/** + * Rank threshold relative to a predictor's original centered norm. + * Sampled Excel Online cases with proportional predictors retain a 1e-5 perturbation and remove 1e-6. + * Keep the threshold separate from roundoff handling for the fitted statistics. + */ +const RANK_TOLERANCE = 1e-6 + +/** Computes a Euclidean norm without squaring large inputs or spreading an unbounded range. */ +function norm(values: number[], startRow = 0): number { + let result = 0 + for (let rowIndex = startRow; rowIndex < values.length; rowIndex++) { + result = Math.hypot(result, values[rowIndex]) + } + return result +} + +/** Computes the mean relative to the first value to preserve small changes on large offsets. */ +function mean(values: number[]): number { + const origin = values[0] + let sum = 0 + for (const value of values) { + sum += (value - origin) / values.length + } + return origin + sum +} + +/** + * Constructs a reflection for a nonzero column tail without modifying the column. + * + * Let x = values.slice(startRow), r = magnitude, and e0 = [1, 0, ...]. + * Choose s = 1 when x[0] >= 0, otherwise -1. + * The vector v = x / r + s * e0 reflects x to [-s * r, 0, ...]. Adding the same + * sign avoids cancellation when constructing v's first entry. + * Its squared length is v^T * v = 2 * (1 + abs(x[0]) / r), so the reflection + * factor 2 / (v^T * v) simplifies to 1 / (1 + abs(x[0]) / r). + * + * @param values - selected column, including any already processed rows + * @param startRow - first row of the column tail to reflect + * @param magnitude - nonzero norm of values starting at startRow, already computed by the caller + */ +function createHouseholderReflection(values: number[], startRow: number, magnitude: number): HouseholderReflection { + const sign = values[startRow] >= 0 ? 1 : -1 + const vector = values.slice(startRow).map(value => value / magnitude) + vector[0] += sign + const factor = 1 / (1 + Math.abs(values[startRow]) / magnitude) + return {vector, factor, diagonalValue: -sign * magnitude} +} + +/** + * Applies a Householder reflection to values starting at startRow, in place. + * + * For the active tail w and reflection vector v, with ^T denoting transpose: + * dotProduct = v^T * w + * reflectionScale = factor * dotProduct, where factor = 2 / (v^T * v) + * w_new = w - reflectionScale * v + * The subtracted vector is twice the projection onto v: this reverses that + * component while preserving the perpendicular component and the total length. + * + * startRow stays fixed during this call; rowIndex - startRow indexes the shorter v. + * Earlier rows remain unchanged. Compute the full dot product before overwriting w. + */ +function applyHouseholderReflection(values: number[], vector: number[], startRow: number, factor: number): void { + let dotProduct = 0 + for (let rowIndex = startRow; rowIndex < values.length; rowIndex++) { + dotProduct += vector[rowIndex - startRow] * values[rowIndex] + } + const reflectionScale = dotProduct * factor + for (let rowIndex = startRow; rowIndex < values.length; rowIndex++) { + values[rowIndex] -= reflectionScale * vector[rowIndex - startRow] + } +} + +/** + * Selects the largest remaining column tail without changing column order. + * diagonalIndex is both the first candidate column and the first active row. + * Strict comparison keeps the first candidate on ties, preserving Excel's duplicate selection. + */ +function selectPivotColumn(columns: RegressionColumn[], diagonalIndex: number, activeColumnCount: number): number { + let selectedColumnIndex = diagonalIndex + let largestNorm = -1 + for (let columnIndex = diagonalIndex; columnIndex < activeColumnCount; columnIndex++) { + const trailingNorm = norm(columns[columnIndex].values, diagonalIndex) + if (trailingNorm > largestNorm) { + largestNorm = trailingNorm + selectedColumnIndex = columnIndex + } + } + return selectedColumnIndex +} + +/** + * Solves R * solution = rightHandSide from the bottom row upward. + * R[row, column] is stored as columns[column].values[row]. Subtract the already + * solved terms to the right of the diagonal, then divide by the diagonal value. + * Exactly zero diagonals leave a zero coefficient for Excel-compatible zero columns. + */ +function backSubstitute(columns: RegressionColumn[], processedColumnCount: number, rightHandSide: number[]): number[] { + const solution = Array(processedColumnCount).fill(0) + for (let rowIndex = processedColumnCount - 1; rowIndex >= 0; rowIndex--) { + if (columns[rowIndex].values[rowIndex] === 0) { + continue + } + let remainingValue = rightHandSide[rowIndex] + for (let columnIndex = rowIndex + 1; columnIndex < processedColumnCount; columnIndex++) { + remainingValue -= columns[columnIndex].values[rowIndex] * solution[columnIndex] + } + solution[rowIndex] = remainingValue / columns[rowIndex].values[rowIndex] + } + return solution +} + +/** + * Computes standard errors from the retained triangular system, in original predictor order. + * + * For nonsingular R, covariance = residualVariance * inverseR * transpose(inverseR). + * A coefficient's standard error is therefore sqrt(residualVariance) times the norm + * of its row of inverseR. Solve R * inverseColumn = e_j for each unit vector e_j + * to accumulate those row norms without building a full inverse or covariance matrix. + * + * With an intercept, centering changes it by -sum(coefficient * predictorMean), so each + * inverse column contributes interceptWeight = inverseColumn[0] - sum(mean * inverseEntry). + * Zero diagonals retain backSubstitute's Excel-compatible zero-coefficient behavior. + */ +function computeCoefficientStandardErrors( + columns: RegressionColumn[], + processedColumnCount: number, + predictorMeans: number[], + fitIntercept: boolean, + residualVariance: number, +): CoefficientStandardErrors { + const standardErrors = Array(predictorMeans.length).fill(0) + let interceptError = 0 + if (residualVariance === 0) { + return {standardErrors, interceptError} + } + + const interceptColumnCount = fitIntercept ? 1 : 0 + const residualStandardError = Math.sqrt(residualVariance) + for (let inverseColumnIndex = 0; inverseColumnIndex < processedColumnCount; inverseColumnIndex++) { + const rightHandSide = Array(processedColumnCount).fill(0) + rightHandSide[inverseColumnIndex] = 1 + const inverseColumn = backSubstitute(columns, processedColumnCount, rightHandSide) + let interceptWeight = fitIntercept ? inverseColumn[0] : 0 + for (let columnIndex = interceptColumnCount; columnIndex < processedColumnCount; columnIndex++) { + const predictorIndex = columns[columnIndex].predictorIndex + // Scale before hypot: squaring inverse entries first can overflow or underflow. + standardErrors[predictorIndex] = Math.hypot(standardErrors[predictorIndex], residualStandardError * inverseColumn[columnIndex]) + interceptWeight -= predictorMeans[predictorIndex] * inverseColumn[columnIndex] + } + interceptError = Math.hypot(interceptError, residualStandardError * interceptWeight) + } + return {standardErrors, interceptError} +} + +/** + * Copies observations and predictors into working arrays without modifying the inputs. + * With an intercept, uses centeredX = x - mean(x) and shiftedY = y - y[0] to preserve + * small variations on large offsets. Without an intercept, leaves both coordinates unchanged. + * Predictor means retain their input order for recovering b = mean(y) - sum(m * mean(x)). + * When no intercept is fitted, predictorMeans contains zero offsets instead. + * + * @param predictors - observation rows with predictor columns in their original order + * @param observations - numeric response for each observation + * @param fitIntercept - whether to add a column of ones and shift the coordinates + * @returns Working columns in Excel tie order, predictor centering offsets, and the shifted response. + */ +function prepareRegression(predictors: number[][], observations: number[], fitIntercept: boolean): PreparedRegression { + const observationCount = observations.length + const predictorCount = predictors[0].length + const predictorMeans = Array(predictorCount).fill(0) + const columns: RegressionColumn[] = [] + if (fitIntercept) { + columns.push({values: Array(observationCount).fill(1), predictorIndex: -1, originalNorm: Math.sqrt(observationCount)}) + } + for (let predictorIndex = 0; predictorIndex < predictorCount; predictorIndex++) { + const values = predictors.map(row => row[predictorIndex]) + predictorMeans[predictorIndex] = fitIntercept ? mean(values) : 0 + const centeredValues = values.map(value => value - predictorMeans[predictorIndex]) + columns.push({values: centeredValues, predictorIndex, originalNorm: norm(centeredValues)}) + } + + // Excel tie ordering: [intercept, x1, x2, ...] -> [intercept, x2, ..., x1]. + // This matches swapping a trailing intercept into the first input position. + if (fitIntercept && predictorCount > 1) { + const firstPredictor = columns.splice(1, 1)[0] + columns.push(firstPredictor) + } + // Remove a common response offset before reflection to preserve small variations. + const responseOrigin = fitIntercept ? observations[0] : 0 + const transformedResponse = observations.map(value => value - responseOrigin) + return {columns, predictorMeans, transformedResponse} +} + +/** + * Overwrites working columns with a pivoted triangular system and applies the same + * Householder reflections to the response. Retained columns store R; the transformed + * response supplies the right-hand side of R * solution = response and its residual tail. + * Original predictor positions remain available through each column's predictorIndex. + * + * Exactly zero predictors intentionally consume an available slot. Excel Online returns + * a zero coefficient but still excludes that response entry from residual statistics. + * A dependent nonzero predictor is removed without consuming a slot. + * + * @param columns - mutable working columns, including a fixed leading intercept when present + * @param transformedResponse - mutable shifted response, transformed in place alongside the columns + * @param interceptColumnCount - one when fitting an intercept, otherwise zero + * @returns Number of processed columns, including zero columns; this is not mathematical rank. + */ +function factorizeRegressionInPlace(columns: RegressionColumn[], transformedResponse: number[], interceptColumnCount: number): number { + const observationCount = transformedResponse.length + // Columns are partitioned as [processed | candidates | rejected]: + // [0, processedColumnCount), [processedColumnCount, activeColumnCount), and the rest. + // processedColumnCount is also the current diagonal index, not the mathematical rank: + // exactly zero columns still consume a slot for Excel-compatible residual statistics. + let activeColumnCount = columns.length + let processedColumnCount = 0 + while (processedColumnCount < Math.min(observationCount, activeColumnCount)) { + if (processedColumnCount >= interceptColumnCount) { + const pivotColumnIndex = selectPivotColumn(columns, processedColumnCount, activeColumnCount) + const previousColumn = columns[processedColumnCount] + columns[processedColumnCount] = columns[pivotColumnIndex] + columns[pivotColumnIndex] = previousColumn + } + const column = columns[processedColumnCount] + const trailingNorm = norm(column.values, processedColumnCount) + // Rejected predictors do not advance the diagonal. An originally zero column + // gives 0 < 0 here, so it is retained and advances processedColumnCount below. + if (processedColumnCount >= interceptColumnCount && trailingNorm < RANK_TOLERANCE * column.originalNorm) { + activeColumnCount-- + columns[processedColumnCount] = columns[activeColumnCount] + columns[activeColumnCount] = column + continue + } + if (trailingNorm !== 0) { + const {vector, factor, diagonalValue} = createHouseholderReflection(column.values, processedColumnCount, trailingNorm) + for (let columnIndex = processedColumnCount + 1; columnIndex < activeColumnCount; columnIndex++) { + applyHouseholderReflection(columns[columnIndex].values, vector, processedColumnCount, factor) + } + applyHouseholderReflection(transformedResponse, vector, processedColumnCount, factor) + column.values[processedColumnCount] = diagonalValue + column.values.fill(0, processedColumnCount + 1) + } + processedColumnCount++ + } + return processedColumnCount +} + +/** + * Computes response variation and residual error without modifying either response array. + * Total variation is sum((y - mean(y))^2) with an intercept, otherwise sum(y^2). + * Residual error is the squared norm below the processed rows of the transformed response. + * Uses the factorization's Excel-compatible row boundary, including consumed zero columns. + * + * @param observations - response values in their original coordinates + * @param transformedResponse - shifted response after the factorization's reflections + * @param processedColumnCount - number of response entries assigned to the triangular solve + * @param predictorCount - original predictor count, used for the simple-fit roundoff rule + * @param fitIntercept - whether total variation is measured around the response mean + * @returns Sums of squares, remaining degrees of freedom, and residual variance (zero when no degrees remain). + */ +function computeResidualStatistics( + observations: number[], + transformedResponse: number[], + processedColumnCount: number, + predictorCount: number, + fitIntercept: boolean, +): ResidualStatistics { + const observationCount = observations.length + const responseMean = fitIntercept ? mean(observations) : 0 + const centeredResponse = observations.map(value => value - responseMean) + const totalSumSquares = norm(centeredResponse) ** 2 + // Reflections preserve squared error; entries below the solved rows are the residual tail. + let residualSumSquares = norm(transformedResponse, processedColumnCount) ** 2 + // Simple exact fits have zero residual statistics in Excel. Discard only + // reflection roundoff, whose squared scale is O((n * machine epsilon)^2). + const simpleFitRoundoff = predictorCount === 1 && residualSumSquares <= totalSumSquares * (observationCount * Number.EPSILON) ** 2 + if (totalSumSquares === 0 || processedColumnCount === observationCount || simpleFitRoundoff) { + residualSumSquares = 0 + } + const degreesOfFreedom = observationCount - processedColumnCount + const residualVariance = degreesOfFreedom === 0 ? 0 : residualSumSquares / degreesOfFreedom + return {totalSumSquares, residualSumSquares, degreesOfFreedom, residualVariance} +} + +/** + * Fits observations using pivoted Householder QR, with a fixed leading intercept column. + * Prepares columns, factorizes them, restores coefficients, then computes fit statistics. + * Predictor centering preserves accuracy for large offsets. For n observations and k predictors, + * the factorization uses O(n*k + k*k) storage and costs O(n*k*k) for n >= k; Q is never materialized. + * + * @param predictors - observation rows, with predictor columns in their original order + * @param observations - numeric response for each observation + * @param fitIntercept - whether to include a constant term + * @param statistics - whether to compute coefficient uncertainties + */ +export function fitLinearRegression(predictors: number[][], observations: number[], fitIntercept: boolean, statistics: boolean): LinearRegressionResult { + const predictorCount = predictors[0].length + const interceptColumnCount = fitIntercept ? 1 : 0 + const {columns, predictorMeans, transformedResponse} = prepareRegression(predictors, observations, fitIntercept) + const processedColumnCount = factorizeRegressionInPlace(columns, transformedResponse, interceptColumnCount) + + // Solve in pivot order, then restore the caller's predictor order. Rejected slots stay zero. + const solution = backSubstitute(columns, processedColumnCount, transformedResponse) + const coefficients = Array(predictorCount).fill(0) + for (let columnIndex = interceptColumnCount; columnIndex < processedColumnCount; columnIndex++) { + coefficients[columns[columnIndex].predictorIndex] = solution[columnIndex] + } + // Recover the intercept in original coordinates: b = mean(y) - sum(coefficient * mean(x)). + const intercept = fitIntercept + ? mean(observations) - coefficients.reduce((sum, coefficient, predictorIndex) => sum + coefficient * predictorMeans[predictorIndex], 0) + : 0 + + const {totalSumSquares, residualSumSquares, degreesOfFreedom, residualVariance} = computeResidualStatistics( + observations, transformedResponse, processedColumnCount, predictorCount, fitIntercept, + ) + const {standardErrors, interceptError} = statistics + ? computeCoefficientStandardErrors(columns, processedColumnCount, predictorMeans, fitIntercept, residualVariance) + : {standardErrors: Array(predictorCount).fill(0), interceptError: 0} + return { + coefficients, + intercept, + standardErrors, + interceptError, + residualSumSquares, + totalSumSquares, + degreesOfFreedom, + retainedPredictorCount: processedColumnCount - interceptColumnCount, + } +}