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))
Posted in R | Tagged

R语言做聚类分析

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))

ClusterAnalysisH.png

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))

ClusterAnalysisK.png

Posted in R | Tagged

R语言做简单线性规划

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

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,问:应如何下料,可使所用原料最省?(可选方案如下)

可选方案 方案1 方案2 方案3 方案4 方案5
2.9m 1 2 0 1 0
2.1m 0 0 2 2 1
1.5m 3 1 2 0 3
合计 7.4 7.3 7.2 7.1 6.6
剩余料头 0 0.1 0.2 0.3 0.8

解:

# 优化方向: 
# 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
Posted in R | Tagged

R从文本文件读取矩阵

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")
Posted in R | Tagged

R语言矩阵运算

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定义数组

##向量转数组,按列转
##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))
Posted in R | Tagged

Heart Curve

今天看了一下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")
Posted in R | Tagged