[GenABEL-dev] [Genabel-commits] r1698 - pkg/ProbABEL/src
L.C. Karssen
lennart at karssen.org
Sun Apr 27 21:48:07 CEST 2014
Thanks for splitting this into separate functions, Maarten.
Could you try to add some basic Doxygen documentation for these
functions? That would be a great help (even though the names of the
functions are already explaining a lot) towards getting a
well-documented code base for ProbABEL.
Lennart.
On 24-04-14 20:50, noreply at r-forge.r-project.org wrote:
> Author: maartenk
> Date: 2014-04-24 20:50:51 +0200 (Thu, 24 Apr 2014)
> New Revision: 1698
>
> Modified:
> pkg/ProbABEL/src/reg1.cpp
> pkg/ProbABEL/src/reg1.h
> Log:
> -refactored linear_reg::estimate a bit to make it more readable.
>
> Modified: pkg/ProbABEL/src/reg1.cpp
> ===================================================================
> --- pkg/ProbABEL/src/reg1.cpp 2014-04-24 16:50:53 UTC (rev 1697)
> +++ pkg/ProbABEL/src/reg1.cpp 2014-04-24 18:50:51 UTC (rev 1698)
> @@ -327,6 +327,14 @@
> sigma2 = (Y - tXW * beta_vec).squaredNorm();
> beta.data = beta_vec;
> }
> +void linear_reg::LeastSquaredRegression(mematrix<double> X,LDLT<MatrixXd>& Ch) {
> + int m = X.ncol;
> + MatrixXd txx = MatrixXd(m, m).setZero().selfadjointView<Lower>().rankUpdate(
> + X.data.adjoint());
> + Ch = LDLT < MatrixXd > (txx.selfadjointView<Lower>());
> + beta.data = Ch.solve(X.data.adjoint() * reg_data.Y.data);
> + sigma2 = (reg_data.Y.data - (X.data * beta.data)).squaredNorm();
> +}
>
> void linear_reg::logLikelihood(const mematrix<double>& X) {
> /*
> @@ -364,6 +372,27 @@
> }
>
>
> +
> +void linear_reg::RobustSEandCovariance(mematrix<double> X, mematrix<double> robust_sigma2,
> + MatrixXd tXX_inv, int offset) {
> + MatrixXd Xresiduals = X.data.array().colwise()
> + * residuals.data.col(0).array();
> + MatrixXd XbyR =
> + MatrixXd(X.ncol, X.ncol).setZero().selfadjointView<Lower>().rankUpdate(
> + Xresiduals.adjoint());
> + robust_sigma2.data = tXX_inv * XbyR * tXX_inv;
> + sebeta.data = robust_sigma2.data.diagonal().array().sqrt();
> + covariance.data =
> + robust_sigma2.data.bottomLeftCorner(offset, offset).diagonal();
> +}
> +
> +void linear_reg::PlainSEandCovariance(double sigma2_internal, MatrixXd tXX_inv,
> + int offset) {
> + sebeta.data = (sigma2_internal * tXX_inv.diagonal().array()).sqrt();
> + covariance.data = sigma2_internal
> + * tXX_inv.bottomLeftCorner(offset, offset).diagonal().array();
> +}
> +
> void linear_reg::estimate(int verbose, double tol_chol,
> int model, int interaction, int ngpreds, masked_matrix& invvarmatrixin,
> int robust, int nullmodel) {
> @@ -415,13 +444,6 @@
> {
> //retrieve masked data W
> invvarmatrixin.update_mask(reg_data.masked_data);
> -
> - // This regression is Weighted Least Square: used for mmscore :
> - // FLOPS count are calculated for 3*1000 matrix as follow:
> - //C=AB (m X n matrix A and n x P matrix B)
> - //flops=mp(2n-1) (when n is big enough flops=mpn2)
> - //Oct 26, 2009
> -
> mmscore_regression(X, invvarmatrixin, Ch);
> double N = X.nrow;
> //sigma2_internal = sigma2 / (N - static_cast<double>(length_beta));
> @@ -434,14 +456,7 @@
> else // NO mm-score regression : normal least square regression
> {
>
> - int m = X.ncol;
> - MatrixXd txx = MatrixXd(m, m).setZero().selfadjointView<Lower>().\
> - rankUpdate(X.data.adjoint());
> - Ch = LDLT <MatrixXd>(txx.selfadjointView<Lower>());
> - beta.data = Ch.solve(X.data.adjoint() * reg_data.Y.data);
> - sigma2 = (reg_data.Y.data - (X.data * beta.data)).squaredNorm();
> -
> -
> + LeastSquaredRegression(X,Ch);
> double N = static_cast<double>(X.nrow);
> double P = static_cast<double>(length_beta);
> sigma2_internal = sigma2 / (N - P);
> @@ -468,43 +483,22 @@
> Identity(length_beta, length_beta));
>
> mematrix<double> robust_sigma2(X.ncol, X.ncol);
> - if (robust)
> - {
> - MatrixXd Xresiduals = X.data.array().colwise()\
> - *residuals.data.col(0).array();
> - MatrixXd XbyR = MatrixXd(X.ncol, X.ncol).setZero()\
> - .selfadjointView<Lower>().rankUpdate(Xresiduals.adjoint());
> - robust_sigma2.data = tXX_inv * XbyR * tXX_inv;
> - }
> - //cout << "estimate 0\n";
> - if (robust)
> - {
> - sebeta.data = robust_sigma2.data.diagonal().array().sqrt();
> - }
> - else
> - {
> - sebeta.data =
> - (sigma2_internal
> - * tXX_inv.diagonal().array()).sqrt();
> - }
> - int offset = X.ncol- 1;
> - //if additive and interaction and 2 predictors and more then 2 betas
>
> - if (model == 0 && interaction != 0 && ngpreds == 2 && length_beta > 2){
> - offset = X.ncol - 2;
> - }
>
> + int offset = X.ncol- 1;
> + //if additive and interaction and 2 predictors and more then 2 betas
> + if (model == 0 && interaction != 0 && ngpreds == 2 && length_beta > 2){
> + offset = X.ncol - 2;
> + }
> +
> if (robust)
> {
> - covariance.data = robust_sigma2.data.bottomLeftCorner(
> - offset, offset).diagonal();
> + RobustSEandCovariance(X, robust_sigma2, tXX_inv, offset);
> }
> else
> {
> - covariance.data = sigma2_internal
> - * tXX_inv.bottomLeftCorner(offset,
> - offset).diagonal().array();
> - }
> + PlainSEandCovariance(sigma2_internal, tXX_inv, offset);
> + }
>
> }
>
>
> Modified: pkg/ProbABEL/src/reg1.h
> ===================================================================
> --- pkg/ProbABEL/src/reg1.h 2014-04-24 16:50:53 UTC (rev 1697)
> +++ pkg/ProbABEL/src/reg1.h 2014-04-24 18:50:51 UTC (rev 1698)
> @@ -104,6 +104,11 @@
> void mmscore_regression(const mematrix<double>& X,
> const masked_matrix& W_masked, LDLT<MatrixXd>& Ch);
> void logLikelihood(const mematrix<double>& X);
> + void LeastSquaredRegression(mematrix<double> X,LDLT<MatrixXd>& Ch);
> + void RobustSEandCovariance(mematrix<double> X, mematrix<double> robust_sigma2,
> + MatrixXd tXX_inv, int offset);
> + void PlainSEandCovariance(double sigma2_internal, MatrixXd tXX_inv,
> + int offset);
> };
>
> class logistic_reg: public base_reg {
>
> _______________________________________________
> Genabel-commits mailing list
> Genabel-commits at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-commits
>
--
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
L.C. Karssen
Utrecht
The Netherlands
lennart at karssen.org
http://blog.karssen.org
GPG key ID: A88F554A
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 213 bytes
Desc: OpenPGP digital signature
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20140427/7ea3adb9/attachment.sig>
More information about the genabel-devel
mailing list