R语言系列：简单作图

以下函数只为满足常用的若干作图需求。

基本作图：
plot(x)、plot(x, y)    #散点图，最多两个变量
    #可使用参数type生成不同的效果图。常用‘l’、‘o’、‘h’，分别为折线图，点线图，垂线图。
    #'s'和'S'是折线图，前者是先水平后垂直，后者是先垂直后水平；'n'是不显示，用于画空白图
    #若需对三个以上变量两两作图，可先合并在一个数据框，再对数据框使用plot
    #例：r=data.frame(x,y,z); plot(r);
boxplot(x)、boxplot(x,y,z,...)    #箱式图，可使用多个参数将多个箱式图做在一起
coplot(x~y|z)    #在z的每个值或每个区间上做x与y的散点图
pairs(x)    #当x为矩阵时，做x各列之间的散点图
hist(x, freq)    #直方图。参数freq默认为TRUE，根据频数作图；若为FALSE，则据构成比（总和为1）作图。
barplot(table(x,y), beside=FALSE)    #对定性变量x、y做条图，默认堆积条图，使用beside=T则为并列条图。
qqnorm(x)    #QQ图（正态分位数-分位数图）
qqplot(x,y)    #y对x的分位数-分位数图
pie(table(x))    #对定性变量x做饼图
arrows    #在两点之间画箭头线，箭头与线段之间的夹角可调
segments    #在两点之间画线段

加参考线：    #低级绘图参数，直接在原图上修改
lines(x)、lines(x, y)    #添加折线
abline(lm(y~x))    #添加y对x的回归直线
abline(a,b)    #a为截距，b为斜率
abline(v=)    #添加垂直线
abline(h=)    #添加水平线

添加点：
points(x,y)    #低级绘图参数，直接在原图上修改

画曲线： 
1、curve(expr, from=0, to=1, n=101, add=FALSE)    
    例：curve(qnorm); curve(log); curve(x^2)    #第一个参数可以是函数，也可以是含x的表达式。 
    例：curve(x^2, 0, 100);    #from和to规定表达式或函数的自变量范围，默认0～1。 
    n：为自变量范围内取多少个点进行描图，n越大，曲线越平滑 
    add：当取值为TURE时，该曲线添加于前一图上。当两图的自变量取值范围不重合时，无法叠加。   
2、plot(fun, from=0, to=1) 
    用法和curve相似，但第一个参数不能使用含x的表达式。 

    注意在plot中尽量避免使用add参数。

使用参数pch设定点的样式
1、数字0～20，表示21种不同的符号
    curve(x/20,0,21,lty=0); for (i in 0:20) {points(i,0.5,pch=i);}
2、8个字符，表示8种点的样式
    curve(x/7,0,7,lty=0); chr=c('*','.','o','O','0','+','-','|');
    for (i in 0:7) {points(i,0.5,pch=chr[i+1])}
3、数字21～25，表示可以填充背景色的五种符号，需和参数bg联用
    curve(x/4,0,4,lty=0); for (i in 0:4) {points(i,0.5,pch=21+i,bg=sample(colors(),1))}

使用参数col设定颜色
    例：plot(dnorm, col='red')
    使用sort(colors())可以查看所有已命名的颜色   

使用参数lty设定线型
1、简单取值可以使用0～6
    curve(x*8,0,1,lty=0); for (i in 1:16) {abline(h=i/2,lty=i);}
    #0表示空白，随取值增大，6种线型循环出现
2、使用十六进制数字组成的字符串
    长度只能是偶数位，最长8位；奇数位为表示实线长度，偶数位表示空白长度
    注意：需要引号；不能有0。
    curve(x^2,0,1,lty='32';)    #3单位实线和2单位空白循环
    curve(x^2,0,1,lty='32AA';)    #3单位实线、2单位空白、10单位实线、10单位空白循环

使用参数bty设定边框
    bty='o'    #默认显示四条边框
    bty='l'    #不显示右上两条边框，例：curve(x^2,0,1,bty='l')
    bty='n'    #不显示边框

使用参数xlab和ylab修改坐标轴的意义
    curve(x*5,0,1,xlab='',ylab='',lty=0)    #不显示坐标轴的符号

使用xlim和ylim限定坐标轴的范围
    plot(rnorm,xlim=c(0,100),ylim=c(0,1000),lty=0)

使用title加图的标题
    title(main = NULL, sub = NULL, xlab = NULL, ylab = NULL)
    #main是写在图上方的标题，sub是写在图下方的标题

更多参数请查看帮助：?par

添加图例
1、图例在作图区域内
     legend(x,y,legend,pch,col,lty,ncol=1,bty='o')
          #x,y为图例左上角的坐标，也可以使用以下字符来标记特殊位置："bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center"
          #legend是图例的文字，一般是一个字符向量
          #pch,col,lty是图中所用的样式
          #ncol是图例的列数，bty指明图例的边框显示
    例：
    curve(x*8,0,1,lty=0); for (i in 1:6) {abline(h=i/2,lty=i);}
    legend(0,8,legend=letters[c(1:6)],lty=1:6,bty='n');
