Linear Regression

regression.txt

// slope of model y~x
// inputs:
//    xs: the input data
//    ys: the response data
// return:
//    the slope
slope = lambda(xs, ys,
    let(
        x_bar, average(xs),
        y_bar, average(ys),
        numerator, sum((xs - x_bar) * (ys - y_bar)),
        denominator, sum((xs - x_bar)^2),
        numerator/denominator
    )
);

// intercept of model y ~ x
// inputs:
//     xs: the input data
//     ys: the response data
// return:
//     the intercept    
intercept = lambda(xs, ys,
    let(
        y_bar, average(ys),
        x_bar, average(xs),
        y_bar - jp.slope(xs, ys)*x_bar
    )
);


// least squares fit for a model y ~ x
// This repicates the excel LINEST function
// however the results are returned in an
// order suitable for calculating the modelled
// reponse values from input.
// inputs:
//      xs: the inputs
//      ys: the resonses
// return:
//      a stack (2 x 1) [[b0], [b1]]
least_squares_est =  lambda(xs, ys,
    let(
        x_bar, average(xs),
        y_bar, average(ys),
        numerator, sum((xs - x_bar) * (ys - y_bar)),
        denominator, sum((xs - x_bar)^2),
        b1_, numerator/denominator,
        b0_, y_bar - b1_ * x_bar,
        vstack(b0_, b1_)
    )
);

// calculate predicted response from set of samples, using a linear
// model developed from a series of calibration standards
// inputs:
//      xs: known inputs from calibration curve
//      ys: known responses from calibration curve
//      samples: sample inputs
// return:
//      sequence of predicted responses from the samples
// assume:
//      count(xs) == count(ys)
//      xs, ys, samples are numeric
y_bar =  lambda(xs, ys, samples,
    let(
        n_, count(samples),
        seq_, sequence(n_, 1, 1, 0),
        design_mat, hstack(seq_, samples),
        bs_, jp.least_squares_est(xs, ys),
        mmult(design_mat, bs_)
    )
);