root/trunk/regress.d

/**A module for performing linear regression.  This module has an unusual
 * interface, as it is range-based instead of matrix based. Values for
 * independent variables are provided as either a tuple or a range of ranges.
 * This means that one can use, for example, map, to fit high order models and
 * lazily evaluate certain values.  (For details, see examples below.)
 *
 * Author:  David Simcha*/
 /*
 * License:
 * Boost Software License - Version 1.0 - August 17th, 2003
 *
 * Permission is hereby granted, free of charge, to any person or organization
 * obtaining a copy of the software and accompanying documentation covered by
 * this license (the "Software") to use, reproduce, display, distribute,
 * execute, and transmit the Software, and to prepare derivative works of the
 * Software, and to permit third-parties to whom the Software is furnished to
 * do so, all subject to the following:
 *
 * The copyright notices in the Software and this entire statement, including
 * the above license grant, this restriction and the following disclaimer,
 * must be included in all copies of the Software, in whole or in part, and
 * all derivative works of the Software, unless such copies or derivative
 * works are solely in the form of machine-executable object code generated by
 * a source language processor.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 */
module dstats.regress;

import std.math, std.algorithm, std.traits, std.array, std.traits, std.exception,
    std.typetuple, std.typecons, std.numeric;

import dstats.alloc, std.range, std.conv, dstats.distrib, dstats.cor,
    dstats.base, dstats.summary, dstats.sort;

version(unittest) {
    import std.stdio, dstats.random, std.functional;
    void main(){}
}

///
struct PowMap(ExpType, T)
if(isForwardRange!(T)) {
    Unqual!T range;
    Unqual!ExpType exponent;
    double cache;

    this(T range, ExpType exponent) {
        this.range = range;
        this.exponent = exponent;

        static if(isIntegral!ExpType) {
            cache = pow(cast(double) range.front, exponent);
        } else {
            cache = pow(cast(ExpType) range.front, exponent);
        }
    }

    @property double front() const pure nothrow {
        return cache;
    }

    void popFront() {
        range.popFront;
        if(!range.empty) {
            cache = pow(cast(double) range.front, exponent);
        }
    }

    @property typeof(this) save() {
        return this;
    }

    @property bool empty() {
        return range.empty;
    }
}

/**Maps a forward range to a power determined at runtime.  ExpType is the type
 * of the exponent.  Using an int is faster than using a double, but obviously
 * less flexible.*/
PowMap!(ExpType, T) powMap(ExpType, T)(T range, ExpType exponent) {
    alias PowMap!(ExpType, T) RT;
    return RT(range, exponent);
}

// Very ad-hoc, does a bunch of matrix ops for linearRegress and
// linearRegressBeta.  Specifically, computes xTx and xTy.
// Written specifically to be efficient in the context used here.
private void rangeMatrixMulTrans(U, T...)
(ref double[] xTy, out double[][] xTx, U y, ref T matIn) {
    static if(isArray!(T[0]) &&
        isInputRange!(typeof(matIn[0][0])) && matIn.length == 1) {
        alias matIn[0] mat;
    } else {
        alias matIn mat;
    }

    bool someEmpty() {
        if(y.empty) {
            return true;
        }
        foreach(range; mat) {
            if(range.empty) {
                return true;
            }
        }
        return false;
    }

    void popAll() {
        foreach(ti, range; mat) {
            mat[ti].popFront;
        }
        y.popFront;
    }

    xTy[] = 0;

    xTx = newStack!(double[])(mat.length);
    foreach(ref elem; xTx) {
        elem = newStack!double(mat.length);
    }

    foreach(row; xTx) {
        row[] = 0;
    }

    mixin(newFrame);
    auto deltas = newStack!double(mat.length);

    // Using an algorithm similar to the one for Pearson cor to improve
    // numerical stability.  Calculate means and covariances, then
    // combine them:  Sum of squares = mean1 * N * mean2 + N * cov.
    double k = 0;
    double[] means = newStack!(double)(mat.length);
    means[] = 0;
    double yMean = 0;

    while(!someEmpty) {
        foreach(i, elem; mat) {
            deltas[i] = cast(double) elem.front;
        }

        immutable kMinus1 = k;
        immutable kNeg1 = 1 / ++k;
        deltas[] -= means[];
        means[] += deltas[] * kNeg1;

        immutable yDelta = cast(double) y.front - yMean;
        yMean += yDelta * kNeg1;

        foreach(i, delta1; deltas) {
            xTy[i] += kMinus1 * delta1 * kNeg1 * yDelta;
            xTx[i][i] += kMinus1 * delta1 * kNeg1 * delta1;
            auto xTxi = xTx[i];

            foreach(j, delta2; deltas[0..i]) {
                xTxi[j] += kMinus1 * delta1 * kNeg1 * delta2;
            }
        }
        popAll();
    }

    // k is n now that we're done looping over the data.
    alias k n;

    // mat now consists of covariance * n.    Divide by n and add mean1 * mean2
    // to get sum of products.
    foreach(i; 0..xTx.length) foreach(j; 0..i + 1) {
        xTx[i][j] /= n;
        xTx[i][j] += means[i] * means[j];
    }

    // Similarly for the xTy vector
    foreach(i, ref elem; xTy) {
        xTy[i] /= n;
        elem += yMean * means[i];
    }
    symmetrize(xTx);
}

// Copies values from lower triangle to upper triangle.
private void symmetrize(double[][] mat) pure nothrow {
    foreach(i; 1..mat.length) {
        foreach(j; 0..i) {
            mat[j][i] = mat[i][j];
        }
    }
}

// Uses Gauss-Jordan elim. w/ row pivoting to invert from.  Stores the results
// in to and leaves from in an undefined state.
package void invert(double[][] from, double[][] to) {
    // Normalize.
    foreach(i, row; from) {
        double absMax = 1.0 / reduce!(max)(map!(abs)(row[0..from.length]));
        row[] *= absMax;
        to[i][] = 0;
        to[i][i] = absMax;
    }

    foreach(col; 0..from.length) {
        size_t bestRow;
        double biggest = 0;
        foreach(row; col..from.length) {
            if(abs(from[row][col]) > biggest) {
                bestRow = row;
                biggest = abs(from[row][col]);
            }
        }

        swap(from[col], from[bestRow]);
        swap(to[col], to[bestRow]);
        immutable pivotFactor = from[col][col];

        foreach(row; 0..from.length) if(row != col) {
            immutable ratio = from[row][col] / pivotFactor;

            // If you're ever looking to optimize this code, good luck.  The
            // bottleneck is almost ENTIRELY this one line:
            from[row][] -= from[col][] * ratio;
            to[row][] -= to[col][] * ratio;
        }
    }

    foreach(i; 0..from.length) {
        immutable diagVal = from[i][i];
        from[i][] /= diagVal;
        to[i][] /= diagVal;
    }
}

unittest {
    auto mat = [[1.0, 2, 3], [4.0, 7, 6], [7.0, 8, 9]];
    auto toMat = [new double[3], new double[3], new double[3]];
    invert(mat, toMat);
    assert(approxEqual(toMat[0], [-0.625, -0.25, 0.375]));
    assert(approxEqual(toMat[1], [-0.25, 0.5, -0.25]));
    assert(approxEqual(toMat[2], [0.708333, -0.25, 0.041667]));
}

// Solve a system of linear equations mat * b = vec for b.  The result is
// stored in vec, and mat is left in an undefined state. Uses Gaussian
// elimination w/ row pivoting.  Roughly 4x faster than using inversion to
// solve the system, and uses roughly half the memory.
private void solve(double[][] mat, double[] vec) {
    // Normalize.
    foreach(i, row; mat) {
        double absMax = 1.0 / reduce!(max)(map!(abs)(row[0..mat.length]));
        row[] *= absMax;
        vec[i] *= absMax;
    }

    foreach(col; 0..mat.length - 1) {
        size_t bestRow;
        double biggest = 0;
        foreach(row; col..mat.length) {
            if(abs(mat[row][col]) > biggest) {
                bestRow = row;
                biggest = abs(mat[row][col]);
            }
        }

        swap(mat[col], mat[bestRow]);
        swap(vec[col], vec[bestRow]);
        immutable pivotFactor = mat[col][col];

        foreach(row; col + 1..mat.length) {
            immutable ratio = mat[row][col] / pivotFactor;

            // If you're ever looking to optimize this code, good luck.  The
            // bottleneck is almost ENTIRELY this one line:
            mat[row][col..$] -= mat[col][col..$] * ratio;
            vec[row] -= vec[col] * ratio;
        }
    }

    foreach(i; 0..mat.length) {
        double diagVal = mat[i][i];
        mat[i][] /= diagVal;
        vec[i] /= diagVal;
    }

    // Do back substitution.
    for(size_t row = mat.length - 1; row != size_t.max; row--) {
        auto v1 = vec[row + 1..$];
        auto v2 = mat[row][$ - v1.length..$];
        vec[row] -= dotProduct(v1, v2);
    }
}

unittest {
    auto mat = [[2.0, 1, -1], [-3.0, -1, 2], [-2.0, 1, 2]];
    auto vec = [8.0, -11, -3];
    solve(mat, vec);
    assert(approxEqual(vec, [2, 3, -1]));

    auto mat2 = [[1.0, 2, 3], [4.0, 7, 6], [7.0, 8, 9]];
    auto vec2 = [8.0, 6, 7];
    solve(mat2, vec2);
    assert(approxEqual(vec2, [-3.875, -0.75, 4.45833]));
}

