Tag Archives: R

1、线性拟合

#生成测试数据
x = seq(-5,5,0.1)
y = 3*x^2+6*x+9+rnorm(length(x))*3;

#把x^2用I来标记成一个变量
#然后进行线性拟合
z=lm(y~I(x^2)+x)

#绘制数据点
#及拟合曲线
plot(x,y)
lines(x,fitted(z))

2、局部多项式回归拟合

#生成测试数据
x = seq(-5,5,0.1)
y = 3*x^2+6*x+9+rnorm(length(x))*3;

#局部多项式回归拟合
z=predict(loess(y~x))

#绘制数据点
#及拟合曲线
plot(x,y)
lines(x,z)

#lowes默认使用局部多项式回归拟合
#z1=lowess(x,y)
#lines(z1)

3、非线性最小二乘拟合

#生成测试数据
x = seq(-5,5,0.1)
y = 3*x^2+6*x+9+rnorm(length(x))*3
ds <- data.frame(x=x, y=y)

#进行非线性最小二乘拟合
f=function(x, a, b, c, d) {a+b*x+c*x^2}
z=nls(y~f(x, a, b, c), data=ds, start=list(a=9, b=6, c=3))

#输出拟合结果
summary(z);

#绘制数据点
#及拟合曲线
plot(x,y)
lines(x,fitted(z))

1、分层聚类

#载入R自带的测试数据
#data(iris)
#attach(iris)
inData = iris[,1:4]

#计算距离矩阵,并绘图
inData.dist = dist(inData)
#inData.dist = dist(inData,method='euclidean')
heatmap(as.matrix(inData.dist), labRow = F, labCol = F)

#进行分层聚类
#并绘图
#inData.hc <- hclust(inData.dist)
#inData.hc <- hclust(inData.dist,method='ward')
plot(inData.hc, labels = FALSE, hang = -1)

#标识聚类结果，结果设为3类
rect.hclust(inData.hc, k = 3)

#将Tree进行分组
inData.groups <- cutree(inData.hc, 3)
#输出结果表格
table(inData.groups, Species)

#进行降维处理
#绘图对比结果
#形状是正确的数据
#颜色为聚类后的数据
mds=cmdscale(inData.dist,k=2,eig=T)
x = mds$points[,1]
y = mds$points[,2]
library(ggplot2)
p=ggplot(data.frame(x,y),aes(x,y))
p+geom_point(size=3,alpha=0.8,aes(colour=factor(inData.groups),shape=iris$Species))

2、K值聚类

#载入R自带的测试数据
#data(iris)
#attach(iris)
inData = iris[,1:4]

#计算距离矩阵,并绘图
inData.dist = dist(inData)
#inData.dist = dist(inData,method='euclidean')
heatmap(as.matrix(inData.dist), labRow = F, labCol = F)

#进行K值聚类
inData.kc <- kmeans(inData.dist,centers=3)

#绘图对比结果
#形状是正确的数据
#颜色为聚类后的数据
library(ggplot2)
x=inData[c("Sepal.Length")]
y=inData[c("Sepal.Width")]
p=ggplot(data.frame(x,y),aes(x,y))
p+geom_point(size=3,alpha=0.8,aes(colour=factor(inData.kc$cluster),shape=iris$Species))

首先要安装线性规划扩展包

install.packages("lpSolve")

例1：
某工厂甲、乙两种产品，每件甲产品要耗钢材2kg、煤2kg、产值为120元；每件乙产品要耗钢材3kg，煤1kg，产值为100元。现钢厂有钢材600kg，煤400kg，试确定甲、乙两种产品各生产多少件，才能使该厂的总产值最大？

解：

# 优化方向: 
# maximize
# 目标函数：
# f(X1，X2)＝120X1＋100X2
# 约束条件：
# 2X1＋3X2<=600
# X1＋X2<=400

#引用lpSolve
library (lpSolve)

#求解
f.obj <- c (120, 100)
f.con <- matrix ( c (2, 3, 1, 1), nrow = 2, byrow = TRUE )
f.dir <- c ( "<=" , "<=" )
f.rhs <- c (600, 400)
lp.result <- lp ( "max" , f.obj, f.con, f.dir, f.rhs)

#输出结果
lp.result
#Success: the objective function is 36000 

#输出解
lp.result$solution
#[1] 300   0

例2：
某工厂要做100套钢架，每套用长为2.9m，2.1m，1.5m的圆钢各一根。已知原料每根长7.4m，问：应如何下料，可使所用原料最省？(可选方案如下)

解：

# 优化方向: 
# minimize
# 目标函数：
# f(X1,X2,X3,X4,X5)＝0X1＋0.1X2+0.2X3+0.3X4+0.8X5
# 约束条件：
# X1＋2X2 + 0X3 + X4 + 0X5>=100
# 0X1＋0X2 + 2X3 + 2X4 + 1X5>=100
# 3X1＋X2 + 2X3 + 0X4 + 3X5>=100