2、图例在作图区域外
    需使用par修改参数mar和pty。同时还需要修改参数xpd。
    默认mar=c(5,4,4,2)+0.1，即图形下左上右四个边界的宽度分别为5.1,4.1,4.1,2.1厘米。
    默认pty='m'，即最大化作图区域，取值为's'则限制作图区域为方形。
    默认xpd=FALSE，即不允许在作图区域外作图，改为TRUE即可。
    例：
    op=par();     #保存par的原值
    par(mar=c(5,4,4,5),pty='s');    #准备在图形右边添加图例
    curve(x*8,0,1,lty=0); for (i in 1:16) {abline(h=i/2,lty=i);}
    legend(1.02,4,legend=letters[c(1:6)],lty=1:6,bty='n',xpd=T);
    par(op);    #恢复par的原值

在任意位置添加文字
    curve(x*1,0,1,lty=0)    #做一张空白图
    text(0.5,0.5,'中心')    #在(0.5,0.5)标上‘中心’
    text(locator(1),'任意')    #在鼠标点击的位置标上‘任意’
        locator(n)    #返回鼠标点击n次的坐标
    text(locator(1),'任意',cex=2)    #利用cex参数将字号放大一倍
    text(locator(1),expression(sqrt(x)));    #标示公式，不能使用font和vfont修改字体
    text(locator(1),'x',vfont=c('serif','italic'))    #设定字体
        #字体设定参看：?Hershey。vfont的值为(typeface, fontindex)。
    注：text只能在作图区域内添加，而利用mtext可以在外面添加

做一个空白图（什么都没有）
    plot(rnorm,lty=0,bty='n',xaxt='n',yaxt='n',xlab='',ylab='',xlim=c(-1,1),ylim=c(-1,1));

曲线拟合：（线性回归方法：lm）
1、x排序
2、求线性回归方程并赋予一个新变量
    z=lm(y~x+I(x^2)+...)
3、plot(x,y)    #做y对x的散点图
4、lines(x,fitted(z))    #添加拟合值对x的散点图并连线

曲线拟合：（曲线回归，nls）
lm是将曲线直线化再做回归，nls是直接拟合曲线。
需要三个条件：曲线方程、数据位置、系数的估计值。
如果曲线方程比较复杂，可以先命名一个自定义函数。
例：
    f=function(x1, x2, a, b) {a+x1+x2^b};
    result=nls(x$y~f(x$x1, x$x2, a, b), data=x, start=list(a=1, b=2));
        #x可以是数据框或列表，但不能是矩阵
        #对系数的估计要尽量接近真实值，如果相差太远会报错：“奇异梯度”
    summary(result);    #结果包含对系数的估计和p值
根据估计的系数直接在散点图上使用lines加曲线即可。

曲线拟合：（局部回归）
lowess(x, y=NULL, f = 2/3, iter = 3)
    #可以只包含x，也可使用x、y两个变量
    #f为窗宽参数，越大越平滑
    #iter为迭代次数，越大计算越慢
loess(y~x, data, span=0.75, degree=2)
    #data为包含x、y的数据集；span为窗宽参数
    #degree默认为二次回归
    #该方法计算1000个数据点约占10M内存
举例：
x=seq(0, 10, 0.1); y=sin(x)+rnorm(101)    #x的值必须排序
plot(x,y);    #做散点图
lines(lowess(x,y));    #利用lowess做回归曲线
lines(x,predict(loess(y~x)));    #利用loess做回归曲线，predict是取回归预测值
z=loess(y~x); lines(x, z$fit);    #利用loess做回归曲线的另一种做法

核密度估计曲线
1、hist(x, freq=FALSE)    #根据构成比做直方图
2、核密度估计
density(x,
    bw='nrd0',    #设置窗宽，默认为‘nrd0’（只是为了兼容，并不是推荐数值），可尝试不同数字选择最合适的。
    kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine"),
    #选择列表中的一种计算方法，默认为第一种。方法名称可使用首字母代替。
    weights)    #给不同的x值赋予权重，长度和x相同。默认权重相同。
例：a=density(x,bw=0.5, kernel='c');
3、lines(a)    #添加核密度曲线

绘制多边形
    curve(x*1,0,1,lty=0);
    polygon(c(0.2,0.8,0.7,0.1),c(0.3,0.5,0.8,0.7));
    polygon(c(0.2,0.8,0.7,0.1),c(0.3,0.5,0.8,0.7),col='blue');    #颜色填充
    polygon(c(0.2,0.8,0.7,0.1),c(0.3,0.5,0.8,0.7),density=10);    #斜线填充

切割作图区域
    layout
    par(mfcol, mfrow)
    split.screen