// Cholesky decomposition functions adapted from Don Clugston's MathExtra
// lib, used w/ permission.
private void choleskyDecompose(double[][] a, double[] diag) {
    immutable N = diag.length;

    foreach(i; 0..N) {
        const ai = a[i];
        double sum = ai[i];

        for(sizediff_t k = i - 1; k >= 0; --k) {
            sum -= ai[k] * ai[k];
        }

        if (sum > 0.0) {
            diag[i] = sqrt(sum);

            foreach(j; i + 1..N) {
                sum = ai[j] - dotProduct(ai[0..i], a[j][0..i]);
                a[j][i] = sum / diag[i];
            }
        } else {
            // not positive definite (could be caused by rounding errors)
            diag[i] = 0;
            // make this whole row zero so they have no further effect
            foreach(j; i + 1..N) a[j][i] = 0;
        }
    }
}

private void choleskySolve(double[][] a, double[] diag, double[] b, double[] x) {
    immutable N = x.length;

    foreach(i; 0..N) {
        const ai = a[i];

        if(diag[i] > 0)  {
            double sum = b[i];
            sum -= dotProduct(ai[0..i], x[0..i]);
            x[i] = sum / diag[i];
        } else x[i] = 0; // skip pos definite rows
    }

    for(sizediff_t i = N - 1; i >= 0; --i) {
        if (diag[i] > 0) {
            double sum = x[i];
            for(sizediff_t k = i + 1; k < N; ++k) sum -= a[k][i] * x[k];
            x[i] = sum / diag[i];
        } else x[i] = 0; // skip pos definite rows
    }
    // Convert failed columns in solution to NANs if required.
    foreach(i; 0..N) {
        if(diag[i] !> 0) x[i] = double.nan;
    }
}

private void choleskySolve(double[][] a, double[] b, double[] x) {
    mixin(newFrame);
    auto diag = newStack!double(x.length);
    choleskyDecompose(a, diag);
    choleskySolve(a, diag, b, x);
}

/**Struct that holds the results of a linear regression.  It's a plain old
 * data struct.*/
struct RegressRes {
    /**The coefficients, one for each range in X.  These will be in the order
     * that the X ranges were passed in.*/
    double[] betas;

    /**The standard error terms of the X ranges passed in.*/
    double[] stdErr;

    /**The lower confidence bounds of the beta terms, at the confidence level
     * specificied.  (Default 0.95).*/
    double[] lowerBound;

    /**The upper confidence bounds of the beta terms, at the confidence level
     * specificied.  (Default 0.95).*/
    double[] upperBound;

    /**The P-value for the alternative that the corresponding beta value is
     * different from zero against the null that it is equal to zero.*/
    double[] p;

    /**The coefficient of determination.*/
    double R2;

    /**The adjusted coefficient of determination.*/
    double adjustedR2;

    /**The root mean square of the residuals.*/
    double residualError;

    /**The P-value for the model as a whole.  Based on an F-statistic.  The
     * null here is that the model has no predictive value, the alternative
     * is that it does.*/
    double overallP;

    // Just used internally.
    private static string arrStr(T)(T arr) {
        return text(arr)[1..$ - 1];
    }

    /**Print out the results in the default format.*/
    string toString() {
        return "Betas:  " ~ arrStr(betas) ~ "\nLower Conf. Int.:  " ~
            arrStr(lowerBound) ~ "\nUpper Conf. Int.:  " ~ arrStr(upperBound) ~
            "\nStd. Err:  " ~ arrStr(stdErr) ~ "\nP Values:  " ~ arrStr(p) ~
            "\nR^2:  " ~ text(R2) ~
            "\nAdjusted R^2:  " ~ text(adjustedR2) ~
            "\nStd. Residual Error:  " ~ text(residualError)
            ~ "\nOverall P:  " ~ text(overallP);
    }
}

/**Forward Range for holding the residuals from a regression analysis.*/
struct Residuals(F, U, T...) {
    static if(T.length == 1 && isForwardRange!(typeof(T[0].front()))) {
        alias T[0] R;
        alias typeof(array(R.init)) XType;
        enum bool needDup = true;
    } else {
        alias T R;
        alias staticMap!(Unqual, R) XType;
        enum bool needDup = false;
    }

    Unqual!U Y;
    XType X;
    F[] betas;
    double residual;
    bool _empty;

    void nextResidual() {
        double sum = 0;
        size_t i = 0;
        foreach(elem; X) {
            double frnt = elem.front;
            sum += frnt * betas[i];
            i++;
        }
        residual = Y.front - sum;
    }

    this(F[] betas, U Y, R X) {
        static if(is(typeof(X.length))) {
            dstatsEnforce(X.length == betas.length,
                "Betas and X must have same length for residuals.");
        } else {
            dstatsEnforce(walkLength(X) == betas.length,
                "Betas and X must have same length for residuals.");
        }

        static if(needDup) {
            this.X = array(X);
        } else {
            this.X = X;
        }

        foreach(i, elem; this.X) {
            static if(isForwardRange!(typeof(elem))) {
                this.X[i] = this.X[i].save;
            }
        }

        this.Y = Y;
        this.betas = betas;
        if(Y.empty) {
            _empty = true;
            return;
        }
        foreach(elem; X) {
            if(elem.empty) {
                _empty = true;
                return;
            }
        }
        nextResidual;
    }

    @property double front() const pure nothrow {
        return residual;
    }

    void popFront() {
        Y.popFront;
        if(Y.empty) {
            _empty = true;
            return;
        }
        foreach(ti, elem; X) {
            X[ti].popFront;
            if(X[ti].empty) {
                _empty = true;
                return;
            }
        }
        nextResidual;
    }

    @property bool empty() const pure nothrow {
        return _empty;
    }

    @property typeof(this) save() {
        auto ret = this;
        ret.Y = ret.Y.save;
        foreach(ti, xElem; ret.X) {
            ret.X[ti] = ret.X[ti].save;
        }

        return ret;
    }
}

/**Given the beta coefficients from a linear regression, and X and Y values,
 * returns a range that lazily computes the residuals.
 */
Residuals!(F, U, T) residuals(F, U, T...)(F[] betas, U Y, T X)
if(isFloatingPoint!F && isForwardRange!U && allSatisfy!(isForwardRange, T)) {
    alias Residuals!(F, U, T) RT;
    return RT(betas, Y, X);
}

// Compiles summary statistics while iterating, to allow ridge regression over
// input ranges.
private struct SummaryIter(R) {
    R range;
    MeanSD summ;

    this(R range) {
        this.range = range;
    }

    double front() @property {
        return range.front;
    }

    void popFront() {
        summ.put(range.front);
        range.popFront();
    }

    bool empty() @property {
        return range.empty;
    }

    double mse() @property const pure nothrow { return summ.mse; }

    double N() @property const pure nothrow { return summ.N; }
}

private template SummaryType(R) {
    alias SummaryIter!R SummaryType;
}

/**
Perform a linear regression and return just the beta values.  The advantages
to just returning the beta values are that it's faster and that each range
needs to be iterated over only once, and thus can be just an input range.
The beta values are returned such that the smallest index corresponds to
the leftmost element of X.  X can be either a tuple or a range of input
ranges.  Y must be an input range.

If, after all X variables are passed in, a numeric type is passed as the last
parameter, this is treated as a ridge parameter and ridge regression is
performed.  Ridge regression is a form of regression that penalizes the L2 norm
of the beta vector and therefore results in more parsimonious models.
However, it makes statistical inference such as that supported by
linearRegress() difficult to impossible.  Therefore, linearRegress() doesn't
support ridges.

If no ridge parameter is passed, or equivalently if the ridge parameter is
zero, then ordinary least squares regression is performed.

Notes:  The X ranges are traversed in lockstep, but the traversal is stopped
at the end of the shortest one.  Therefore, using infinite ranges is safe.
For example, using repeat(1) to get an intercept term works.

References:

http://www.mathworks.com/help/toolbox/stats/ridge.html

Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S.
Fourth Edition. Springer, New York. ISBN 0-387-95457-0
(This is the citation for the MASS R package.)

Examples:
---
int[] nBeers = [8,6,7,5,3,0,9];
int[] nCoffees = [3,6,2,4,3,6,8];
int[] musicVolume = [3,1,4,1,5,9,2];
int[] programmingSkill = [2,7,1,8,2,8,1];
double[] betas = linearRegressBeta(programmingSkill, repeat(1), nBeers, nCoffees,
    musicVolume, map!"a * a"(musicVolume));

// Now throw in a ridge parameter of 2.5.
double[] ridgeBetas = linearRegressBeta(programmingSkill, repeat(1), nBeers,
    nCoffees, musicVolume, map!"a * a"(musicVolume), 2.5);
---
 */
double[] linearRegressBeta(U, T...)(U Y, T XIn)
if(doubleInput!(U)) {
    double[] dummy;
    return linearRegressBetaBuf!(U, T)(dummy, Y, XIn);
}

/**
Same as linearRegressBeta, but allows the user to specify a buffer for
the beta terms.  If the buffer is too short, a new one is allocated.
Otherwise, the results are returned in the user-provided buffer.
 */
