Blabla
##通过在命令、变量、逗号等附近添加空格来提高R的可读性.
##可使用control+tab键来在控制台console和图形设备graphic device间切换. >help.start()：得到HTML 格式的帮助(等价于:帮助->html帮助) >?xxx
or
>help(x)： 得到任何特定名字的函数的帮助
##对于关键字和运算符，与函数的帮助类似，但是需要加上引号， 如：> ?'+' 或>?”+” #等价于 help('+') help(“+”)
>??xx :在help.start()启动的浏览页上,”Search Engine &Keywords”
>objects() or >ls():用来显示当前保存在R 环境中的对象名字
>ls(pat=“x”):只显示在名称中带有某个指定字符的对象 >ls(pat=“^x”):只显示在名称中以某个指定字符开头的对象 >ls()不能列出名称以点号开头的向量 可使用ls(all = T) >history():看到之前保存的数据和命令。 >rm(x,y): 删除对象
>rm(list=ls(all=TRUE)) : 删除内存中所有对象,all=TRUE可以省略 >dir.create(\"c:/Users/Mr.Young/Desktop/R\") 创建新的工作目录 >getwd()和>setwd():获取/设置工作空间目录 ##临时修改，只针对当前文件
>setwd(file=\"c:/users/Mr.Young/Documents\")
file可以省略，且该操作等价于“文件”菜单“change dir” >getwd()
[1] \"c:/users/Mr.Young/Documents\" （目录的分隔符用“/”(slash)或“\\\\”） 永久修改，R右键->属性->快捷方式->起始位置 Windows下默认的分隔符为 \\ (backslash),注意更改. >example():显示函数的例子（不加引号） >list.files():查看当前目录下的文件 >save(x,y,file=”d:/xy.Rdata”); #保存变量 >load(”d:/xy.Rdata”); #载入变量
1
%%:求余数 %/%:整除
逻辑与: &(拉丁文为et)（前后都要判断）or && ( “&&”前为F就结束，结果为F) 逻辑非:
| (前后都要判断) or || （“||”前为T就结束，结果为T）
##同一行中输入多条命令，用“;”隔开，否则另起一行。
注释：一行中，从“#”或“##(可有多个#)”开始到句子收尾之间的语句 >log():默认e为底
logbase(x) or >log(x,base):以base为底数的对数 >e必须用 exp(1) ##数据的浏览 >head(x):浏览前6个 >tail(x):浏览后6个 ##数据的编辑
>data.entry():打开数据编辑器
>edit():打开R editor，但是不能直接修改，直接修改必须使前后的数据集同名，或者对原来的数据集修改后赋给一个新的变量，而原来的数据集不变 >fix():打开R editor，可直接修改，fix(a)等价于a=edit(a) A = edit(data.frame()) %in%
1:10 %in% c(1,3,5,9)
[1] TRUE FALSE TRUE FALSE TRUE FALSE FALSE FALSE TRUE FALSE all.equal identical >ceiling():上取整 >floor():下取整 >round(x, digits = 0) >install.packages(“”)
>detach(“package : xxxx”) 删除包
>library()：加载包，用“”，‘’，也可都不用 >detach(package:xxx)关闭包
2
>apropos(“lm”)找到包含“lm”字符的函数或数据 Mathematical Annotation in R
http://127.0.0.1:19703/library/grDevices/html/plotmath.html library(psych) describe(cars)
R中利用sd()或describe()计算出的标准差与公式计算出的不同，因为在R中默认的是将标准差用于推理性统计，分母上为自由度N-1而非N >citation(\"pkgname\"):给出pskname包的相关信息，便于引用 To cite R and R packages in publications
>(A=sqrt(25)):通过为表达式加括号可以直接得到结果. >q():退出R q(save = ”no”):R不做保存直接退出
出现的图形可以用Ctrl+W和Ctrl+C来复制和粘贴（前者像素高） 在运行时按Esc会终止运行
复数可以直接表示为a+bi，注：1+i必须写为1+1i
> X = c(TRUE,FALSE) > X
[1] TRUE FALSE > X*1 [1] 1 0
3
R的数据结构
>mode():对象的类型(数值型numeric，复数型complex，逻辑型logical，字符型character) ##实际上，R 已经有自己独立的函数typeof()，仍然保留mode概念主要是为了和S 兼容。数值型模式numeric实际上是两种独立模式的混合模式，即整数型（integer）和双精度型（double） >class():表示类别
##typeof()处理对象内元素的类型，而class()处理对象本身的类
>x [,1] [,2] [1,] 1 4 [2,] 2 5 [3,] 3 6 > typeof(x) [1] \"integer\" > class(x) [1] \"matrix\" ##通过class还可以更改对象的类: > x = 1:6 > class(x) [1] \"integer\" > class(x) = \"matrix\" 错误于除非维度的长度为二(目前是 0)，否则不能设为矩阵类别
##输入字符型的值时需加上双引号“”，如果恰好需要引用双引号的话，可以在双引号引用前加上反斜杠“\\”, 也可以直接用单引号‘’来界定，此时不需要加“\\”来引用双引号 > a=\"r\\\"data\" > a [1] \"r\\\"data\" >length(xx):对象的长度
缺损值：NA( not available) NaN（Not a Number） e.g.> Inf – Inf >0/0 ##空的字符串向量将会被显示character(0)和空的数值向量显示为numeric(0)
4
> class(x) = \"logical\" > x [1] TRUE TRUE TRUE TRUE TRUE TRUE > b='r\"data' > b [1] \"r\\\"data\" ##向量由一系列类型相同的有序元素构成；矩阵是数组(array)的一个特例：维数为2的数组；而数组又是增加了维度(dim)属性的向量。 向量的建立：
数值型： c()：用c函数连接数据（Concatenate）
seq() rep()
> 1:10 > rep(2:5,2) [1] 2 3 4 5 2 3 4 5 [1] 1 2 3 4 5 6 7 8 9 10 > rep(1:3,times=3,each=2) > 10:1 [1] 10 9 8 7 6 5 4 3 2 1 [1] 1 1 2 2 3 3 1 1 2 2 3 3 1 1 2 2 3 3 > 1:10-1 [1] 0 1 2 3 4 5 6 7 8 9 > seq(1,5,by=0.5) [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 > seq(1,10,length=3) [1] 1.0 5.5 10.0 注: scan() 从键盘逐个输入 > f=scan() # 1: 1 2 3 4 5 6: Read 5 items > > f [1] 1 2 3 4 5 length处正确应该为length.out 实际leng也行(模糊识别) R的参数为模糊匹配
A=vector(length=6)：事先定义向量长度，再对A[1]、A[2]等赋值，循环时有用 逻辑型：
> x=c(10,9,8,7,6,5) > y=x>7 > y [1] TRUE TRUE TRUE FALSE FALSE FALSE >sort(x) 返回一个和x长度一样但元素以升序排列的向量 >sort(x,decreasing=TRUE) 降序
>order(x):返回秩（升序）decreasing=TRUE则降序 >prod(x)：（product）得到向量x各元素的乘积
>range(x) ：得到的是一个长度为2的向量，即c(min(x), max(x))
矩阵的建立：
① x=c(10,9,8,7,6,5);dim(x)=c(2,3)
5
② matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL) 通过nrow(nr),ncol(nc)控制为维度，dimnames包含了可选的、以字符型向量表示的行名和列名
> x=matrix(1:6,2,3) > x
[,1] [,2] [,3] [1,] 1 3 5 [2,] 2 4 6 > x[1,3] [1] 5 > x[2,] [1] 2 4 6 > x[,2] [1] 3 4 > x[,2:3] [,1] [,2] [1,] 3 5 [2,] 4 6 > x[3] [1] 3 （##按列找的第三个）
> y=matrix(1:6,2,3,byrow=TRUE) ##按列填充 > y [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 > dim(y) [1] 2 3 > z=diag(3) > z [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] 0 0 1 > z=diag(c(3,4,5)) > z [,1] [,2] [,3] [1,] 3 0 0 [2,] 0 4 0 [3,] 0 0 5 用两种方式改名 r=c(\"r1\ c=c(\"c1\ a=matrix(1:12,nc=3,dimnames=list(r,c)) b=matrix(1:12,nc=3) rownames(b)=r colnames(b)=c >crossprod(X,Y) 和t(X) %*% y 等价,在运算上更为高效 >nrow(A) 和>ncol(A)分别返回矩阵A 的行数和列数 >t(X) 矩阵的转置 transpose
>solve(X) 矩阵的求逆 solve(A,b) AX=b,求X X%*%Y 矩阵相乘 X * Y 两个矩阵对应元素乘积
6
数组的建立:
array(data = NA, dim = length(data), dimnames = NULL) dim1 <- c(\"A1\A2\") dim2 <- c(\"B1\ dim3 <- c(\"C1\ z <- array(1:12, c(2, 3, 2), dimnames = list(dim1, dim2, dim3)) a[-4]
>drop():Delete the dimensions of an array which have only one level. x <- 1:4 z <- x %*% x # scalar (\"inner\") product z drop(z) A=array(1:27,c(3,3,3)) B=array(27:1,c(3,3,3)) A%*%B=? 答案是：A*B的结果中所有元素的相加 即sum(A*B[,,]) list:包含任何类型的对象的列表 > L1=list(1:6,matrix(1:4,2)) > L1 [[1]] [1] 1 2 3 4 5 6 [[2]] [,1] [,2] [1,] 1 3
7
a[-length(a)]:删除元素
> z , , C1 B1 B2 B3 A1 1 3 5 A2 2 4 6 , , C2 B1 B2 B3 A1 7 9 11 A2 8 10 12 > z [,1] [1,] 30 > drop(z) [1] 30 引用： > L1[[1]] [1] 1 2 3 4 5 6 [2,] 2 4
> L2=list(x=1:6,y=matrix(1:4,2)) > L2 $x
[1] 1 2 3 4 5 6 $y
[,1] [,2] [1,] 1 3 [2,] 2 4
删除list 中的元素，使用NULL 表示无效的对象。
lst[[”a”]][”b”]=NULL # or lst $a$b < NULL
> a=list(2,3) > a [[1]] [1] 2 [[2]] [1] 3 > unlist(a) [1] 2 3
data frame（数据框）
可以理解是一个松散的数据集, 是用于存储数据的一种结构。它可以是由不同类型的列（数字、因子、字符等）组成的类矩阵（matrix-like）列表示变量，行表示观测。在同一个数据框中可以存储不同类型（如数值型、字符型）的变量。
引用： > L2$x [1] 1 2 3 4 5 6
8
ID <- c(1, 2, 3, 4) age <- c(25, 34, 28, 52) type<-c(\"Type1\ 3\") status<-c(\"Poor\\ data<-data.frame(ID,age,type,status) data
数据编辑器
X1=c(1,2,3,4,5) X2=c(\"a\ data.frame(X1,X2) Y=data.frame(数字=X1,字母=X2) Y 可以使用下标操作对象，编辑对象中的数据元素。但是R提供的一个可视化的工具能够带来更多的便利，这就是数据编辑器。使用data.entry()函数可以打开数据编辑器 > x = array(6:1,2:3) > data.entry(x)
apply(X,MARGIN,FUN)：对一个对象施加某种运算 MARGIN=1:按行计算；MARGIN=2：按列计算； MARGIN=c(1,2)：按行列计算(三维以上使用)
-------------------------------------------------------------------------------------------------------- A=matrix(1:24, ,6)
apply(A,2,mean) = colMeans(A) colMeans(A)[c(1,3)]= apply(A[,c(1,3)],2,mean) apply(A,1,mean) A=array(1:24,c(4,3,2)) A , , 1
[,1] [,2] [,3] [1,] 1 5 9 [2,] 2 6 10
, , 2 [,1] [,2] [,3]
[1,] 13 17 21 [2,] 14 18 22 [3,] 15 19 23 9 16 20 24 [4,] [3,] 3 7 11 [4,] 4 8 12 apply(A,3,mean) [1] 6.5 18.5 apply(A,c(1,3),sum) apply(A,c(1,2),sum) [,1] [,2] [,1] [,2] [,3] [1,] 15 51 [1,] 14 22 30 [2,] 18 54 [2,] 16 24 32 [3,] 21 57 [3,] 18 26 34 [4,] 24 60 [4,] 20 28 36 mapply (FUN, SIMPLIFY = TRUE,) > p=list(A=1:10,B=10:20) > p $A
[1] 0 1 2 3 4 5 6 7 8 9 10 $B
[1] 10 11 12 13 14 15 16 17 18 19 20 > f=list(A=50:60,B=60:70) > f $A
[1] 50 51 52 53 54 55 56 57 58 59 60 $B
[1] 60 61 62 63 64 65 66 67 68 69 70 > mapply(sum,p,f) A B 660 880 > sum(0:10,50:60) [1] 660
> myadd=function(x,y) x+y
10
apply(A,c(2,3),sum) [,1] [,2] [1,] 10 58 [2,] 26 74 [3,] 42 90 > mapply(myadd,p,f) A B [1,] 50 70 [2,] 52 72 [3,] 54 74 [4,] 56 76 [5,] 58 78 [6,] 60 80 [7,] 62 82 [8,] 64 84 [9,] 66 86 [10,] 68 88 [11,] 70 90
> mapply(myadd,p,f,SIMPLIFY=FALSE) $A
[1] 50 52 54 56 58 60 62 64 66 68 70 $B
[1] 70 72 74 76 78 80 82 84 86 88 90
> mapply(\"+\",p,f,SIMPLIFY=FALSE) $A
[1] 50 52 54 56 58 60 62 64 66 68 70 $B
[1] 70 72 74 76 78 80 82 84 86 88 90 > att=list(A=rep(1,11),B=rep(1,11)) > att $A
[1] 1 1 1 1 1 1 1 1 1 1 1
11
$B
[1] 1 1 1 1 1 1 1 1 1 1 1
> myadd=function(x,y,z) x+y+z
> mapply(myadd,p,f,att,SIMPLIFY=FALSE) $A
[1] 51 53 55 57 59 61 63 65 67 69 71 $B
[1] 71 73 75 77 79 81 83 85 87 89 91
sapply(X, FUN):对x的每一个变量使用FUN，计算的数据是x的全部值，以向量的形式输出结果.
lapply(X, FUN): 对x的每一个变量使用FUN，计算的数据是x的全部值，以列表的形式输出结果.
tapply(X,INDEX,FUN=)：根据index的不同水平对x进行FUN运算,即计算的数据时x的一个子集
lapply返回一个列表 用“l”标示
sapply尽可能将结果简化成向量/矩阵 用“s”标示 tapply创建表格 用“t”标示
R编程
循环结构重复地执行一个或一系列语句，直到某个条件不为真为止。循环结构包括for和while结构。
for(知道终止条件)循环重复地执行一个语句，直到某个变量的值不再包含在序列seq中为止 for(var in seq) statement
e.g. for(i in 1:5) print(1:i)
while(无法知道运行次数)循环重复地执行一个语句，直到条件不为真为止(用i时要先定义哦)
i=1;while(i<=5) {print(1:i) ; i=i+1}
12
R自编函数结构形式如下:
myfunction<-function(arg1,arg2,…) { statements ; return(object) } 注意:>for( i in 1:10 ){ d[i]=i}
错误于d[i] = i : 找不到对象'd' 应该先设定d=0 d=NULL
R中内嵌的分布
d 密度函数PDF 即f(x) probability density function
p（累积）分布函数CDF 即F(x) cumulative distribution function q 分位数函数quantile r 随机数函数random dnorm(x, mean = 0, sd = 1, log = FALSE) ##概率密度函数f(x) dnorm(x, mean = 0, sd = 1, log=TRUE)= log(dnorm(x, mean = 0, sd = 1) > x=1e500 > log(x) [1] Inf > 500*log(10)
[1] 1151.293 (先取对数的重要性! more stable for very large values！) pnorm(q, mean = 0, sd = 1)
qnorm (p, mean = 0, sd = 1) (Given CDF, what is x?), 分位函数 rnorm(x,mean=0,sd=1) 随机模拟或者随机数发生器 e.g.
pbinom是求二项分布的分布函数值。
比如X服从B(100,0.5) 那么P(X<=55)怎么求？用pbinom(55,100,0.5) 那P(X>55)呢，当然你可以用1-pbinom(55,100,0.5) 或者 pbinom(55,100,0.5,lower.tail=TRUE)
##随机序列用于产生符合一定分布规则的数据。有大量的函数用于产生随机序列，这里只列出一些函数的名称：
13
set.seed(1)
##计算机并不能产生真正的随机数，如果你不设种子，计算机会用系统时钟来作为种子，如果你要模拟什么的话，每次的随机数都是不一样的，这样就不方便你研究，如果你事先设置了种子，这样每次的随机数都是一样的，便于重现你的研究，也便于其他人检验你的分析结果。
--------------------------------------------------------------------------------------------------------
runif(5) set.seed(1234) runif(5) set.seed(1234) runif(5)
-------------------------------------------------------------------------------------------------------- 对数正态分布lnorm(μ,𝜎2)
ln(X)服从参数为(μ,σ^2)的正态分布，则X的分布就是对数正态分布
求似然函数
14
写dnorm函数
随机数的产生
Direct methods
Probability Integral Transform
考虑用均匀分布的随机数产生指数分布的随机数，前提是知道指数分布的CDF -------------------------------------------------------------------------------------------------------- u2x=function(u,lambda) { } n=100 u=runif(n) lambda=0.5 x=u2x(u,lambda)
hist(x,breaks=20,freq=FALSE) mean(x)
xax=seq(0,max(x),0.01)
xdens=lambda*exp(-lambda*xax) lines(xax,xdens,col=\"blue\windows()
plot(xax,xdens,type=\"l\
--------------------------------------------------------------------------------------------------------
x=-log(1-u)/lambda return(x)
15
-------------------------------------------------------------------------------------------------------- 用rexp和diy的结果作对比 n=1000 u=runif(n) lambda=0.5
u2x=function(u,lambda) { }
myrexp=function(n,lambda) {
u=runif(n,0,1) x=u2x(u,lambda) x=-log(1-u)/lambda return(x)
return(x) }
par(mfrow=c(1,2))
16
Rrexp=rexp(n,rate=lambda) x=u2x(u,lambda) xax=seq(0,max(x),0.01)
xdens=lambda*exp(-lambda*xax)
hist(Rrexp,breaks=20,freq=FALSE,main=\"Rexp\") lines(xax,xdens,col=\"blue\
hist(myrexp(1000,0.5),breaks=20,freq=FALSE,main=\"myrexp\lines(xax,xdens,col=\"blue\
--------------------------------------------------------------------------------------------------------
Indirect Method
Reject & accept method
Generate random variable from N(μ,𝜎2)
-------------------------------------------------------------------------------------------------------- i=1
y=rep(NA,100) while(i<=100) {
17
u=runif(1,0,1)
v=runif(1,-2,-2) M=1/(sqrt(2*pi)*0.25) test=dnorm(v,0,1)/(0.25*M) if(u y[i]=v i=i+1 } } hist(y) -------------------------------------------------------------------------------------------------------- Optimization Newton-Raphson Algorithm ##单参数场合(一维) --------------------------------------------------------------------- 日期4.18 x=runif(1000,1,1) y=-2*x+ rnorm(1000,0,1) df=function(y,x,beta) { dlogDens=sum((y-x*beta)*x) } d2f=function(x) { d2logDens=sum(-x^2) } n=1 converge=FALSE while(converge==FALSE) { beta[1]=2 beta[n+1]=beta[n]-1/d2f(x)* df(y,x,beta[n]) if(abs(beta[n+1]-beta[n])<0.0000001) 18 { converge=TRUE } n=n+1 } beta --------------------------------------------------------------------- >optimize():optimize(f = , interval =c() , ..., lower = min(interval), upper = max(interval), maximum = FALSE,…） f=function(x) x^3+2*x+1 optimize(f,c(-2,2)) >uniroot():基于二分法来计算方程根，当初始区间不能满足要求 时，会返回错误信息 f=function(x) x^3-2*x-1 uniroot(f,c(-10,10)) ##多参数场合(多维) >optim(par,fn):par为函数初值(要优化的参数组成的向量),fn是求最值的函数 默认求最小值，通过control=list(fnscale=-1)可求最大值，fnscale可对结果进行缩放e.g.fnscale=-2 极大似然估计（Maximum likelihood）：这种方法估计的参数将会使得对数似然值最大或者负的对数似然值最小 --------------------------------------------------------------------- ## Generate some data beta0 <- 1 beta1 <- 3 sigma <- 1 n <- 1000 x <- rnorm(n, 3, 1) 19 y <- beta0 +x*beta1 + rnorm(n, mean = 0, sd = sigma) plot(x, y, col = \"blue\##编写似然函数如下： logNormLikelihood=function(par,y,x) { } ##Optimizing process optimOut<-optim(c(1,3,1),logNormLikelihood,control=list(fnscale=-1),method=\"BFGS\ beta0Hat <- optimOut$par[1] beta1Hat <- optimOut$par[2] sigmaHat <- optimOut$par[3] yHat <- beta0Hat + beta1Hat*x ##Use c(0,-1,0.1) $par [1] 1.184566 2.936066 170.831848 $value [1] -6059.635 $counts function gradient 109 100 $convergence ##Comparison with OLS [1] 1 $message NULL 20 beta0=par[1] beta1=par[2] sigma=par[3] mean=beta0+x*beta1 logDens=dnorm(x=y,mean=mean,sd=sigma,log=TRUE) loglikelihood=sum(logDens) return(loglikelihood) ##Use c(1,1,1) $par [1] 1.0802672 2.9678528 0.9769752 $value [1] -1395.644 $counts function gradient 142 75 $convergence [1] 0 ##0 indicates successful completion $message NULL lm(formula,data):结果是一个包含回归信息的列表，它包含以下信息： coefficients：回归系数（矩阵） residuals：返回模型残差（矩阵）等 ##Use OLS method to estimate a linear model given y and x myLM <- lm(y~x) myLMCoef <- myLM$coefficients yHatOLS <- myLMCoef[1] + myLMCoef[2]*x ## Plot the OLS’s result: points(sort(x), yHat[order(x)], type = \"l\ ## Plot the MLE’s result: points(sort(x), yHatOLS[order(x)], type = \"l\pch = 20) --------------------------------------------------------------------- 导入数据 read.table() :导入数据后生成数据框 read.table(“clipboard”,header = T) 注意只能传输显示在屏幕上的数据 容易损失精度 names()：查看变量 str()：查看数据框中每个变量的属性，及时发现错误的变量类型 > str(Squid) 'data.frame': 2644 obs. of 6 variables: $ Sample : int 1 2 3 4 5 6 7 8 9 10 ... $ Year : int 1 1 1 1 1 1 1 1 1 1 ... $ Month : int 1 1 1 1 1 1 1 1 1 2 ... $ Location : int 1 3 1 1 1 1 1 3 3 1 ... $ Sex : int 2 2 2 2 2 2 2 2 2 2 ... $ GSI : num 10.44 9.83 9.74 9.31 8.99 ... 实际中，三者应该结合一起使用. library(XLConnect) 21 Namelist=loadWorkbook() readWorksheet() 读取一个Excel文件的最好方式，就是在Excel中将其导出为一个逗号分隔文件(csv)，并使用read.csv()导入R中。也可以使用RODBC包来访问Excel文件。电子表格的第一行应当包含变量/列的名称 --------------------------------------------------------------------- det(x, ...) determinant(x, logarithm = TRUE, ...) x <- matrix(1:4, ncol = 2) > unlist(determinant(x)) modulus sign 0.6931472 -1.0000000 eigen():A的特征值和特征向量。 A=matrix(1:4,2,2) 若y <- eigen(A)，那么：y$val是A的特征值；y$vec是A的特征向量 prod(y$val) =det(A) prod(y$val)==det(A) abs (prod(y$val)-det(A)) all.equal() B=matrix(1:25,5,n) upper.tri(B) upper.tri(B,diag=TRUE) B[upper.tri(B,diag=TRUE)] whichcols=ceiling((1:25)/5) whichrows=(1:25)%%5 whichrows[whichrows==0]=5 B[whichrows>=whichcols] 22 kronecker(X, Y, FUN = \"*\"...)or X%x%Y Kronecker Products Stacking Vectors and Matrixes Description Stacks either a lower triangle matrix or a matrix. Usage vec(x) vech(x) Arguments x a numeric matrix qr分解 SVD分解 ------------------------------------------------------------------------------------------------------ New1 ------------------------------------------------------------------------------------------------------ Am+n的广义逆A+: A+A=I A+=(A’A)-1A’ β=(X’X)-1X’y= X+y=R-1QTy 23 三种方法比较如下： ------------------------------------------------------------------------------------------------------ y=matrix(1:5) X=cbind(1,rnorm(5),runif(5)) betaHat1=solve(t(X)%*%X)%*%t(X)%*%y betaHat1 XQR=qr(X) Q=qr.Q(XQR) R=qr.R(XQR) XPlus=solve(R)%*%t(Q) XPlus%*%X ##应该为单位阵 betaHat2=XPlus%*%y XSVD=svd(X) U=XSVD$u V=XSVD$v Sigma=diag(XSVD$d) Sigma1=diag(1/XSVD$d) XPlus1=V%*%Sigma1%*%t(U) betaHat3=XPlus1%*%y cbind(betaHat1,betaHat2,betaHat3) Review1 ------------------------------------------------------------------------------------------------------ mycbind=function(...) ##...表示参数的集合 { args=list(...) print(args) 24 nargs=length(args) out=NULL for (i in 1:nargs) { out=cbind(out,args[[i]]) } return(out) } mycbind(1:3,2:4,3:5) ------------------------------------------------------------------------------------------------------ A=matrix(1:24,4,6) A[,1,drop=FALSE] A[2:3,3:4] A[c(1,4),c(3,4)] B=cbind(rep(rep(TRUE,4),3),rep(rep(FALSE,4),3)) B=cbind(matrix(rep(rep(TRUE,4),3), ,3),matrix(rep(rep(FALSE,4),3), ,3)) A[B] A>20 A[A>20] A[A>10 & A<15] ##中括号用来取元素，用于访问数据子集 A[c(T,T,T)] A[c(T,T,F)] which函数 A=c(2,3,6,1,7,4,2,7,8,5) >A[which.max(A)] [1] 8 >A[which.min(A)] [1] 1 > which(A==2) [1] 1 7 > A[which(A>6)] [1] 7 7 8 ------------------------------------------------------------------------------------------------------ B=array(1:24,c(2,3,4)) 25 A=list(1,2,3) B=list(2,4,8) A+B mapply(sum,A,B) lapply(A,function(x) x/2) A=matrix(1:24, ,6) A apply(A,2,mean) ------------------------------------------------------------------------------------------------------ Debugging in R interactive debugging e.g. for (i in 1:n) control flow e.g. if integration test e.g. fun1=function{ fun2=function{} } log files Statistical Process Control 检验产生的随机数是否符合特定分布，histogram or 矩E(x) D(x) ------------------------------------------------------------------------------------------------------ logsum=function(a,b) { if(a>0 & b>0) { log(a)+log(b) } else {stop(\"a and b must be positive!\")} } j=NA for(i in 1:10) { browser() j=j+1 print(j) Sys.sleep(1) 26 } Browse[1]>n ##看下一步 Browse[1]>c ##不停下来,运行直到再次遇到browse()或直接Browse[1]>回车 Browse[1]>Q ##跳出调试 if(i>8) browser() ##当条件满足才进行browser j=NA for(i in 1:10) { Q=1 c=2 n=3 browser() j=j+1 print(j) } Browse[1]>print(Q) ##当变量中有n,c,Q时，打印的方法 option ##全局变量设置 options(browserNLdisablesd=FALSE) ##browserNLdisabled: logical: whether newline is disabled as a synonym for \"n\" is the browser. TRUE回车不再执行continue options(digits=3) a=1/3 try() is a wrapper to run an expression that might fail and allow the user's code to handle error-recovery. proc.time() system.time(expr): timing an R expression,实质是调用proc.time(),运行expr, 再次调用proc.time, 最后返回两次proc.time之差 数据的标准化与中心化以及scale详解 1.数据的中心化 所谓数据的中心化是指数据集中的各项数据减去数据集的均值。 例如有数据集1, 2, 3, 6, 3，其均值为3,那么中心化之后的数据集为 27 1-3,2-3,3-3,6-3,3-3,即：-2,-1,0,3,0 2.数据的标准化 所谓数据的标准化是指中心化之后的数据在除以数据集的标准差，即数据集中的各项数据减去数据集的均值再除以数据集的标准差。 例如有数据集1, 2, 3, 6, 3，其均值为3,其标准差为1.87，那么标准化之后的数据集为(1-3)/1.87,(2-3)/1.87,(3-3)/1.87,(6-3)/1.87,(3-3)/1.87, 即：-1.069,-0.535,0,1.604,0 数据中心化和标准化的意义是一样的，为了消除量纲对数据结构的影响。 在R语言中可以使用scale方法来对数据进行中心化和标准化： #限定输出小数点后数字的位数为3位 > options(digits=3) > data <- c(1, 2, 3, 6, 3) #数据中心化 > scale(data, center=T,scale=F) [,1] [1,] -2 [2,] -1 [3,] 0 [4,] 3 [5,] 0 attr(,\"scaled:center\") [1] 3 #数据标准化 > scale(data, center=T,scale=T) [,1] [1,] -1.06904 [2,] -0.53452 [3,] 0.00000 [4,] 1.60357 [5,] 0.00000 attr(,\"scaled:center\") [1] 3 attr(,\"scaled:scale\") [1] 1.8708 scale方法中的两个参数center和scale的解释： 1.center和scale默认为真,即T或者TRUE 2.center为真表示数据中心化 3.scale为真表示数据标准化 数据的标准化和正态化 log transformation in R 标准化standardization 为什么要进行数据的标准化standardization？ http://www.biomedware.com/files/documentation/boundaryseer/Preparing_data/Why_standar 28 dize_variables.htm 正态化normalization 以下是正态化的一种算法 http://stn.spotfire.com/spotfire_client_help/norm/norm_scale_between_0_and_1.htm 在R中使用这个函数可以实现 range01 <- function(x){(x-min(x))/(max(x)-min(x))} 参考： http://r.789695.n4.nabble.com/Changing-data-scales-td3759180.html http://stackoverflow.com/questions/5700907/auto-scale-data-in-r http://stackoverflow.com/questions/5468280/scale-a-series-between-two-points-in-r http://stackoverflow.com/questions/15468866/scaling-a-numeric-matrix-in-r-with-values-0-to-1 标准化和正态化的异同 http://www.dataminingblog.com/standardization-vs-normalization/ 在R中对数据标准化并作图： http://rss.acs.unt.edu/Rdoc/library/base/html/scale.html http://rss.acs.unt.edu/Rdoc/library/pls/html/stdize.html http://www.pmc.ucsc.edu/~mclapham/Rtips/ordination.htm 多元统计分析中的主成分分析： http://www.math.pku.edu.cn/teachers/lidf/docs/statsoft/html/sas-5/multivar.html 本来想在中文网页中找一些有用的信息，结果搜到的内容是下面这个链接： http://wenwen.soso.com/z/q406168868.htm 表示很无语。 log transformation in R in Ecology - Community data: data<-read.csv(\"File.csv\ trans<-log(data) attach(trans) 参考： http://www.sci.muni.cz/botany/zeleny/wiki/anadat-r/doku.php?id=en:pca http://msenux.redwoods.edu/math/R/TransformingData.php http://cc.oulu.fi/~jarioksa/softhelp/vegan/html/decostand.html http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/ape/html/pcoa.html http://biol09.biol.umontreal.ca/PLcourses/ http://biol09.biol.umontreal.ca/PLcourses/Section_7.7_Transformations.pdf 29 因篇幅问题不能全部显示，请点此查看更多更全内容