制作一个太极图
    par(pty='s');
    ##########################################
    #中间线
    zjx=function() {
        pz=NULL;
        temp=-1*sqrt(1-seq(0,-1,length=100)^2)-1;
        pz=rbind(pz,cbind(seq(0,-1,length=100),temp));
        temp=sqrt(1-seq(-1,0,length=100)^2)-1;
        pz=rbind(pz,cbind(seq(-1,0,length=100),temp));
        temp=-1*sqrt(1-seq(0,1,length=100)^2)+1;
        pz=rbind(pz,cbind(seq(0,1,length=100),temp));
        temp=sqrt(1-seq(1,0,length=100)^2)+1;
        pz=rbind(pz,cbind(seq(1,0,length=100),temp));
        return(pz);
    }
    #######################################
    #画布
    plot(rnorm,lty=0,bty='n',xaxt='n',yaxt='n',
         xlab='',ylab='',xlim=c(-2,2),ylim=c(-2,2));
  
    #左图
    pz=NULL;
    temp=sqrt(4-seq(0,-2,length=100)^2);
    pz=rbind(pz,cbind(seq(0,-2,length=100),temp));
    temp=-1*sqrt(4-seq(-2,0,length=100)^2);
    pz=rbind(pz,cbind(seq(-2,0,length=100),temp));
    temp=-1*sqrt(1-seq(0,-1,length=100)^2)-1;
    pz=rbind(pz,zjx());
  
    polygon(pz,col='black');
    points(0,1,pch=21,bg='white',cex=10);
  
    #右图
    pz=NULL;
    temp=sqrt(4-seq(0,2,length=100)^2);
    pz=rbind(pz,cbind(seq(0,2,length=100),temp));
    temp=-1*sqrt(4-seq(2,0,length=100)^2);
    pz=rbind(pz,cbind(seq(2,0,length=100),temp));
    pz=rbind(pz,zjx());
  
    polygon(pz,col='white');
    points(0,-1,pch=21,bg='black',cex=10);

R语言矩阵运算

主要包括以下内容：
创建矩阵向量；矩阵加减，乘积；矩阵的逆；行列式的值；特征值与特征向量；QR分解；奇异值分解；广义逆；backsolve与fowardsolve函数；取矩阵的上下三角元素；向量化算子等.