double[] linearRegressBetaBuf(U, TRidge...)(double[] buf, U Y, TRidge XRidge)
if(doubleInput!(U)) {
    mixin(newFrame);

    static if(isFloatingPoint!(TRidge[$ - 1]) || isIntegral!(TRidge[$ - 1])) {
        // ridge param.
        alias XRidge[$ - 1] ridgeParam;
        alias TRidge[0..$ - 1] T;
        alias XRidge[0..$ - 1] XIn;
        enum bool ridge = true;
        dstatsEnforce(ridgeParam >= 0,
            "Cannot do ridge regerssion with ridge param <= 0.");

        static SummaryIter!R summaryIter(R)(R range) {
            return typeof(return)(range);
        }

    } else {
        enum bool ridge = false;
        enum ridgeParam = 0;
        alias TRidge T;
        alias XRidge XIn;

        static R summaryIter(R)(R range) {
            return range;
        }
    }

    static if(isArray!(T[0]) && isInputRange!(typeof(XIn[0][0])) &&
        T.length == 1) {
        auto X = tempdup(map!(summaryIter)(XIn[0]));
        alias typeof(X[0]) E;
    } else {
        static if(ridge) {
            alias staticMap!(SummaryType, T) XType;
            XType X;

            foreach(ti, elem; XIn) {
                X[ti] = summaryIter(elem);
            }
        } else {
            alias XIn X;
        }
    }

    double[][] xTx;
    double[] xTy = newStack!double(X.length);
    double[] ret;
    if(buf.length < X.length) {
        ret = new double[X.length];
    } else {
        ret = buf[0..X.length];
    }

    rangeMatrixMulTrans(xTy, xTx, Y, X);

    static if(ridge) {
        if(ridgeParam > 0) {
            foreach(i, range; X) {
                // The whole matrix is scaled by N, as well as the xTy vector,
                // so scale the ridge param similarly.
                xTx[i][i] += ridgeParam * range.mse / range.N;
            }
        }
    }

    choleskySolve(xTx, xTy, ret);
    return ret;
}

/**
Perform a linear regression as in linearRegressBeta, but return a
RegressRes with useful stuff for statistical inference.  If the last element
of input is a real, this is used to specify the confidence intervals to
be calculated.  Otherwise, the default of 0.95 is used.  The rest of input
should be the elements of X.

When using this function, which provides several useful statistics useful
for inference, each range must be traversed twice.  This means:

1.  They have to be forward ranges, not input ranges.

2.  If you have a large amount of data and you're mapping it to some
    expensive function, you may want to do this eagerly instead of lazily.

Notes:

The X ranges are traversed in lockstep, but the traversal is stopped
at the end of the shortest one.  Therefore, using infinite ranges is safe.
For example, using repeat(1) to get an intercept term works.

If the confidence interval specified is exactly 0, this is treated as a
special case and confidence interval calculation is skipped.  This can speed
things up significantly and therefore can be useful in monte carlo and possibly
data mining contexts.

Bugs:  The statistical tests performed in this function assume that an
intercept term is included in your regression model.  If no intercept term
is included, the P-values, confidence intervals and adjusted R^2 values
calculated by this function will be wrong.

Examples:
---
int[] nBeers = [8,6,7,5,3,0,9];
int[] nCoffees = [3,6,2,4,3,6,8];
int[] musicVolume = [3,1,4,1,5,9,2];
int[] programmingSkill = [2,7,1,8,2,8,1];

// Using default confidence interval:
auto results = linearRegress(programmingSkill, repeat(1), nBeers, nCoffees,
    musicVolume, map!"a * a"(musicVolume));

// Using user-specified confidence interval:
auto results = linearRegress(programmingSkill, repeat(1), nBeers, nCoffees,
    musicVolume, map!"a * a"(musicVolume), 0.8675309);
---
*/
RegressRes linearRegress(U, TC...)(U Y, TC input) {
    static if(is(TC[$ - 1] : double)) {
        double confLvl = input[$ - 1];
        enforceConfidence(confLvl);
        alias TC[0..$ - 1] T;
        alias input[0..$ - 1] XIn;
    } else {
        double confLvl = 0.95; // Default;
        alias TC T;
        alias input XIn;
    }

    mixin(newFrame);
    static if(isForwardRange!(T[0]) && isForwardRange!(typeof(XIn[0].front())) &&
        T.length == 1) {

        enum bool arrayX = true;
        alias typeof(XIn[0].front) E;
        E[] X = tempdup(XIn[0]);
    } else static if(allSatisfy!(isForwardRange, T)) {
        enum bool arrayX = false;
        alias XIn X;
    } else {
        static assert(0, "Linear regression can only be performed with " ~
            "tuples of forward ranges or ranges of forward ranges.");
    }

    double[][] xTx = newStack!(double[])(X.length),
        xTxNeg1 = newStack!(double[])(X.length);

    foreach(i; 0..X.length) {
        xTx[i] = newStack!double(X.length);
    }

    double[] xTy = newStack!double(X.length);

    foreach(i; 0..X.length) {
        xTxNeg1[i] = newStack!double(X.length);
    }

    static if(arrayX) {
        auto xSaved = X.tempdup;
        foreach(ref elem; xSaved) {
            elem = elem.save;
        }
    } else {
        auto xSaved = X;
        foreach(ti, Type; X) {
            xSaved[ti] = X[ti].save;
        }
    }

    rangeMatrixMulTrans(xTy, xTx, Y.save, X);
    invert(xTx, xTxNeg1);
    double[] betas = new double[X.length];
    foreach(i; 0..betas.length) {
        betas[i] = 0;
        foreach(j; 0..betas.length) {
            betas[i] += xTxNeg1[i][j] * xTy[j];
        }
    }

    X = xSaved;
    auto residuals = residuals(betas, Y, X);
    double S = 0;
    ulong n = 0;
    PearsonCor R2Calc;
    for(; !residuals.empty; residuals.popFront) {
        double residual = residuals.front;
        S += residual * residual;
        double Yfront = residuals.Y.front;
        double predicted = Yfront - residual;
        R2Calc.put(predicted, Yfront);
        n++;
    }
    immutable ulong df =  n - X.length;
    immutable double R2 = R2Calc.cor ^^ 2;
    immutable double adjustedR2 = 1.0L - (1.0L - R2) * ((n - 1.0L) / df);

    immutable double sigma2 = S / (n - X.length);

    double[] stdErr = new double[betas.length];
    foreach(i, ref elem; stdErr) {
        elem = sqrt( S * xTxNeg1[i][i] / df / n);
    }

    double[] lowerBound, upperBound;
    if(confLvl == 0) {
        // Then we're going to skip the computation to save time.  (See below.)
        lowerBound = betas;
        upperBound = betas;
    } else {
        lowerBound = new double[betas.length];
        upperBound = new double[betas.length];
    }
    auto p = new double[betas.length];

    foreach(i, beta; betas) {
        try {
            p[i] = 2 * min(studentsTCDF(beta / stdErr[i], df),
                           studentsTCDFR(beta / stdErr[i], df));
        } catch(DstatsArgumentException) {
            // Leave it as a NaN.
        }

        if(confLvl > 0) {
            // Skip confidence level computation if level is zero, to save
            // computation time.  This is important in monte carlo and possibly
            // data mining contexts.
            try {
                double delta = invStudentsTCDF(0.5 * (1 - confLvl), df) *
                     stdErr[i];
                upperBound[i] = beta - delta;
                lowerBound[i] = beta + delta;
            } catch(DstatsArgumentException) {
                // Leave confidence bounds as NaNs.
            }
        }
    }

    double F = (R2 / (X.length - 1)) / ((1 - R2) / (n - X.length));
    double overallP;
    try {
        overallP = fisherCDFR(F, X.length - 1, n - X.length);
    } catch(DstatsArgumentException) {
        // Leave it as a NaN.
    }

    return RegressRes(betas, stdErr, lowerBound, upperBound, p, R2,
        adjustedR2, sqrt(sigma2), overallP);
}


/**Struct returned by polyFit.*/
struct PolyFitRes(T) {

    /**The array of PowMap ranges created by polyFit.*/
    T X;

    /**The rest of the results.  This is alias this'd.*/
    RegressRes regressRes;
    alias regressRes this;
}

/**Convenience function that takes a forward range X and a forward range Y,
 * creates an array of PowMap structs for integer powers from 0 through N,
 * and calls linearRegressBeta.
 *
 * Returns:  An array of doubles.  The index of each element corresponds to
 * the exponent.  For example, the X<sup>2</sup> term will have an index of
 * 2.
 */
double[] polyFitBeta(T, U)(U Y, T X, uint N, double ridge = 0) {
    double[] dummy;
    return polyFitBetaBuf!(T, U)(dummy, Y, X, N);
}

/**Same as polyFitBeta, but allows the caller to provide an explicit buffer
 * to return the coefficients in.  If it's too short, a new one will be
 * allocated.  Otherwise, results will be returned in the user-provided buffer.
 */
double[] polyFitBetaBuf(T, U)(double[] buf, U Y, T X, uint N, double ridge = 0) {
    mixin(newFrame);
    auto pows = newStack!(PowMap!(uint, T))(N + 1);
    foreach(exponent; 0..N + 1) {
        pows[exponent] = powMap(X, exponent);
    }

    if(ridge == 0) {
        return linearRegressBetaBuf(buf, Y, pows);
    } else {
        return linearRegressBetaBuf(buf, Y, pows, ridge);
    }
}

/**Convenience function that takes a forward range X and a forward range Y,
 * creates an array of PowMap structs for integer powers 0 through N,
 * and calls linearRegress.
 *
 * Returns:  A PolyFitRes containing the array of PowMap structs created and
 * a RegressRes.  The PolyFitRes is alias this'd to the RegressRes.*/
PolyFitRes!(PowMap!(uint, T)[])
polyFit(T, U)(U Y, T X, uint N, double confInt = 0.95) {
    enforceConfidence(confInt);
    auto pows = new PowMap!(uint, T)[N + 1];
    foreach(exponent; 0..N + 1) {
        pows[exponent] = powMap(X, exponent);
    }
    alias PolyFitRes!(typeof(pows)) RT;
    RT ret;
    ret.X = pows;
    ret.regressRes = linearRegress(Y, pows, confInt);
    return ret;
}