#引用lpSolve
library (lpSolve)

#求解
f.obj <- c (0, 0.1, 0.2, 0.3, 0.8)
f.con <- matrix ( c (1, 2, 0, 1, 0, 0, 0, 2, 2, 1, 3, 1, 2, 0, 3), nrow = 3, byrow = TRUE )
f.dir <- c ( ">=", ">=", ">=" )
f.rhs <- c (100, 100, 100)
lp.result <- lp ( "min" , f.obj, f.con, f.dir, f.rhs)

#输出结果
lp.result
#Success: the objective function is 10

#输出解
lp.result$solution
#[1] 100   0  50   0   0

1、数据文件
1.1、InA.txt

1,2,3
4,5,6
7,8,9

1.2、InB.txt

1
2
3
4
5
6
7
8
9

1.3、InC.csv

col1,col2,col3
1,2,3
4,5,6
7,8,9

3、数据读取
3.1读行取文本转为矩阵

filename="C:\\TestData\\InA.txt"
matrixA <- as.matrix(read.table(filename, header=FALSE, sep = ",",as.is=TRUE))

3.2按Cell读取文本转为矩阵

filename="C:\\TestData\\InB.txt"
matrixB <- matrix(scan(filename),ncol=3,byrow=TRUE)

3.3直接读取CSV

filename="C:\\TestData\\InC.csv"
matrixC =read.csv(filename,header=FALSE)

4、数据输出
4.1输出为文本文件

filename="C:\\TestData\\resultA.txt"
write.table(resultA, file = resultA, row.names = F, quote = F)

4.2输出为CSV文件

filename="C:\\TestData\\resultB.csv"
write.csv(matrixResult,filename)

5、设置工作路径

setwd("C:\\TestData")

1向量
1.1定义向量

##定义向量，下面两种方式结果相同
##c(..., recursive = FALSE) 
x=1:12
x=c(1,2,3,4,5,6,7,b,9.10)

1.2向量与数字四则运算

##+-×/
x=1:12
y=10
z=x*y

1.3向量的内积

##%*% crossprod
x=1:5
y=1:5
z=x%*%y
z1=crossprod(x)
z2=crossprod(x,y)

1.4向量的外积

##x%o%y tcrossprod
x=1:5
y=1:5
z=x%o%y
z1=tcrossprod(x)
z2=tcrossprod(x,y)
z3=outer(x,y)

1.5矩阵转向量

