Fast implementation of Gaussian log-likelihood for regression

[alashworth]

Issue by bgoodri Thursday Feb 19, 2015 at 02:46 GMT Originally opened as https://github.com/stan-dev/stan/issues/1305

---

@rtrangucci Can you implement something like the function below as a native Stan function, which should be equivalent to

y ~ normal(X * beta, sigma);

but faster with big and / or sparse X matrices. We should use it in the .stan program for a stan_lm() R function.

/**
 * Increments the log_prob with the logarithm of a multivariate normal 
 * distribution with scalar standard deviation of the errors
 * @param beta vector of typically unknown coefficients
 * @param b vector of OLS coefficients that should be precomputed
 * @param middle matrix typically precomputed as crossprod(X)
 * @param SSR positive precomputed value of the sum of squared OLS residuals
 * @param sigma positive unknown value for the standard deviation of the errors
 * @param N integer equal to the number of observations
 * @param constant can be 0 or 0.5 * N * log(2 * pi())
 */
void mvn_ols_lp(vector beta, vector b, matrix middle, 
                real SSR, real sigma, int N, real constant) {
  increment_log_prob(-0.5 * (quad_form(middle, beta - b) + SSR) / sigma^2 - 
                     N * log(sigma) - constant);
}

[alashworth]

Comment by bgoodri Tuesday Feb 24, 2015 at 05:27 GMT

---

This is related to #1226

[alashworth]

Comment by bgoodri Tuesday Feb 24, 2015 at 05:29 GMT

---

And to #1203