unittest {
    // These are a bunch of values gleaned from various examples on the Web.
    double[] heights = [1.47,1.5,1.52,1.55,1.57,1.60,1.63,1.65,1.68,1.7,1.73,1.75,
        1.78,1.8,1.83];
    double[] weights = [52.21,53.12,54.48,55.84,57.2,58.57,59.93,61.29,63.11,64.47,
        66.28,68.1,69.92,72.19,74.46];
    float[] diseaseSev = [1.9, 3.1, 3.3, 4.8, 5.3, 6.1, 6.4, 7.6, 9.8, 12.4];
    ubyte[] temperature = [2,1,5,5,20,20,23,10,30,25];

    // Values from R.
    auto res1 = polyFit(diseaseSev, temperature, 1);
    assert(approxEqual(res1.betas[0], 2.6623));
    assert(approxEqual(res1.betas[1], 0.2417));
    assert(approxEqual(res1.stdErr[0], 1.1008));
    assert(approxEqual(res1.stdErr[1], 0.0635));
    assert(approxEqual(res1.p[0], 0.0419));
    assert(approxEqual(res1.p[1], 0.0052));
    assert(approxEqual(res1.R2, 0.644));
    assert(approxEqual(res1.adjustedR2, 0.6001));
    assert(approxEqual(res1.residualError, 2.03));
    assert(approxEqual(res1.overallP, 0.00518));


    auto res2 = polyFit(weights, heights, 2);
    assert(approxEqual(res2.betas[0], 128.813));
    assert(approxEqual(res2.betas[1], -143.162));
    assert(approxEqual(res2.betas[2], 61.960));

    assert(approxEqual(res2.stdErr[0], 16.308));
    assert(approxEqual(res2.stdErr[1], 19.833));
    assert(approxEqual(res2.stdErr[2], 6.008));

    assert(approxEqual(res2.p[0], 4.28e-6));
    assert(approxEqual(res2.p[1], 1.06e-5));
    assert(approxEqual(res2.p[2], 2.57e-7));

    assert(approxEqual(res2.R2, 0.9989, 0.0001));
    assert(approxEqual(res2.adjustedR2, 0.9987, 0.0001));

    assert(approxEqual(res2.lowerBound[0], 92.9, 0.01));
    assert(approxEqual(res2.lowerBound[1], -186.8, 0.01));
    assert(approxEqual(res2.lowerBound[2], 48.7, 0.01));
    assert(approxEqual(res2.upperBound[0], 164.7, 0.01));
    assert(approxEqual(res2.upperBound[1], -99.5, 0.01));
    assert(approxEqual(res2.upperBound[2], 75.2, 0.01));

    auto res3 = linearRegress(weights, repeat(1), heights, map!"a * a"(heights));
    assert(res2.betas == res3.betas);

    double[2] beta1Buf;
    auto beta1 = linearRegressBetaBuf
        (beta1Buf[], diseaseSev, repeat(1), temperature);
    assert(beta1Buf.ptr == beta1.ptr);
    assert(beta1Buf[] == beta1[]);
    assert(approxEqual(beta1, res1.betas));
    auto beta2 = polyFitBeta(weights, heights, 2);
    assert(approxEqual(beta2, res2.betas));

    auto res4 = linearRegress(weights, repeat(1), heights);
    assert(approxEqual(res4.p, 3.604e-14));
    assert(approxEqual(res4.betas, [-39.062, 61.272]));
    assert(approxEqual(res4.p, [6.05e-9, 3.60e-14]));
    assert(approxEqual(res4.R2, 0.9892));
    assert(approxEqual(res4.adjustedR2, 0.9884));
    assert(approxEqual(res4.residualError, 0.7591));
    assert(approxEqual(res4.lowerBound, [-45.40912, 57.43554]));
    assert(approxEqual(res4.upperBound, [-32.71479, 65.10883]));

    // Test residuals.
    assert(approxEqual(residuals(res4.betas, weights, repeat(1), heights),
        [1.20184170, 0.27367611,  0.40823237, -0.06993322,  0.06462305,
         -0.40354255, -0.88170814,  -0.74715188, -0.76531747, -0.63076120,
         -0.65892680, -0.06437053, -0.08253613,  0.96202014,  1.39385455]));

    // Test nonzero ridge parameters.
        // Values from R's MASS package.
    auto a = [1, 2, 3, 4, 5, 6, 7];
    auto b = [8, 6, 7, 5, 3, 0, 9];
    auto c = [2, 7, 1, 8, 2, 8, 1];

    // With a ridge param. of zero, ridge regression reduces to regular
    // OLS regression.
    assert(approxEqual(linearRegressBeta(a, repeat(1), b, c, 0),
        linearRegressBeta(a, repeat(1), b, c)));

    // Test the ridge regression. Values from R MASS package.
    auto ridge1 = linearRegressBeta(a, repeat(1), b, c, 1);
    auto ridge2 = linearRegressBeta(a, repeat(1), b, c, 2);
    auto ridge3 = linearRegressBeta(c, [[1,1,1,1,1,1,1], a, b], 10);
    assert(approxEqual(ridge1, [6.0357757, -0.2729671, -0.1337131]));
    assert(approxEqual(ridge2, [5.62367784, -0.22449854, -0.09775174]));
    assert(approxEqual(ridge3, [5.82653624, -0.05197246, -0.27185592 ]));
}

private MeanSD[] calculateSummaries(X...)(X xIn) {
    // This is slightly wasteful because it sticks this shallow dup in
    // an unfreeable pos on TempAlloc.
    static if(X.length == 1 && isRoR!(X[0])) {
        auto ret = newStack!MeanSD(xIn[0].length);
        mixin(newFrame);
        auto x = tempdup(xIn[0]);

        foreach(ref range; x) {
            range = range.save;
        }
    } else {
        auto ret = newStack!MeanSD(xIn.length);
        alias xIn x;

        foreach(ti, R; X) {
            x[ti] = x[ti].save;
        }
    }

    ret[] = MeanSD.init;

    bool someEmpty() {
        foreach(range; x) {
            if(range.empty) return true;
        }

        return false;
    }

    void popAll() {
        foreach(ti, elem; x) {
            x[ti].popFront();
        }
    }

    while(!someEmpty) {
        foreach(i, range; x) {
            ret[i].put(range.front);
        }
        popAll();
    }

    return ret;
}

private double softThresh(double z, double gamma) {
    if(gamma >= abs(z)) {
        return 0;
    } else if(z > 0) {
        return z - gamma;
    } else {
        return z + gamma;
    }
}

private struct PreprocessedData {
    MeanSD[] xSumm;
    MeanSD ySumm;
    double[] y;
    double[][] x;
}

private PreprocessedData preprocessStandardize(Y, X...)
(Y yIn, X xIn) {
    static if(X.length == 1 && isRoR!(X[0])) {
        auto xRaw = tempdup(xIn[0]);
    } else {
        alias xIn xRaw;
    }

    auto summaries = calculateSummaries(xRaw);
    immutable minLen = to!size_t(
        reduce!min(
            map!"a.N"(summaries)
        )
    );

    auto x = newStack!(double[])(summaries.length);
    foreach(i, range; xRaw) {
        x[i] = tempdup(map!(to!double)(take(range, minLen)));
        x[i][] -= summaries[i].mean;
        x[i][] /= sqrt(summaries[i].mse);
    }

    double[] y;
    MeanSD ySumm;
    if(yIn.length) {
        y = tempdup(map!(to!double)(take(yIn, minLen)));
        ySumm = meanStdev(y);
        y[] -= ySumm.mean;
    }

    return PreprocessedData(summaries, ySumm, y, x);
}

/**
Performs lasso (L1) and/or ridge (L2) penalized linear regression.  Due to the
way the data is standardized, no intercept term should be included in x
(unlike linearRegress and linearRegressBeta).  The intercept coefficient is
implicitly included and returned in the first element of the returned array.
Usage is otherwise identical.

Note:  Setting lasso equal to zero is equivalent to performing ridge regression.
       This can also be done with linearRegressBeta.  However, the
       linearRegressBeta algorithm is optimized for memory efficiency and
       large samples.  This algorithm is optimized for large feature sets.

Returns:  The beta coefficients for the regression model.

References:

Friedman J, et al Pathwise coordinate optimization. Ann. Appl. Stat.
2007;2:302-332.

Goeman, J. J., L1 penalized estimation in the Cox proportional hazards model.
Biometrical Journal 52(1), 70{84.

Eilers, P., Boer, J., Van Ommen, G., Van Houwelingen, H. 2001 Classification of
microarray data with penalized logistic regression. Proceedings of SPIE.
Progress in Biomedical Optics and Images vol. 4266, pp. 187-198
*/
double[] linearRegressPenalized(Y, X...)
(Y yIn, X xIn, double lasso, double ridge) {
    mixin(newFrame);

    auto preproc = preprocessStandardize(yIn, xIn);

    auto summaries = preproc.xSumm;
    auto ySumm = preproc.ySumm;
    auto x = preproc.x;
    auto y = preproc.y;

    auto betasFull = new double[x.length + 1];
    betasFull[] = 0;
    auto betas = betasFull[1..$];

    if(lasso > 0) {
        coordDescent(y, x, betas, lasso, ridge, null);
    } else if(y.length > x.length) {
        // Correct for different )#*$# scaling conventions.
        foreach(i, feature; x) {
            feature[] /= sqrt(summaries[i].mse);
        }

        linearRegressBetaBuf(betas, y, x, ridge);

        // More correction for diff. scaling factors.
        foreach(i, ref b; betas) {
            b /= sqrt(summaries[i].mse);
        }

    } else {
        ridgeLargeP(y, x, ridge, betas, null);
    }

    foreach(i, ref elem; betas) {
        elem /= sqrt(summaries[i].mse);
    }

    betasFull[0] = ySumm.mean;
    foreach(i, beta; betas) {
        betasFull[0] -= beta * summaries[i].mean;
    }

    return betasFull;
}

