// Multivariate linear regression
function lm() {
  
  //make science.lin function easier to reference
  var m = science.lin.multiply,
      t = science.lin.transpose,
      inv = science.lin.inverse;
      
  var _y, _X, _data;
  
  function model() {
    var y = vectorToMatrix(_data.map(_y));  // dependent variable (n x 1 matrix)
    var X = _data.map(_X);                  // independent variables (n x k matrix)
    var n = X.length;                       // n = # of observations
    var k = X[0].length;                    // k = # of model parameters
    var invXX = inv(m(t(X), X));            // (X'X)^1  (k x k matrix)
    var B = m(invXX, m(t(X), y));           // model parameters (k x 1 matrix)
    var y_pred = m(X, B);                   // in-sample model predictions
    var e = d3.range(n).map(function(i) {   // prediction error
      return [y[i][0] - y_pred[i][0]];
    });
    
    // Variance-covariance matrix using scaled White's formula
    // (see chapter 4 from http://www.ssc.wisc.edu/~bhansen/econometrics/)
    var u = X.map(function(x, i) {
      return x.map(function(d) { return d * e[i][0]; });
    });
    var whites_scaler = n/(n-k);
    var V = m(invXX, m(m(t(u), u), invXX))
      .map(function(row) {
        return row.map(function(d) { return whites_scaler * d; });
      });
    
    // parameter standard errors
    var se = diag(V).map(function(d) { return [Math.sqrt(d[0])]; }); 
     
    // t-statistics
    var tstat = B.map(function(b, i) { return [b[0] / se[i][0]]; });
    
    return {
      B: B,
      se: se,
      V: V,
      t: tstat,
      getRegressionInterval: getRegressionInterval
    };
    
    // get asymptotic 95% confidence regression interval 
    // at a given data point
    function getRegressionInterval(d) {
      var x = [_X(d)]; // 1 x k
      var y_pred = m(x, B)[0][0];
      var interval = 1.96 * Math.sqrt(m(m(x, V), t(x))[0][0]);
      return {
        y_lower: y_pred - interval,
        y_pred: y_pred,
        y_upper: y_pred + interval
      };
    }
  }
  
  model.y = function(_) {
    if (!arguments.length) return _y;
    _y = _;
    return model;
  };
  
  model.X = function(_) {
    if (!arguments.length) return _X;
    _X = _;
    return model;
  };
  
  model.data = function(_) {
    if (!arguments.length) return _data;
    _data = _;
    return model;
  };
  
  return model;
  
  // convert to/from "array of numbers" from/to "1D matrix"
  function vectorToMatrix(_) { return _.map(function(d) { return [d]; }); }
  function matrixToVector(_) { return _.map(function(d) { return d[0]; }); }
  
  // create matrix of size r x c with all 1's
  function ones(r, c) {
    return d3.range(r).map(function() {
      return d3.range(c).map(function() { return 1; });
    });
  }
  
  // get diagonal vector from square matrix
  function diag(x) { return x.map(function(d, j) { return [d[j]]; }); }
}