辅导STAT - 611、讲解C/C++程序、辅导Gram-Schmidt、C/C++设计辅导

STAT - 611 Assignment 6

Due Tuesday Dec 4th, 2018, in class.

Remember (from syllabus): Solutions must be neatly word-processed (using Latex or a word processor)

and stapled. Please annotate your work with brief, clear sentences explaining your approach and interpreting

your results (you are not expected to write full-blown data analysis reports however). For assignments

involving mathematical manipulations, students can write answers by hand provided penmanship is neat;

illegible answers will be marked as incorrect.

*All these questions require you write computer programs. Please write a report with your code, and an

explanation of what youve done. Additionally upload the files containing your computer code to Canvas.

Include instructions detailing how to use them in your report.

1. The files “QR.cpp” and “QR.h” (available in Canvas) contain an implementation of the QR decomposition

via the modified Gram-Schmidt procedure. The prototype of this function, in the file QR.h,

is:

void QR(double *At_raw, double *Rt_raw, int nrow_A, int ncol_A);

Where the arguments are

At raw: The data of the matrix A, of dimensions I ×J, to decompose. These data are represented

as a flat array in column-major-order (or equivalently, a flat array containing the transpose of

such matrix in row-major-order).

Rt raw: space allocated externally for the R matrix, as a flat array (of length J 2).

nrow A: the number of rows of A (I).

ncol A: the number of columns of A (J).

The resulting matrix Q substitutes the data in At raw and is returned in column-major-order (or as

a flat representation of QT

in row-major-order). The matrix R is returned in the space pointed by

Rt raw, again as a flat array in column-major-order.

Your task is:

(a) Write a C function callable from R that invokes the QR() function, so that we can use it from R.

The prototype should be

extern"C" void QR_R(double *At_raw, double *Rt_raw, int *I, int *J);

Write your function at the bottom of the QR.cpp file, and the prototype in the file QR.h. Note

that by including extern"C" in the prototype of the function you don’t need to include the “#ifdef

cplusplus...” clause, but compilation could fail if you try to compile it with a C compiler (not

C++). Compile the source code so that the binary file is callable from R.

Load your module into R and write a wrapper function of the form

my_QR <- function(X){

#your code here

}

where X is an R matrix. The function should return a named list with the matrices Q and R.

Test your interface by decomposing a matrix X and testing that is the case X = QR and that

QtQ = I. Comment your results (do you note anything strange in the product QT Q? Explain)

The following functions/operators might be useful:

1

t(a): Compute the transpose of the matrix a.

a %*% b: multiply matrices a and b.

ncol(a), nrow(b): get the number of columns and rows of matrix a, respectively.

Check the help files of those functions for more details.

(b) Write an R function for estimating the coefficients of a multiple linear regression

y = Xβ + 

using your my QR R function. Your function must have the form:

lin_regr <- function(y, X){

....

}

where X is a n × p design matrix, and y is a vector containing n response observations. Your

function must return a vector containing the estimated regression coefficients.

Test your work with the following code:

> set.seed(1)

> X <- matrix(rnorm(3*7),ncol = 3)

> y <- X %*% 1:3 + rnorm(7)

> lin_regr(y, X)

and compare your output to R’s implementation of least squares regression (the “-1” is to avoid

fitting an intercept term):

> lm(y~X - 1)