private void coordDescent
(double[] y, double[][] x, double[] betas, double lasso, double ridge, double[] w) {
    mixin(newFrame);

    auto predictions = newStack!double(y.length);
    predictions[] = 0;

    void makePredictions() {
        foreach(j, beta; betas) {
            predictions[] += x[j][] * beta;
        }
    }

    if(reduce!max(0.0, map!abs(betas)) > 0) {
        makePredictions();
    }

    auto residuals = newStack!double(y.length);

    uint iter = 0;
    enum maxIter = 10_000;
    enum relEpsilon = 1e-5;
    enum absEpsilon = 1e-10;
    immutable n = cast(double) y.length;
    auto perm = tempdup(iota(0U, x.length));

    auto weightDots = newStack!double(x.length);
    if(w.length == 0) {
        ridge /= n;
        lasso /= n;
        weightDots[] = 1;
    } else {
        foreach(j, col; x) {
            weightDots[j] = dotProduct(w, map!"a * a"(col));
        }
    }

    double doIter(double[] betas, double[][] x, double mul) {
//        stderr.writeln("ITER:  ", betas, '\t', predictions);
        double maxRelError = 0;
        foreach(j, ref b; betas) {
            if(b == 0) {
                residuals[] = (-predictions[] + y[]);
            } else {
                residuals[] = (-predictions[] + x[j][] * b + y[]);
                predictions[] -= x[j][] * b;
            }

            double z;
            if(w.length) {
                z = 0;
                foreach(i, weight; w) {
                    z += weight * x[j][i] * residuals[i];
                }
            } else {
                residuals[] /= n;
                z = dotProduct(residuals, x[j]);
            }

            auto newB = softThresh(z, lasso * mul) /
                (weightDots[j] + ridge * mul);

            immutable absErr = abs(b - newB);
            immutable err = abs(b - newB) / max(abs(b), abs(newB));

            if(absErr > absEpsilon) {
                maxRelError = max(maxRelError, err);
            }

            b = newB;
            if(b != 0) {
                predictions[] += x[j][] * b;
            }
        }

        return maxRelError;
    }

    void toConvergence(double mul) {
        double maxRelErr = doIter(betas, x, mul);
        iter++;
        if(maxRelErr < relEpsilon) return;

        static bool absGreater(double x, double y) { return abs(x) > abs(y); }

        while(iter < maxIter) {
            size_t split = 0;
            while(split < betas.length && abs(betas[split]) > 0) {
                split++;
            }

            try {
                qsort!absGreater(betas, x, perm, weightDots);
            } catch(SortException) {
                betas[] = double.nan;
                break;
            }

            maxRelErr = double.infinity;
            for(; !(maxRelErr < relEpsilon) && split < betas.length
            && iter < maxIter; iter++) {
                maxRelErr = doIter(betas[0..split], x[0..split], mul);
            }


            maxRelErr = doIter(betas, x, mul);
            iter++;
            if(maxRelErr < relEpsilon) break;
        }

        //stderr.writefln("Converged in %s iterations for %s mult.", iter, mul);
    }

    toConvergence(1);
    try {
        qsort(perm, x, betas);
    } catch(SortException) {
        return;
    }
}

// Compute(X X')' = C.
double[][] makeC(double[][] x) {
    if(x.length == 0) return null;
    immutable n = x[0].length;
    auto c = newStack!(double[])(n);

    // Only need lower half.
    foreach(i, ref elem; c) {
        elem = newStack!double(i + 1);
        elem[] = 0;
    }

    foreach(col; x) {
        foreach(i; 0..n) {
            foreach(j; i..n) {
                c[j][i] += col[i] * col[j];
            }
        }
    }

    return c;
}

private double[][] doCTWC(double[][] c, double[] w = null) {
    auto ret = newStack!(double[])(c.length);
    foreach(i, ref elem; ret) elem = newStack!double(c.length);

    double getElem(size_t i, size_t j) {
        return (i >= j) ? c[i][j] : c[j][i];
    }

    foreach(i; 0..c.length) foreach(j; 0..i + 1) {
        ret[i][j] = 0;

        foreach(k; 0..c.length) {
            auto toAdd = getElem(i, k) * getElem(j, k);
            if(w.length) toAdd *= w[k];
            ret[i][j] += toAdd;
        }
    }

    symmetrize(ret);
    return ret;
}

// An implementation of ridge regression for large dimension.
private void ridgeLargeP
(double[] yIn, double[][] x, double lambdaIn, double[] betas, double[] w = null) {
    {
        mixin(newFrame);
        auto y = tempdup(yIn);
        auto c = makeC(x);

        // This scaling preserves numerical stability when the distrib. of
        // x has a long tail.
        auto maxElem = reduce!max(0.0, map!abs(joiner(c)));
        foreach(col; c) col[] /= maxElem;
        y[] /= (maxElem);
        auto lambda = lambdaIn / maxElem;
        auto cwc = doCTWC(c, w);

        foreach(i, row; c) foreach(j, elem; row) {
            cwc[j][i] += lambda * elem;
            if(i != j) cwc[i][j] += lambda * elem;
        }

        // Multiply c' * y.  c is symmetric, so it doesn't matter that I'm really
        // multiplying the transpose.
        if(w.length) y[] *= w[];
        auto cTy = newStack!double(y.length);
        cTy[] = 0;

        foreach(i; 0..c.length) foreach(j; 0..c.length) {
            if(i >= j) {
                cTy[j] += y[i] * c[i][j];
            } else {
                cTy[j] += y[i] * c[j][i];
            }
        }

        solve(cwc, cTy);

        // Multiply X' * csi to get answer.
        betas[] = 0;
        foreach(i, col; x) {
            betas[i] = dotProduct(cTy, col);
        }
    }

    static bool notFinite(double num) { return !isFinite(num); }

    // In a few pathological cases this algo ends up singular to machine
    // precision, and NaNs or infs up in betas.  Fall back to solving the
    // normal equations directly.
    if(filter!notFinite(betas).empty) return;
    version(unittest) stderr.writeln("SINGULAR");
    ridgeFallback(yIn, x, lambdaIn, betas, w);
}

// Used if we end up singular in ridgeLargeP.
private void ridgeFallback
(double[] yIn, double[][] x, double lambda, double[] betas, double[] w) {
    mixin(newFrame);

    // Compute xTx or xT * w * x depending on whether w is null.
    auto xTx = newStack!(double[])(x.length);
    foreach(ref row; xTx) row = newStack!double(x.length);

    foreach(i, xi; x) foreach(j, xj; x[0..i + 1]) {
        xTx[i][j] = 0;
        foreach(k; 0..xi.length) {
            immutable weight = (w.length) ? w[k] : 1.0;
            xTx[i][j] += xi[k] * weight * xj[k];
        }
    }

    symmetrize(xTx);

    double[] y;
    if(w.length) {
        y = tempdup(yIn);
        y[] *= w[];
    } else {
        y = yIn;
    }

    auto xTy = newStack!double(x.length);
    foreach(i, ref val; xTy) {
        val = dotProduct(y, x[i]);
    }

    // Add in the penalty.
    foreach(i; 0..xTx.length) {
        xTx[i][i] += lambda;
    }

    choleskySolve(xTx, xTy, betas);
}

unittest {
    // Test ridge regression.  We have three impls for all kinds of diff.
    // scenarios.  See if they all agree.  Note that the ridiculously small but
    // nonzero lasso param is to force the use of the coord descent algo.
    auto y = new double[12];
    auto x = new double[][16];
    foreach(ref elem; x) elem = new double[12];
    x[0][] = 1;
    auto gen = Random(31415);  // For random but repeatable results.

    foreach(iter; 0..1000) {
        foreach(col; x[1..$]) foreach(ref elem; col) elem = rNorm(0, 1, gen);
        foreach(ref elem; y) elem = rNorm(0, 1, gen);
        immutable ridge = uniform(0.1, 10.0, gen);

        auto normalEq = linearRegressBeta(y, x, ridge);
        auto coordDescent = linearRegressPenalized(
            y, x[1..$], double.min_normal, ridge);
        auto linalgTrick = linearRegressPenalized(y, x[1..$], 0, ridge);

        // Every once in a blue moon coordinate descent doesn't converge that
        // well.  These small errors are of no practical significance, hence
        // the wide tolerance.  However, if the direct normal equations
        // and linalg trick don't agree extremely closely, then something's
        // fundamentally wrong.
        assert(approxEqual(normalEq, coordDescent, 0.02, 1e-4), text(
            normalEq, coordDescent));
        assert(approxEqual(linalgTrick, coordDescent, 0.02, 1e-4), text(
            linalgTrick, coordDescent));
        assert(approxEqual(normalEq, linalgTrick, 1e-6, 1e-8), text(
            normalEq, linalgTrick));
    }

    // Make sure fallback and largeP algo get same answers.
    auto fb = new double[x.length];
    auto lp = new double[x.length];
    ridgeFallback(y, x, 1, fb, null);
    ridgeLargeP(y, x, 1, lp);
    assert(approxEqual(fb, lp));

    // With weights.
    auto w = randArray!uniform(y.length, 0.1, 1, gen);
    ridgeFallback(y, x, 1, fb, w);
    ridgeLargeP(y, x, 1, lp, w);
    assert(approxEqual(fb, lp));

    // Test stuff that's got some lasso in it.  Values from R's Penalized
    // package.
    y = [1.0, 2.0, 3, 4, 5, 6, 7];
    x = [[8.0, 6, 7, 5, 3, 0, 9],
         [3.0, 6, 2, 4, 3, 6, 8],
         [3.0, 1, 4, 1, 5, 9, 2],
         [2.0, 7, 1, 8, 2, 8, 1]];

    assert(approxEqual(linearRegressPenalized(y, x, 1, 0),
        [4.16316, -0.3603197, 0.6308278, 0, -0.2633263]));
    assert(approxEqual(linearRegressPenalized(y, x, 1, 3),
        [2.519590, -0.09116883, 0.38067757, 0.13122413, -0.05637939]));
    assert(approxEqual(linearRegressPenalized(y, x, 2, 0.1),
        [1.247235, 0, 0.4440735, 0.2023602, 0]));
    assert(approxEqual(linearRegressPenalized(y, x, 5, 7),
        [3.453787, 0, 0.10968736, 0.01253992, 0]));
}