vec<-function (x){
          t(t(as.vector(x)))
}
vech<-function (x){
          t(x&#91;lower.tri(x,diag=T)&#93;)
}
x=matrix(1:16,nrow=4,ncol=4,byrow=TRUE)
vec(x)
vech(x)
&#91;/code&#93;

1.6生成滞后序列
&#91;code lang="R"&#93;
library(fMultivar)
x=1:20
tslag(x,1:4,trim=F)
tslag(x,1:4,trim=T)
&#91;/code&#93;

<strong>2数组</strong>
2.1定义数组
[code lang="R"]
##向量转数组，按列转
##dim(x)
dim(x)=c(3,4);
##或者直接定义多维数组
x=array(1:12,dim=c(3,4))

3矩阵
3.1定义矩阵

##matrix(data=NA，nrow=1，ncol=1，byrow=FALSE，dimnames=NULL)
##data项为必要的矩阵元素，nrow为行数，ncol为列数，nrow与ncol的乘积应为矩阵元素个数，byrow项控制排列元素时是否按行进行，dimnames给定行和列的名称
x=matrix(1:16,nrow=4,ncol=4,byrow=TRUE)
##行数
nrow(x)
##列数
ncol(x)
##上三角阵
C=x
C[lower.tri(A)]=0
##下三角阵
C=x
C[upper.tri(A)]=0
##元素的行与列
row(x)
col(x)
D=x
D[row(x)<col(x)&#93;=0
D=x
D&#91;row(x)>col(x)]=0

3.2定义对角矩阵

##diag
A=c(1:5)
##对角阵
B=diag(A)
##获取对角元素
C=diag(B)
##生成单位对角阵
D=diag(5)

3.3矩阵与数字四则运算

##+-*/
A=matrix(1:9,nrow=3,ncol=3,byrow=TRUE)
B=10
C=A+B
D=A-B
E=A*B
F=A/B
##行和
rowSums(A)
apply(A,1,sum)
##行平均值
rowMeans(A)
apply(A,1,mean)
##行方差
apply(A,1,var)
##列和
colSums(A)
apply(A,2,sum)
##列平均值
colMeans(A)
apply(A,2,mean)
##列方差
apply(A,2,var)

3.4矩阵与矩阵加减法

##+-
A=matrix(1:9,nrow=3,ncol=3,byrow=TRUE)
B=matrix(10:18,nrow=3,ncol=3,byrow=TRUE)
C=A+B
D=A-B

3.5矩阵与矩阵乘法

##%*%
A=matrix(1:9,nrow=3,ncol=3,byrow=TRUE)
B=matrix(10:18,nrow=3,ncol=3,byrow=TRUE)
C=A%*%B
##t(A)%*%B;
C1=crossprod(A,B)
##A%*%t(B)
C2=tcrossprod(A,B)
##Hadamard积
D=A*A

3.6矩阵转置

##t(x)
A=matrix(1:6,nrow=2);
B=t(A)

3.7矩阵行列式

##det(x, ...)
A=matrix(1:9,nrow=3);
B=det(A)

3.8矩阵求逆

##solve
A=matrix(c(1,2,3,2,2,1,3,4,3),nrow=3,byrow=TRUE);
B=det(A)
C=solve(A)

##需要strucchange包,A'A的逆
library(strucchange)
C1=solveCrossprod(A,method="qr")
C2=solveCrossprod(A,method="chol")
C3=solveCrossprod(A,method="solve")
C4=solve(crossprod(A,A))

3.9求解线性方程

##solve
A=matrix(c(1,2,3,2,2,1,3,4,3),nrow=3,byrow=TRUE);
B=c(1,2,8)
C=solve(A,B)

##对于系数为三角阵的情况，可以用 backsolve及fowardsolve函数求解
##这一段有些错误，时间不够了，下次修复
A=matrix(1:9,3,3)
V=c(1,2,3)
B=A
B[upper.tri(B)]=0
backsolve(A,V,upper.tri=F,transpose=F)
forwardsolve(A,V,upper.tri=F,transpose=F)
##why not the same?
##something is wrong
solve(B,V)

3.10特征值及特征向量

##eigen
A=matrix(c(1,-2,2,-2,-2,4,2,4,-2),nrow=3,byrow=TRUE);
B=eigen(A)
##结果应该为对角阵
C=solve(B$vectors)%*%A%*%B$vectors

4.矩阵分解
4.1奇异值分解

##A=U%*%D%*%t(V)，其中U, V是正交阵，D为对角阵，也就是矩阵A的奇异值
##svd
A=matrix(c(1,0,1,0,0,1),nrow=3,byrow=TRUE);
B=svd(A)
C=B$u%*%diag(B$d)%*%t(B$v)

4.2QR分解

##设A为m*n矩阵，如果存在m*m酉矩阵Q（即Q(H)Q=QQ(H)=I）和m*n阶梯形矩阵R，使得A=QR，那么此分解称为QR分解
##qr
A=matrix(c(1,1,1,2,-1,-1,2,-4,5),nrow=3,byrow=TRUE);
B=qr(A)
Q=qr.Q(B)
R=qr.R(B)
C=Q%*%R

4.3Schur分解

##A=USU(H)，其中U是标准的正交矩阵（即满足UU(H)=I），S为上三角矩阵，并且对角线上的元素为A的特征值。
##Schur
library(Matrix);
A=matrix(c(1,1,1,2,-1,-1,2,-4,5),nrow=3,byrow=TRUE);
B=Schur(A, vectors=TRUE);
C=B$Q%*%B$T%*%t(B$Q)

4.4Cholesky分解

##对任意的正定矩阵A，存在上三角矩阵R，使A=t(R)%*%R，则称为A的Cholesky分解
##chol
A=matrix(c(4,-1,1,-1,4.25,2.75,1,2.75,3.5),nrow=3,byrow=TRUE);
B=chol(A);
C=t(B)%*%B;
##求行列式
D=prod(diag(chol(A))^2)
##求逆
E=chol2inv(chol(A))

今天看了一下R语言，感觉还是很强大的，做了几个小例子，比如这几个画心形曲线的：

dat<- data.frame(t=seq(0, 2*pi, by=0.1) )
xhrt <- function(t) 16*sin(t)^3
yhrt <- function(t) 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t)
dat$y=yhrt(dat$t)
dat$x=xhrt(dat$t)
with(dat, plot(x,y, type="l"))
with(dat, polygon(x,y, col="hotpink"))

MASS::eqscplot(0:1,0:1,type="n",xlim=c(-1,1),ylim=c(-0.8,1.5))
curve(4/5*sqrt(1-x^2)+sqrt(abs(x)),from=-1,to=1,add=TRUE,col=2)
curve(4/5*-sqrt(1-x^2)+sqrt(abs(x)),from=-1,to=1,add=TRUE,col=2)

这两个的区别，主要是公式不一样，还有其他的公式，可以查看这里：
Heart Curve

还有用符号直接画的：

hearts <- expression(symbol('\251'))
plot( 0, xlim=c(0,2), ylim=c(0,2), type="n" )
text( 1, 1, hearts , col="red", cex=30 )  

plot(1, 1, pch = "♥", cex = 20, xlab = "", ylab = "", col = "red")