Skip to content

_PolyFit

Description

_PolyFit performs a (weighted) polynomial least-squares fit and returns the pieces the Poly* functions need, packed into one array: the coefficients (ascending order), the matrix (X′WX)−1(X^\prime WX)^{-1}, and the residual standard deviation and degrees of freedom. The weights are wi=1/si2w_i = 1/s_i^2; a single number for S gives an ordinary (unweighted) fit.

Used by: PolyConf, PolyHat, PolyInv, PolyPred, PolySlope

Calls: Ones

Source Code

_PolyFit = LAMBDA(
    Y, // dependent values as row or column vector
    X, // independent values as row or column vector
    p, // polynomial order
    S, // std devs of Y (same size as Y), or a single number for equal weights
    LET(
        // filter input data for missing/non-numeric values
        NumFilter, ISNUMBER(X) * ISNUMBER(Y) * ISNUMBER(S),
        Xn, TOCOL(FILTER(X, NumFilter)),
        Yn, TOCOL(FILTER(Y, NumFilter)),
        // IF broadcasts a single S value to every row
        sdn, TOCOL(FILTER(IF(NumFilter, S, 1), NumFilter)),
        n, ROWS(Yn),
        wt, 1 / sdn ^ 2,
        rw, SQRT(wt),
        // X matrix including the column of ones
        Xm, HSTACK(Ones(n), Xn ^ SEQUENCE(, p)),
        // LINEST on the weighted system (every row scaled by SQRT(w))
        m, DROP(LINEST(Yn * rw, Xm * rw, FALSE, TRUE), , -1),
        b, CHOOSECOLS(TAKE(m, 1), SEQUENCE(, p + 1, p + 1, -1)),
        sres, INDEX(m, 3, 2),
        DOF, INDEX(m, 4, 2),
        inner, MINVERSE(MMULT(TRANSPOSE(Xm * wt), Xm)),
        IF(
            OR(sdn <= 0),
            "ERROR: std devs must be positive",
            VSTACK(b, inner, EXPAND(HSTACK(sres, DOF), 1, p + 1, ""))
        )
    )
);