/**
Computes a logistic regression using a maximum likelihood estimator
and returns the beta coefficients.  This is a generalized linear model with
the link function f(XB) = 1 / (1 + exp(XB)). This is generally used to model
the probability that a binary Y variable is 1 given a set of X variables.

For the purpose of this function, Y variables are interpreted as Booleans,
regardless of their type.  X may be either a range of ranges or a tuple of
ranges.  However, note that unlike in linearRegress, they are copied to an
array if they are not random access ranges.  Note that each value is accessed
several times, so if your range is a map to something expensive, you may
want to evaluate it eagerly.

If the last parameter passed in is a numeric value instead of a range,
it is interpreted as a ridge parameter and ridge regression is performed.  This
penalizes the L2 norm of the beta vector (in a scaled space) and results
in more parsimonious models.  It limits the usefulness of inference techniques
(p-values, confidence intervals), however, and is therefore not offered
in logisticRegres().

If no ridge parameter is passed, or equivalenty if the ridge parameter is
zero, then ordinary maximum likelihood regression is performed.

Note that, while this implementation of ridge regression was tested against
the R Design Package implementation, it uses slightly different conventions
that make the results not comparable without transformation.  dstats uses a
biased estimate of the variance to scale the beta vector penalties, while
Design uses an unbiased estimate.  Furthermore, Design penalizes by 1/2 of the
L2 norm, whereas dstats penalizes by the L2 norm.  Therefore, if n is the
sample size, and lambda is the penalty used with dstats, the proper penalty
to use in Design to get the same results is 2 * (n - 1) * lambda / n.

Also note that, as in linearRegress, repeat(1) can be used for the intercept
term.

Returns:  The beta coefficients for the regression model.

References:

http://en.wikipedia.org/wiki/Logistic_regression

http://socserv.mcmaster.ca/jfox/Courses/UCLA/logistic-regression-notes.pdf

S. Le Cessie and J. C. Van Houwelingen.  Ridge Estimators in Logistic
Regression.  Journal of the Royal Statistical Society. Series C
(Applied Statistics), Vol. 41, No. 1(1992), pp. 191-201

Frank E Harrell Jr (2009). Design: Design Package. R package version 2.3-0.
http://CRAN.R-project.org/package=Design
 */
double[] logisticRegressBeta(T, U...)(T yIn, U xRidge) {
    static if(isFloatingPoint!(U[$ - 1]) || isIntegral!(U[$ - 1])) {
        alias xRidge[$ - 1] ridge;
        alias xRidge[0..$ - 1] xIn;
    } else {
        enum double ridge = 0.0;
        alias xRidge xIn;
    }

    return logisticRegressImpl(false, ridge, yIn, xIn).betas;
}

/**
Plain old data struct to hold the results of a logistic regression.
*/
struct LogisticRes {
    /**The coefficients, one for each range in X.  These will be in the order
     * that the X ranges were passed in.*/
    double[] betas;

    /**The standard error terms of the X ranges passed in.*/
    double[] stdErr;

    /**
    The Wald lower confidence bounds of the beta terms, at the confidence level
    specificied.  (Default 0.95).*/
    double[] lowerBound;

    /**
    The Wald upper confidence bounds of the beta terms, at the confidence level
    specificied.  (Default 0.95).*/
    double[] upperBound;

    /**
    The P-value for the alternative that the corresponding beta value is
    different from zero against the null that it is equal to zero.  These
    are calculated using the Wald Test.*/
    double[] p;

    /**
    The log likelihood for the null model.
    */
    double nullLogLikelihood;

    /**
    The log likelihood for the model fit.
    */
    double logLikelihood;

    /**
    Akaike Information Criterion, which is a complexity-penalized goodness-
    of-fit score, equal to 2 * k - 2 log(L) where L is the log likelihood and
    k is the number of parameters.
    */
    double aic() const pure nothrow @property @safe {
        return 2 * (betas.length - logLikelihood);
    }

    /**
    The P-value for the model as a whole, based on the likelihood ratio test.
    The null here is that the model has no predictive value, the alternative
    is that it does have predictive value.*/
    double overallP;

    // Just used internally.
    private static string arrStr(T)(T arr) {
        return text(arr)[1..$ - 1];
    }

    /**Print out the results in the default format.*/
    string toString() {
        return "Betas:  " ~ arrStr(betas) ~ "\nLower Conf. Int.:  " ~
            arrStr(lowerBound) ~ "\nUpper Conf. Int.:  " ~ arrStr(upperBound) ~
            "\nStd. Err:  " ~ arrStr(stdErr) ~ "\nP Values:  " ~ arrStr(p) ~
            "\nNull Log Likelihood:  " ~ text(nullLogLikelihood) ~
            "\nLog Likelihood:  " ~ text(logLikelihood) ~
            "\nAIC:  " ~ text(aic) ~
            "\nOverall P:  " ~ text(overallP);
    }
}

/**
Similar to logisticRegressBeta, but returns a LogisticRes with useful stuff for
statistical inference.  If the last element of input is a floating point
number instead of a range, it is used to specify the confidence interval
calculated.  Otherwise, the default of 0.95 is used.

References:

http://en.wikipedia.org/wiki/Wald_test
http://en.wikipedia.org/wiki/Akaike_information_criterion
*/
LogisticRes logisticRegress(T, V...)(T yIn, V input) {
    return logisticRegressImpl!(T, V)(true, 0, yIn, input);
}