1   创建一个向量
在R中可以用函数c()来创建一个向量，例如：
> x=c(1,2,3,4)
> x
[1] 1 2 3 4 
2   创建一个矩阵
在R中可以用函数matrix()来创建一个矩阵，应用该函数时需要输入必要的参数值。
> args(matrix)
function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL) 
data项为必要的矩阵元素，nrow为行数，ncol为列数，注意nrow与ncol的乘积应为矩阵元素个数，byrow项控制排列元素时是否按行进行，dimnames给定行和列的名称。例如：
> matrix(1:12,nrow=3,ncol=4)
    [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> matrix(1:12,nrow=4,ncol=3)
    [,1] [,2] [,3]
[1,]   1   5   9
[2,]   2   6   10
[3,]   3   7   11
[4,]   4   8   12
> matrix(1:12,nrow=4,ncol=3,byrow=T)
    [,1] [,2] [,3]
[1,]   1   2   3
[2,]   4   5   6
[3,]   7   8   9
[4,]   10   11   12 
> rowname
[1] "r1" "r2" "r3"
> colname=c("c1","c2","c3","c4")
> colname
[1] "c1" "c2" "c3" "c4"
> matrix(1:12,nrow=3,ncol=4,dimnames=list(rowname,colname))
  c1 c2 c3 c4
r1 1 4 7 10
r2 2 5 8 11
3   矩阵转置
A为m×n矩阵，求A'在R中可用函数t()，例如：
> A=matrix(1:12,nrow=3,ncol=4)
> A
   [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> t(A)
   [,1] [,2] [,3]
[1,]   1   2   3
[2,]   4   5   6
[3,]   7   8   9
[4,]   10   11   12若将函数t()作用于一个向量x，则R默认x为列向量，返回结果为一个行向量，例如：
> x
[1] 1 2 3 4 5 6 7 8 9 10
> t(x)
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]   1   2   3   4   5   6   7   8   9   10
> class(x)
[1] "integer"
> class(t(x))
[1] "matrix"
若想得到一个列向量，可用t(t(x))，例如：
> x
[1] 1 2 3 4 5 6 7 8 9 10
> t(t(x))
    [,1]
[1,]   1
[2,]   2
[3,]   3
[4,]   4
[5,]   5
[6,]   6
[7,]   7
[8,]   8
[9,]   9
[10,]  10
> y=t(t(x))
> t(t(y))
    [,1]
[1,]   1
[2,]   2
[3,]   3
[4,]   4
[5,]   5
[6,]   6
[7,]   7
[8,]   8
[9,]   9
[10,]   104   矩阵相加减
在R中对同行同列矩阵相加减，可用符号：“＋”、“－”，例如：
> A=B=matrix(1:12,nrow=3,ncol=4)
> A+B
    [,1] [,2] [,3] [,4]
[1,]   2   8   14   20
[2,]   4   10   16   22
[3,]   6   12   18   24
> A-B
   [,1] [,2] [,3] [,4]
[1,]   0   0   0   0
[2,]   0   0   0   0
[3,]   0   0   0   05   数与矩阵相乘
A为m×n矩阵，c>0，在R中求cA可用符号：“*”，例如：
> c=2
> c*A
    [,1] [,2] [,3] [,4]
[1,]   2   8   14   20
[2,]   4   10  16   22
[3,]   6   12  18   24
6   矩阵相乘
A为m×n矩阵，B为n×k矩阵，在R中求AB可用符号：“％*％”，例如：
> A=matrix(1:12,nrow=3,ncol=4)
> B=matrix(1:12,nrow=4,ncol=3)
> A%*%B
    [,1] [,2] [,3]
[1,]   70  158 246
[2,]   80  184 288
[3,]   90  210 330
若A为n×m矩阵，要得到A'B，可用函数crossprod()，该函数计算结果与t(A)%*%B相同，但是效率更高。例如：
> A=matrix(1:12,nrow=4,ncol=3)
> B=matrix(1:12,nrow=4,ncol=3)
> t(A)%*%B
    [,1] [,2] [,3]
[1,]  30   70 110
[2,]  70  174 278
[3,] 110  278 446
> crossprod(A,B)
    [,1] [,2] [,3]
[1,]  30  70 110
[2,]  70 174 278
[3,] 110 278 446矩阵Hadamard积：若A={aij}m×n, B={bij}m×n, 则矩阵的Hadamard积定义为：
A⊙B={aij bij }m×n,R中Hadamard积可以直接运用运算符“*”例如：
> A=matrix(1:16,4,4)
> A
    [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
> B=A
> A*B
    [,1] [,2] [,3] [,4]
[1,]   1   25   81 169
[2,]   4   36 100 196
[3,]   9   49 121 225
[4,]   16   64 144 256R中这两个运算符的区别区加以注意。7   矩阵对角元素相关运算
例如要取一个方阵的对角元素，
> A=matrix(1:16,nrow=4,ncol=4)
> A
    [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
> diag(A)
[1] 1 6 11 16对一个向量应用diag()函数将产生以这个向量为对角元素的对角矩阵，例如：
> diag(diag(A))
    [,1] [,2] [,3] [,4]
[1,]   1   0   0   0
[2,]   0   6   0   0
[3,]   0   0   11   0
[4,]   0   0   0   16
对一个正整数z应用diag()函数将产生以z维单位矩阵，例如：
> diag(3)
    [,1] [,2] [,3]
[1,]   1   0   0
[2,]   0   1   0
[3,]   0   0   1
8   矩阵求逆
矩阵求逆可用函数solve()，应用solve(a, b)运算结果是解线性方程组ax = b，若b缺省，则系统默认为单位矩阵，因此可用其进行矩阵求逆，例如：
> a=matrix(rnorm(16),4,4)
> a
            [,1]     [,2]     [,3]     [,4]
[1,] 1.6986019   0.5239738 0.2332094 0.3174184
[2,] -0.2010667 1.0913013 -1.2093734   0.8096514
[3,] -0.1797628 -0.7573283 0.2864535 1.3679963
[4,] -0.2217916 -0.3754700 0.1696771 -1.2424030
> solve(a)
              [,1]     [,2]     [,3]     [,4]
[1,] 0.9096360 0.54057479 0.7234861 1.3813059
[2,] -0.6464172 -0.91849017 -1.7546836 -2.6957775
[3,] -0.7841661 -1.78780083 -1.5795262 -3.1046207
[4,] -0.0741260 -0.06308603 0.1854137 -0.6607851
> solve (a) %*%a
                [,1]       [,2]           [,3]       [,4]
[1,] 1.000000e+00 2.748453e-17 -2.787755e-17 -8.023096e-17
[2,] 1.626303e-19 1.000000e+00 -4.960225e-18 6.977925e-16
[3,] 2.135878e-17 -4.629543e-17 1.000000e+00 6.201636e-17
[4,] 1.866183e-17 1.563962e-17 1.183813e-17 1.000000e+009   矩阵的特征值与特征向量
矩阵A的谱分解为A=UΛU',其中Λ是由A的特征值组成的对角矩阵，U的列为A的特征值对应的特征向量，在R中可以用函数eigen()函数得到U和Λ，
> args(eigen)
function (x, symmetric, only.values = FALSE, EISPACK = FALSE)其中：x为矩阵，symmetric项指定矩阵x是否为对称矩阵，若不指定，系统将自动检测x是否为对称矩阵。例如：
> A=diag(4)+1
> A
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> A.eigen=eigen(A,symmetric=T)
> A.eigen
values[1]5111vectors
        [,1]     [,2]       [,3]     [,4]
[1,] 0.5 0.8660254 0.000000e+00 0.0000000
[2,] 0.5 -0.2886751 -6.408849e-17 0.8164966
[3,] 0.5 -0.2886751 -7.071068e-01 -0.4082483
[4,] 0.5 -0.2886751 7.071068e-01 -0.4082483

> A.eigenvectorsvalues)%*%t(A.eigenvectors)[,1][,2][,3][,4][1,]2111[2,]1211[3,]1121[4,]1112>t(A.eigenvectors)%*%A.eigen$vectors
            [,1]       [,2]         [,3]         [,4]
[1,] 1.000000e+00 4.377466e-17 1.626303e-17 -5.095750e-18
[2,] 4.377466e-17 1.000000e+00 -1.694066e-18 6.349359e-18
[3,] 1.626303e-17 -1.694066e-18 1.000000e+00 -1.088268e-16
[4,] -5.095750e-18 6.349359e-18 -1.088268e-16 1.000000e+00
10   矩阵的Choleskey分解
  对于正定矩阵A，可对其进行Choleskey分解，即：A=P'P，其中P为上三角矩阵，在R中可以用函数chol()进行Choleskey分解，例如：
> A
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> chol(A)
        [,1]     [,2]     [,3]     [,4]
[1,] 1.414214 0.7071068 0.7071068 0.7071068
[2,] 0.000000 1.2247449 0.4082483 0.4082483
[3,] 0.000000 0.0000000 1.1547005 0.2886751
[4,] 0.000000 0.0000000 0.0000000 1.1180340
> t(chol(A))%*%chol(A)
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> crossprod(chol(A),chol(A))
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2若矩阵为对称正定矩阵，可以利用Choleskey分解求行列式的值，如：
> prod(diag(chol(A))^2)
[1] 5
> det(A)
[1] 5若矩阵为对称正定矩阵，可以利用Choleskey分解求矩阵的逆，这时用函数chol2inv()，这种用法更有效。如：
> chol2inv(chol(A))
      [,1] [,2] [,3] [,4]
[1,] 0.8 -0.2 -0.2 -0.2
[2,] -0.2 0.8 -0.2 -0.2
[3,] -0.2 -0.2 0.8 -0.2
[4,] -0.2 -0.2 -0.2 0.8
> solve(A)
  [,1] [,2] [,3] [,4]
[1,] 0.8 -0.2 -0.2 -0.2
[2,] -0.2 0.8 -0.2 -0.2
[3,] -0.2 -0.2 0.8 -0.2
[4,] -0.2 -0.2 -0.2 0.811   矩阵奇异值分解
  A为m×n矩阵，rank(A)= r, 可以分解为：A=UDV',其中U'U=V'V=I。在R中可以用函数scd()进行奇异值分解，例如：
> A=matrix(1:18,3,6)
> A
  [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   1   4   7   10   13   16
[2,]   2   5   8   11   14   17
[3,]   3   6   9   12   15   18
> svd(A)
d[1]4.589453e+011.640705e+003.627301e−16u
          [,1]     [,2]     [,3]
[1,] -0.5290354 0.74394551 0.4082483
[2,] -0.5760715 0.03840487 -0.8164966
[3,] -0.6231077 -0.66713577 0.4082483
v[,1][,2][,3][1,]−0.07736219−0.7196003−0.18918124[2,]−0.19033085−0.50893250.42405898[3,]−0.30329950−0.2982646−0.45330031[4,]−0.41626816−0.0875968−0.01637004[5,]−0.529236820.12307110.64231130[6,]−0.642205480.3337389−0.40751869>A.svd=svd(A)>A.svdu%*%diag(A.svdd)v)
  [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   1   4   7   10   13   16
[2,]   2   5   8   11   14   17
[3,]   3   6   9   12   15   18
> t(A.svdu)u
            [,1]       [,2]       [,3]
[1,] 1.000000e+00 -1.169312e-16 -3.016793e-17
[2,] -1.169312e-16 1.000000e+00 -3.678156e-17
[3,] -3.016793e-17 -3.678156e-17 1.000000e+00
> t(A.svdv)v
        [,1]       [,2]       [,3]
[1,] 1.000000e+00 8.248068e-17 -3.903128e-18
[2,] 8.248068e-17 1.000000e+00 -2.103352e-17
[3,] -3.903128e-18 -2.103352e-17 1.000000e+0012   矩阵QR分解
A为m×n矩阵可以进行QR分解，A=QR，其中：Q'Q＝I，在R中可以用函数qr()进行QR分解，例如：
> A=matrix(1:16,4,4)
> qr(A)
qr[,1][,2][,3][,4][1,]−5.4772256−12.7801930−2.008316e+01−2.738613e+01[2,]0.3651484−3.2659863−6.531973e+00−9.797959e+00[3,]0.5477226−0.37816962.641083e−152.056562e−15[4,]0.7302967−0.91247448.583032e−01−2.111449e−16rank
[1] 2

qraux[1]1.182574e+001.156135e+001.513143e+002.111449e−16pivot
[1] 1 2 3 4

attr(,"class")
[1] "qr"rank项返回矩阵的秩，qr项包含了矩阵Q和R的信息，要得到矩阵Q和R，可以用函数qr.Q()和qr.R()作用qr()的返回结果，例如：
> qr.R(qr(A))
      [,1]     [,2]       [,3]       [,4]
[1,] -5.477226 -12.780193 -2.008316e+01 -2.738613e+01
[2,] 0.000000 -3.265986 -6.531973e+00 -9.797959e+00
[3,] 0.000000   0.000000 2.641083e-15 2.056562e-15
[4,] 0.000000   0.000000 0.000000e+00 -2.111449e-16
> qr.Q(qr(A))
      [,1]       [,2]     [,3]     [,4]
[1,] -0.1825742 -8.164966e-01 -0.4000874 -0.37407225
[2,] -0.3651484 -4.082483e-01 0.2546329 0.79697056
[3,] -0.5477226 -8.131516e-19 0.6909965 -0.47172438
[4,] -0.7302967 4.082483e-01 -0.5455419 0.04882607
> qr.Q(qr(A))%*%qr.R(qr(A))
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
> t(qr.Q(qr(A)))%*%qr.Q(qr(A))
        [,1]       [,2]       [,3]       [,4]
[1,] 1.000000e+00 -1.457168e-16 -6.760001e-17 -7.659550e-17
[2,] -1.457168e-16 1.000000e+00 -4.269046e-17 7.011739e-17
[3,] -6.760001e-17 -4.269046e-17 1.000000e+00 -1.596437e-16
[4,] -7.659550e-17 7.011739e-17 -1.596437e-16 1.000000e+00
> qr.X(qr(A))
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   1613   矩阵广义逆(Moore-Penrose)
  n×m矩阵A+称为m×n矩阵A的Moore-Penrose逆，如果它满足下列条件：
①   A A+A=A；②A+A A+= A+；③(A A+)H=A A+；④(A+A)H= A+A
在R的MASS包中的函数ginv()可计算矩阵A的Moore-Penrose逆，例如：
library(“MASS”)
> A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
> ginv(A)
    [,1]   [,2] [,3]   [,4]
[1,] -0.285 -0.1075 0.07 0.2475
[2,] -0.145 -0.0525 0.04 0.1325
[3,] -0.005 0.0025 0.01 0.0175
[4,] 0.135 0.0575 -0.02 -0.0975
验证性质1：
> A%*%ginv(A)%*%A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
验证性质2：
> ginv(A)%*%A%*%ginv(A)
    [,1]   [,2] [,3]   [,4]
[1,] -0.285 -0.1075 0.07 0.2475
[2,] -0.145 -0.0525 0.04 0.1325
[3,] -0.005 0.0025 0.01 0.0175
[4,] 0.135 0.0575 -0.02 -0.0975
验证性质3:
> t(A%*%ginv(A))
  [,1] [,2] [,3] [,4]
[1,] 0.7 0.4 0.1 -0.2
[2,] 0.4 0.3 0.2 0.1
[3,] 0.1 0.2 0.3 0.4
[4,] -0.2 0.1 0.4 0.7
> A%*%ginv(A)
  [,1] [,2] [,3] [,4]
[1,] 0.7 0.4 0.1 -0.2
[2,] 0.4 0.3 0.2 0.1
[3,] 0.1 0.2 0.3 0.4
[4,] -0.2 0.1 0.4 0.7
验证性质4:
> t(ginv(A)%*%A)
  [,1] [,2] [,3] [,4]
[1,] 0.7 0.4 0.1 -0.2
[2,] 0.4 0.3 0.2 0.1
[3,] 0.1 0.2 0.3 0.4
[4,] -0.2 0.1 0.4 0.7
> ginv(A)%*%A
  [,1] [,2] [,3] [,4]
[1,] 0.7 0.4 0.1 -0.2
[2,] 0.4 0.3 0.2 0.1
[3,] 0.1 0.2 0.3 0.4
[4,] -0.2 0.1 0.4 0.714   矩阵Kronecker积
  n×m矩阵A与h×k矩阵B的kronecker积为一个nh×mk维矩阵，
在R中kronecker积可以用函数kronecker()来计算，例如：
> A=matrix(1:4,2,2)
> B=matrix(rep(1,4),2,2)
> A
  [,1] [,2]
[1,]   1   3
[2,]   2   4
> B
  [,1] [,2]
[1,]   1   1
[2,]   1   1
> kronecker(A,B)
  [,1] [,2] [,3] [,4]
[1,]   1   1   3   3
[2,]   1   1   3   3
[3,]   2   2   4   4
[4,]   2   2   4   415   矩阵的维数
  在R中很容易得到一个矩阵的维数，函数dim()将返回一个矩阵的维数，nrow()返回行数，ncol()返回列数，例如：
  > A=matrix(1:12,3,4)
> A
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> nrow(A)
[1] 3
> ncol(A)
[1] 416   矩阵的行和、列和、行平均与列平均
  在R中很容易求得一个矩阵的各行的和、平均数与列的和、平均数，例如：
  > A
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> rowSums(A)
[1] 22 26 30
> rowMeans(A)
[1] 5.5 6.5 7.5
> colSums(A)
[1] 6 15 24 33
> colMeans(A)
[1] 2 5 8 11上述关于矩阵行和列的操作，还可以使用apply()函数实现。
> args(apply)
function (X, MARGIN, FUN, ...)
其中：x为矩阵，MARGIN用来指定是对行运算还是对列运算，MARGIN＝1表示对行运算，MARGIN＝2表示对列运算，FUN用来指定运算函数, ...用来给定FUN中需要的其它的参数，例如：
> apply(A,1,sum)
[1] 22 26 30
> apply(A,1,mean)
[1] 5.5 6.5 7.5
> apply(A,2,sum)
[1] 6 15 24 33
> apply(A,2,mean)
[1] 2 5 8 11apply()函数功能强大，我们可以对矩阵的行或者列进行其它运算，例如：
计算每一列的方差
> A=matrix(rnorm(100),20,5)
> apply(A,2,var)
[1] 0.4641787 1.4331070 0.3186012 1.3042711 0.5238485
> apply(A,2,function(x,a)x*a,a=2)
  [,1] [,2] [,3] [,4]
[1,]   2   8   14   20
[2,]   4   10   16   22
[3,]   6   12   18   24注意：apply(A,2,function(x,a)x*a,a=2)与A*2效果相同，此处旨在说明如何应用alpply函数。17   矩阵X'X的逆
  在统计计算中，我们常常需要计算这样矩阵的逆，如OLS估计中求系数矩阵。R中的包“strucchange”提供了有效的计算方法。
  > args(solveCrossprod)
function (X, method = c("qr", "chol", "solve"))其中：method指定求逆方法，选用“qr”效率最高，选用“chol”精度最高，选用“slove”与slove(crossprod(x,x))效果相同，例如：
> A=matrix(rnorm(16),4,4)
> solveCrossprod(A,method="qr")
      [,1]     [,2]     [,3]     [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
> solveCrossprod(A,method="chol")
      [,1]     [,2]     [,3]     [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
> solveCrossprod(A,method="solve")
      [,1]     [,2]     [,3]     [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
> solve(crossprod(A,A))
      [,1]     [,2]     [,3]     [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.344762718   取矩阵的上、下三角部分
  在R中，我们可以很方便的取到一个矩阵的上、下三角部分的元素，函数lower.tri()和函数upper.tri()提供了有效的方法。
  > args(lower.tri)
function (x, diag = FALSE)函数将返回一个逻辑值矩阵，其中下三角部分为真，上三角部分为假，选项diag为真时包含对角元素，为假时不包含对角元素。upper.tri()的效果与之孑然相反。例如：
> A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
> lower.tri(A)
    [,1] [,2] [,3] [,4]
[1,] FALSE FALSE FALSE FALSE
[2,] TRUE FALSE FALSE FALSE
[3,] TRUE TRUE FALSE FALSE
[4,] TRUE TRUE TRUE FALSE
> lower.tri(A,diag=T)
  [,1] [,2] [,3] [,4]
[1,] TRUE FALSE FALSE FALSE
[2,] TRUE TRUE FALSE FALSE
[3,] TRUE TRUE TRUE FALSE
[4,] TRUE TRUE TRUE TRUE
> upper.tri(A)
    [,1] [,2] [,3] [,4]
[1,] FALSE TRUE TRUE TRUE
[2,] FALSE FALSE TRUE TRUE
[3,] FALSE FALSE FALSE TRUE
[4,] FALSE FALSE FALSE FALSE
> upper.tri(A,diag=T)
    [,1] [,2] [,3] [,4]
[1,] TRUE TRUE TRUE TRUE
[2,] FALSE TRUE TRUE TRUE
[3,] FALSE FALSE TRUE TRUE
[4,] FALSE FALSE FALSE TRUE
> A[lower.tri(A)]=0
> A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   0   6   10   14
[3,]   0   0   11   15
[4,]   0   0   0   16
> A[upper.tri(A)]=0
> A
  [,1] [,2] [,3] [,4]
[1,]   1   0   0   0
[2,]   2   6   0   0
[3,]   3   7   11   0
[4,]   4   8   12   1619   backsolve&fowardsolve函数
这两个函数用于解特殊线性方程组，其特殊之处在于系数矩阵为上或下三角。
> args(backsolve)
function (r, x, k = ncol(r), upper.tri = TRUE, transpose = FALSE)
> args(forwardsolve)
function (l, x, k = ncol(l), upper.tri = FALSE, transpose = FALSE)其中：r或者l为n×n维三角矩阵，x为n×1维向量，对给定不同的upper.tri和transpose的值，方程的形式不同
对于函数backsolve()而言，
例如：
  > A=matrix(1:9,3,3)
> A
  [,1] [,2] [,3]
[1,]   1   4   7
[2,]   2   5   8
[3,]   3   6   9
> x=c(1,2,3)
> x
[1] 1 2 3
> B=A
> B[upper.tri(B)]=0
> B
  [,1] [,2] [,3]
[1,]   1   0   0
[2,]   2   5   0
[3,]   3   6   9
> C=A
> C[lower.tri(C)]=0
> C
  [,1] [,2] [,3]
[1,]   1   4   7
[2,]   0   5   8
[3,]   0   0   9
> backsolve(A,x,upper.tri=T,transpose=T)
[1] 1.00000000 -0.40000000 -0.08888889
> solve(t(C),x)
[1] 1.00000000 -0.40000000 -0.08888889
> backsolve(A,x,upper.tri=T,transpose=F)
[1] -0.8000000 -0.1333333 0.3333333
> solve(C,x)
[1] -0.8000000 -0.1333333 0.3333333
> backsolve(A,x,upper.tri=F,transpose=T)
[1] 1.111307e-17 2.220446e-17 3.333333e-01
> solve(t(B),x)
[1] 1.110223e-17 2.220446e-17 3.333333e-01
> backsolve(A,x,upper.tri=F,transpose=F)
[1] 1 0 0
> solve(B,x)
[1] 1.000000e+00 -1.540744e-33 -1.850372e-17
对于函数forwardsolve()而言，
例如：
  > A
      [,1] [,2] [,3]
[1,]   1   4   7
[2,]   2   5   8
[3,]   3   6   9
> B
  [,1] [,2] [,3]
[1,]   1   0   0
[2,]   2   5   0
[3,]   3   6   9
> C
  [,1] [,2] [,3]
[1,]   1   4   7
[2,]   0   5   8
[3,]   0   0   9
> x
[1] 1 2 3
> forwardsolve(A,x,upper.tri=T,transpose=T)
[1] 1.00000000 -0.40000000 -0.08888889
> solve(t(C),x)
[1] 1.00000000 -0.40000000 -0.08888889
> forwardsolve(A,x,upper.tri=T,transpose=F)
[1] -0.8000000 -0.1333333 0.3333333
> solve(C,x)
[1] -0.8000000 -0.1333333 0.3333333
> forwardsolve(A,x,upper.tri=F,transpose=T)
[1] 1.111307e-17 2.220446e-17 3.333333e-01
> solve(t(B),x)
[1] 1.110223e-17 2.220446e-17 3.333333e-01
> forwardsolve(A,x,upper.tri=F,transpose=F)
[1] 1 0 0
> solve(B,x)
[1] 1.000000e+00 -1.540744e-33 -1.850372e-1720   row()与col()函数
在R中定义了的这两个函数用于取矩阵元素的行或列下标矩阵，例如矩阵A={aij}m×n，
row()函数将返回一个与矩阵A有相同维数的矩阵，该矩阵的第i行第j列元素为i，函数col()类似。例如：
> x=matrix(1:12,3,4)
> row(x)
  [,1] [,2] [,3] [,4]
[1,]   1   1   1   1
[2,]   2   2   2   2
[3,]   3   3   3   3
> col(x)
  [,1] [,2] [,3] [,4]
[1,]   1   2   3   4
[2,]   1   2   3   4
[3,]   1   2   3   4这两个函数同样可以用于取一个矩阵的上下三角矩阵，例如：
> x
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> x[row(x)<col(x)]=0
> x
  [,1] [,2] [,3] [,4]
[1,]   1   0   0   0
[2,]   2   5   0   0
[3,]   3   6   9   0
> x=matrix(1:12,3,4)
> x[row(x)>col(x)]=0
> x
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   0   5   8   11
[3,]   0   0   9   1221   行列式的值
在R中，函数det(x)将计算方阵x的行列式的值，例如：
> x=matrix(rnorm(16),4,4)
> x
      [,1]     [,2]     [,3]     [,4]
[1,] -1.0736375 0.2809563 -1.5796854 0.51810378
[2,] -1.6229898 -0.4175977 1.2038194 -0.06394986
[3,] -0.3989073 -0.8368334 -0.6374909 -0.23657088
[4,] 1.9413061 0.8338065 -1.5877162 -1.30568465
> det(x)
[1] 5.717667
22向量化算子
在R中可以很容易的实现向量化算子，例如：

vec<-function (x){
          t(t(as.vector(x)))
}
vech<-function (x){
          t(x[lower.tri(x,diag=T)])
}
> x=matrix(1:12,3,4)
> x
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> vec(x)
    [,1]
[1,]   1
[2,]   2
[3,]   3
[4,]   4
[5,]   5
[6,]   6
[7,]   7
[8,]   8
[9,]   9
[10,]   10
[11,]   11
[12,]   12
> vech(x)
  [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   1   2   3   5   6   923   时间序列的滞后值
  在时间序列分析中，我们常常要用到一个序列的滞后序列，R中的包“fMultivar”中的函数tslag()提供了这个功能。
  > args(tslag)
function (x, k = 1, trim = FALSE)其中：x为一个向量，k指定滞后阶数，可以是一个自然数列，若trim为假，则返回序列与原序列长度相同，但含有NA值；若trim项为真，则返回序列中不含有NA值，例如：
> x=1:20
> tslag(x,1:4,trim=F)
    [,1] [,2] [,3] [,4]
[1,]   NA   NA   NA   NA
[2,]   1   NA   NA   NA
[3,]   2   1   NA   NA
[4,]   3   2   1   NA
[5,]   4   3   2   1
[6,]   5   4   3   2
[7,]   6   5   4   3
[8,]   7   6   5   4
[9,]   8   7   6   5
[10,]   9   8   7   6
[11,]   10   9   8   7
[12,]   11   10   9   8
[13,]   12   11   10   9
[14,]   13   12   11   10
[15,]   14   13   12   11
[16,]   15   14   13   12
[17,]   16   15   14   13
[18,]   17   16   15   14
[19,]   18   17   16   15
[20,]   19   18   17   16
> tslag(x,1:4,trim=T)
    [,1] [,2] [,3] [,4]
[1,]   4   3   2   1
[2,]   5   4   3   2
[3,]   6   5   4   3
[4,]   7   6   5   4
[5,]   8   7   6   5
[6,]   9   8   7   6
[7,]   10   9   8   7
[8,]   11   10   9   8
[9,]   12   11   10   9
[10,]   13   12   11   10
[11,]   14   13   12   11
[12,]   15   14   13   12
[13,]   16   15   14   13
[14,]   17   16   15   14
[15,]   18   17   16   15
[16,]   19   18   17   16