unittest {
    // Values from R.  Confidence intervals from confint.default().
    // R doesn't automatically calculate likelihood ratio P-value, and reports
    // deviations instead of log likelihoods.  Deviations are just
    // -2 * likelihood.
    alias approxEqual ae;  // Save typing.

    // Start with the basics, with X as a ror.
    auto y1 =  [1,   0, 0, 0, 1, 0, 0];
    auto x1 = [[1.0, 1 ,1 ,1 ,1 ,1 ,1],
              [8.0, 6, 7, 5, 3, 0, 9]];
    auto res1 = logisticRegress(y1, x1);
    assert(ae(res1.betas[0], -0.98273));
    assert(ae(res1.betas[1], 0.01219));
    assert(ae(res1.stdErr[0], 1.80803));
    assert(ae(res1.stdErr[1], 0.29291));
    assert(ae(res1.p[0], 0.587));
    assert(ae(res1.p[1], 0.967));
    assert(ae(res1.aic, 12.374));
    assert(ae(res1.logLikelihood, -0.5 * 8.3758));
    assert(ae(res1.nullLogLikelihood, -0.5 * 8.3740));
    assert(ae(res1.lowerBound[0], -4.5264052));
    assert(ae(res1.lowerBound[1], -0.5618933));
    assert(ae(res1.upperBound[0], 2.560939));
    assert(ae(res1.upperBound[1], 0.586275));

    // Use tuple.
    auto y2   = [1,0,1,1,0,1,0,0,0,1,0,1];
    auto x2_1 = [3,1,4,1,5,9,2,6,5,3,5,8];
    auto x2_2 = [2,7,1,8,2,8,1,8,2,8,4,5];
    auto res2 = logisticRegress(y2, repeat(1), x2_1, x2_2);
    assert(ae(res2.betas[0], -1.1875));
    assert(ae(res2.betas[1], 0.1021));
    assert(ae(res2.betas[2], 0.1603));
    assert(ae(res2.stdErr[0], 1.5430));
    assert(ae(res2.stdErr[1], 0.2507));
    assert(ae(res2.stdErr[2], 0.2108));
    assert(ae(res2.p[0], 0.442));
    assert(ae(res2.p[1], 0.684));
    assert(ae(res2.p[2], 0.447));
    assert(ae(res2.aic, 21.81));
    assert(ae(res2.nullLogLikelihood, -0.5 * 16.636));
    assert(ae(res2.logLikelihood, -0.5 * 15.810));
    assert(ae(res2.lowerBound[0], -4.2116584));
    assert(ae(res2.lowerBound[1], -0.3892603));
    assert(ae(res2.lowerBound[2], -0.2528110));
    assert(ae(res2.upperBound[0], 1.8366823));
    assert(ae(res2.upperBound[1], 0.5934631));
    assert(ae(res2.upperBound[2], 0.5733693));

    auto x2Intercept = [1,1,1,1,1,1,1,1,1,1,1,1];
    auto res2a = logisticRegress(y2,
        filter!"a.length"([x2Intercept, x2_1, x2_2]));
    foreach(ti, elem; res2a.tupleof) {
        assert(ae(elem, res2.tupleof[ti]));
    }

    // Use a huge range of values to test numerical stability.

    // The filter is to make y3 a non-random access range.
    auto y3 = filter!"a < 2"([1,1,1,1,0,0,0,0]);
    auto x3_1 = filter!"a > 0"([1, 1e10, 2, 2e10, 3, 3e15, 4, 4e7]);
    auto x3_2 = [1e8, 1e6, 1e7, 1e5, 1e3, 1e0, 1e9, 1e11];
    auto x3_3 = [-5e12, 5e2, 6e5, 4e3, -999999, -666, -3e10, -2e10];
    auto res3 = logisticRegress(y3, repeat(1), x3_1, x3_2, x3_3, 0.99);
    assert(ae(res3.betas[0], 1.115e0));
    assert(ae(res3.betas[1], -4.674e-15));
    assert(ae(res3.betas[2], -7.026e-9));
    assert(ae(res3.betas[3], -2.109e-12));
    assert(ae(res3.stdErr[0], 1.158));
    assert(ae(res3.stdErr[1], 2.098e-13));
    assert(ae(res3.stdErr[2], 1.878e-8));
    assert(ae(res3.stdErr[3], 4.789e-11));
    assert(ae(res3.p[0], 0.336));
    assert(ae(res3.p[1], 0.982));
    assert(ae(res3.p[2], 0.708));
    assert(ae(res3.p[3], 0.965));
    assert(ae(res3.aic, 12.544));
    assert(ae(res3.nullLogLikelihood, -0.5 * 11.0904));
    assert(ae(res3.logLikelihood, -0.5 * 4.5442));
    // Not testing confidence intervals b/c they'd be so buried in numerical
    // fuzz.


    // Test with a just plain huge dataset that R chokes for several minutes
    // on.  If you think this unittest is slow, try getting the reference
    // values from R.
    auto y4 = chain(
                take(cycle([0,0,0,0,1]), 500_000),
                take(cycle([1,1,1,1,0]), 500_000));
    auto x4_1 = iota(0, 1_000_000);
    auto x4_2 = array(map!(compose!(exp, "a / 1_000_000.0"))(x4_1));
    auto x4_3 = take(cycle([1,2,3,4,5]), 1_000_000);
    auto x4_4 = take(cycle([8,6,7,5,3,0,9]), 1_000_000);
    auto res4 = logisticRegress(y4, repeat(1), x4_1, x4_2, x4_3, x4_4, 0.99);
    assert(ae(res4.betas[0], -1.574));
    assert(ae(res4.betas[1], 5.625e-6));
    assert(ae(res4.betas[2], -7.282e-1));
    assert(ae(res4.betas[3], -4.381e-6));
    assert(ae(res4.betas[4], -8.343e-6));
    assert(ae(res4.stdErr[0], 3.693e-2));
    assert(ae(res4.stdErr[1], 7.188e-8));
    assert(ae(res4.stdErr[2], 4.208e-2));
    assert(ae(res4.stdErr[3], 1.658e-3));
    assert(ae(res4.stdErr[4], 8.164e-4));
    assert(ae(res4.p[0], 0));
    assert(ae(res4.p[1], 0));
    assert(ae(res4.p[2], 0));
    assert(ae(res4.p[3], 0.998));
    assert(ae(res4.p[4], 0.992));
    assert(ae(res4.aic, 1089339));
    assert(ae(res4.nullLogLikelihood, -0.5 * 1386294));
    assert(ae(res4.logLikelihood, -0.5 * 1089329));
    assert(ae(res4.lowerBound[0], -1.668899));
    assert(ae(res4.lowerBound[1], 5.439787e-6));
    assert(ae(res4.lowerBound[2], -0.8366273));
    assert(ae(res4.lowerBound[3], -4.27406e-3));
    assert(ae(res4.lowerBound[4], -2.111240e-3));
    assert(ae(res4.upperBound[0], -1.478623));
    assert(ae(res4.upperBound[1], 5.810089e-6));
    assert(ae(res4.upperBound[2], -6.198418e-1));
    assert(ae(res4.upperBound[3], 4.265302e-3));
    assert(ae(res4.upperBound[4], 2.084554e-3));

    // Test ridge stuff.
    auto ridge2 = logisticRegressBeta(y2, repeat(1), x2_1, x2_2, 3);
    assert(ae(ridge2[0], -0.40279319));
    assert(ae(ridge2[1], 0.03575638));
    assert(ae(ridge2[2], 0.05313875));

    auto ridge2_2 = logisticRegressBeta(y2, repeat(1), x2_1, x2_2, 2);
    assert(ae(ridge2_2[0], -0.51411490));
    assert(ae(ridge2_2[1], 0.04536590));
    assert(ae(ridge2_2[2], 0.06809964));
}

/// The logistic function used in logistic regression.
double logistic(double xb) pure nothrow @safe {
    return 1.0 / (1 + exp(-xb));
}

/**
Performs lasso (L1) and/or ridge (L2) penalized logistic regression.  Due to the
way the data is standardized, no intercept term should be included in x
(unlike logisticRegress and logisticRegressBeta).  The intercept coefficient is
implicitly included and returned in the first element of the returned array.
Usage is otherwise identical.

Note:  Setting lasso equal to zero is equivalent to performing ridge regression.
       This can also be done with logisticRegressBeta.  However, the
       logisticRegressBeta algorithm is optimized for memory efficiency and
       large samples.  This algorithm is optimized for large feature sets.

Returns:  The beta coefficients for the regression model.

References:

Friedman J, et al Pathwise coordinate optimization. Ann. Appl. Stat.
2007;2:302-332.

Goeman, J. J., L1 penalized estimation in the Cox proportional hazards model.
Biometrical Journal 52(1), 70{84.

Eilers, P., Boer, J., Van Ommen, G., Van Houwelingen, H. 2001 Classification of
microarray data with penalized logistic regression. Proceedings of SPIE.
Progress in Biomedical Optics and Images vol. 4266, pp. 187-198
*/
double[] logisticRegressPenalized(Y, X...)
(Y yIn, X xIn, double lasso, double ridge) {
    mixin(newFrame);

    static assert(!isInfinite!Y, "Can't do regression with infinite # of Y's.");
    static if(isRandomAccessRange!Y) {
        alias yIn y;
    } else {
        auto y = toBools(yIn);
    }

    static if(X.length == 1 && isRoR!X) {
        enum bool tupleMode = false;
        static if(isForwardRange!X) {
            auto x = toRandomAccessRoR(y.length, xIn);
        } else {
            auto x = toRandomAccessRoR(y.length, tempdup(xIn));
        }
    } else {
        enum bool tupleMode = true;
        auto x = toRandomAccessTuple(xIn).expand;
    }

    auto betas = new double[x.length + 1];
    if(y.length >= x.length && lasso == 0) {
        // Add intercept term.
        static if(tupleMode) {
            doMLENewton(betas, (double[]).init, ridge, y, repeat(1), x);
        } else static if(is(x == double[][])) {
            auto xInt = newStack!(double[])(betas.length);
            xInt[1..$] = x[];
            xInt[0] = tempdup(replicate(1.0, y.length));
            doMLENewton(betas, (double[]).init, ridge, y, xInt);
        } else {
            // No choice but to dup the whole thing.
            auto xInt = newStack!(double[])(betas.length);
            xInt[0] = tempdup(replicate(1.0, y.length));

            foreach(i, ref arr; xInt[1..$]) {
                arr = tempdup(
                    map!(to!double)(take(x[i], y.length))
                );
            }

            doMLENewton(betas, (double[]).init, ridge, y, xInt);
        }

    } else {
        logisticRegressPenalizedImpl(betas, lasso, ridge, y, x);
    }

    return betas;
}

unittest {
    // Test ridge regression.  We have three impls for all kinds of diff.
    // scenarios.  See if they all agree.  Note that the ridiculously small but
    // nonzero lasso param is to force the use of the coord descent algo.
    auto y = new bool[12];
    auto x = new double[][16];
    foreach(ref elem; x) elem = new double[12];
    x[0][] = 1;
    auto gen = Random(31415);  // For random but repeatable results.

    foreach(iter; 0..1000) {
        foreach(col; x[1..$]) foreach(ref elem; col) elem = rNorm(0, 1, gen);

        // Nothing will converge if y is all true's or all false's.
        size_t trueCount;
        do {
            foreach(ref elem; y) elem = cast(bool) rBernoulli(0.5, gen);
            trueCount = count!"a"(y);
        } while(trueCount == 0 || trueCount == y.length);

        immutable ridge = uniform(0.1, 10.0, gen);

        auto normalEq = logisticRegressBeta(y, x, ridge);
        auto coordDescent = logisticRegressPenalized(
            y, x[1..$], double.min_normal, ridge);
        auto linalgTrick = logisticRegressPenalized(y, x[1..$], 0, ridge);

        // Every once in a blue moon coordinate descent doesn't converge that
        // well.  These small errors are of no practical significance, hence
        // the wide tolerance.  However, if the direct normal equations
        // and linalg trick don't agree extremely closely, then something's
        // fundamentally wrong.
        assert(approxEqual(normalEq, coordDescent, 0.02, 1e-4), text(
            normalEq, coordDescent));
        assert(approxEqual(linalgTrick, coordDescent, 0.02, 1e-4), text(
            linalgTrick, coordDescent));
        assert(approxEqual(normalEq, linalgTrick, 1e-6, 1e-8), text(
            normalEq, linalgTrick));
    }

    assert(approxEqual(logisticRegressBeta(y, x[0], x[1], x[2]),
        logisticRegressPenalized(y, x[1], x[2], 0, 0)));
    assert(approxEqual(logisticRegressBeta(y, [x[0], x[1], x[2]]),
        logisticRegressPenalized(y, [x[1], x[2]], 0, 0)));
    assert(approxEqual(logisticRegressBeta(y, [x[0], x[1], x[2]]),
        logisticRegressPenalized(y,
        [to!(float[])(x[1]), to!(float[])(x[2])], 0, 0)));

    // Make sure the adding intercept stuff is right for the Newton path.
    //assert(logisticRegressBeta(x[0], x[1], x[2]) ==

    // Test stuff that's got some lasso in it.  Values from R's Penalized
    // package.
    y = [1, 0, 0, 1, 1, 1, 0];
    x = [[8.0, 6, 7, 5, 3, 0, 9],
         [3.0, 6, 2, 4, 3, 6, 8],
         [3.0, 1, 4, 1, 5, 9, 2],
         [2.0, 7, 1, 8, 2, 8, 1]];

    // Values from R's Penalized package.  Note that it uses a convention for
    // the ridge parameter such that Penalized ridge = 2 * dstats ridge.
    assert(approxEqual(logisticRegressPenalized(y, x, 1, 0),
        [1.642080, -0.22086515, -0.02587546,  0.00000000, 0.00000000 ]));
    assert(approxEqual(logisticRegressPenalized(y, x, 1, 3),
        [0.5153373, -0.04278257, -0.00888014,  0.01316831,  0.00000000]));
    assert(approxEqual(logisticRegressPenalized(y, x, 2, 0.1),
        [0.2876821, 0, 0., 0., 0]));
    assert(approxEqual(logisticRegressPenalized(y, x, 1.2, 7),
        [0.367613 , -0.017227631, 0.000000000, 0.003875104, 0.000000000]));
}

// Scheduled for deprecation.  This was a terrble name choice.
alias logistic inverseLogit;

private:
double absMax(double a, double b) {
    return max(abs(a), abs(b));
}

LogisticRes logisticRegressImpl(T, V...)
(bool inference, double ridge, T yIn, V input) {
    mixin(newFrame);

    static if(isFloatingPoint!(V[$ - 1])) {
        alias input[$ - 1] conf;
        alias V[0..$ - 1] U;
        alias input[0..$ - 1] xIn;
        enforceConfidence(conf);
    } else {
        alias V U;
        alias input xIn;
        enum conf = 0.95;
    }

    static assert(!isInfinite!T, "Can't do regression with infinite # of Y's.");
    static if(isRandomAccessRange!T) {
        alias yIn y;
    } else {
        auto y = toBools(yIn);
    }

    static if(U.length == 1 && isRoR!U) {
        static if(isForwardRange!U) {
            auto x = toRandomAccessRoR(y.length, xIn);
        } else {
            auto x = toRandomAccessRoR(y.length, tempdup(xIn));
        }
    } else {
        auto x = toRandomAccessTuple(xIn).expand;
    }

    typeof(return) ret;
    ret.betas.length = x.length;
    if(inference) ret.stdErr.length = x.length;
    ret.logLikelihood = doMLENewton(ret.betas, ret.stdErr, ridge, y, x);

    if(!inference) return ret;

    static bool hasNaNs(R)(R range) {
        return !filter!isNaN(range).empty;
    }

    if(isNaN(ret.logLikelihood) || hasNaNs(ret.betas) || hasNaNs(ret.stdErr)) {
        // Then we didn't converge or our data was defective.
        return ret;
    }

    ret.nullLogLikelihood = .priorLikelihood(y);
    double lratio = ret.logLikelihood - ret.nullLogLikelihood;

    // Compensate for numerical fuzz.
    if(lratio < 0 && lratio > -1e-5) lratio = 0;
    if(lratio > 0) {
        ret.overallP = chiSquareCDFR(2 * lratio, x.length - 1);
    }

    ret.p.length = x.length;
    ret.lowerBound.length = x.length;
    ret.upperBound.length = x.length;
    immutable nDev = -invNormalCDF((1 - conf) / 2);
    foreach(i; 0..x.length) {
        ret.p[i] = 2 * normalCDF(-abs(ret.betas[i]) / ret.stdErr[i]);
        ret.lowerBound[i] = ret.betas[i] - nDev * ret.stdErr[i];
        ret.upperBound[i] = ret.betas[i] + nDev * ret.stdErr[i];
    }

    return ret;
}

private void logisticRegressPenalizedImpl(Y, X...)
(double[] betas, double lasso, double ridge, Y y, X xIn) {
    static if(isRoR!(X[0]) && X.length == 1) {
        alias xIn[0] x;
    } else {
        alias xIn x;
    }

    auto ps = newStack!double(y.length);
    betas[] = 0;

    mixin(newFrame);
    auto betasRaw = tempdup(betas[1..$]);

    double oldLikelihood = -double.infinity;
    double oldPenalty2 = double.infinity;
    double oldPenalty1 = double.infinity;
    enum eps = 1e-6;
    enum maxIter = 1000;

    auto weights = newStack!double(y.length);
    auto z = newStack!double(y.length);
    auto xMeans = newStack!double(x.length);
    auto xSds = newStack!double(x.length);
    double zMean = 0;
    auto xCenterScale = newStack!(double[])(x.length);
    foreach(ref col; xCenterScale) col = newStack!double(y.length);

    // Puts x in xCenterScale, with weighted mean subtracted and weighted
    // biased stdev divided.  Also standardizes z similarly.
    //
    // Returns:  true if successful, false if weightSum is so small that the
    //           algorithm has converged for all practical purposes.
    bool doCenterScale() {
        immutable weightSum = sum(weights);
        if(weightSum < eps) return false;

        xMeans[] = 0;
        zMean = 0;
        xSds[] = 0;

        // Do copying.
        foreach(i, col; x) {
            copy(take(col, y.length), xCenterScale[i]);
        }

        // Compute weighted means.
        foreach(j, col; xCenterScale) {
            foreach(i, w; weights) {
                xMeans[j] += w * col[i];
            }
        }

        foreach(i, w; weights) {
            zMean += w * z[i];
        }

        xMeans[] /= weightSum;
        zMean /= weightSum;
        z[] -= zMean;
        foreach(i, col; xCenterScale) col[] -= xMeans[i];

        // Compute biased stdevs.
        foreach(j, ref sd; xSds) {
            sd = sqrt(meanStdev(xCenterScale[j]).mse);
        }

        foreach(i, col; xCenterScale) col[] /= xSds[i];
        return true;
    }

    // Rescales the beta coefficients to undo the effects of standardizing.
    void rescaleBetas() {
        betas[0] = zMean;
        foreach(i, b; betasRaw) {
            betas[i + 1] = b / xSds[i];
            betas[0] -= betas[i + 1] * xMeans[i];
        }
    }

    foreach(iter; 0..maxIter) {
        evalPs(betas[0], ps, betas[1..$], x);
        immutable lh = logLikelihood(ps, y);
        immutable penalty2 = ridge * reduce!"a + b * b"(0.0, betas);
        immutable penalty1 = lasso * reduce!"a + (b < 0) ? -b : b"(0.0, betas);

        if((lh - oldLikelihood) / (absMax(lh, oldLikelihood) + 0.1) < eps
        && (oldPenalty2 - penalty2) / (absMax(penalty2, oldPenalty2) + 0.1) < eps
        && (oldPenalty1 - penalty1) / (absMax(penalty1, oldPenalty1) + 0.1) < eps) {
            return;
        } else if(isNaN(lh) || isNaN(penalty2) || isNaN(penalty1)) {
            betas[] = double.nan;
            return;
        }

        oldPenalty2 = penalty2;
        oldPenalty1 = penalty1;
        oldLikelihood = lh;

        foreach(i, ref w; weights) {
            w = (ps[i] * (1 - ps[i]));
        }

        z[] = betas[0];
        foreach(i, col; x) {
            static if(is(typeof(col) : const(double)[])) {
                z[] += col[] * betas[i + 1];
            } else {
                foreach(j, ref elem; z) {
                    elem += col[j] + betas[i + 1];
                }
            }
        }

        foreach(i, w; weights) {
            if(w == 0){
                z[i] = 0;
            } else {
                immutable double yi = (y[i] == 0) ? 0.0 : 1.0;
                z[i] += (yi - ps[i]) / w;
            }
        }

        immutable centerScaleRes = doCenterScale();

        // If this is false then weightSum is so small that all probabilities
        // are for all practical purposes either 0 or 1.  We can declare
        // convergence and go home.
        if(!centerScaleRes) return;

        if(lasso > 0) {
            // Correct for different conventions in defining ridge params
            // so all functions get the same answer.
            immutable ridgeCorrected = ridge * 2.0;
            coordDescent(z, xCenterScale, betasRaw,
                lasso, ridgeCorrected, weights);
        } else {
            // Correct for different conventions in defining ridge params
            // so all functions get the same answer.
            immutable ridgeCorrected = ridge * 2.0;
            ridgeLargeP(z, xCenterScale, ridgeCorrected, betasRaw, weights);
        }

        rescaleBetas();
    }

    immutable lh = logLikelihood(ps, y);
    immutable penalty2 = ridge * reduce!"a + b * b"(0.0, betas);
    immutable penalty1 = lasso * reduce!"a + (b < 0) ? -b : b"(0.0, betas);

    if((lh - oldLikelihood) / (absMax(lh, oldLikelihood) + 0.1) < eps
    && (oldPenalty2 - penalty2) / (absMax(penalty2, oldPenalty2) + 0.1) < eps
    && (oldPenalty1 - penalty1) / (absMax(penalty1, oldPenalty1) + 0.1) < eps) {
        return;
    } else {
        // If we got here, we haven't converged.  Return NaNs instead of bogus
        // values.