par(family="Sarasa Gothic CL")#这个命令运行后就可以使用中文字体了
a<-3+7
b<-8
data("iris")
as.data.frame()
#加注释
control+enter #运行当前行
setwd('//Users//wangxinran//Desktop//R语言')#将括号内更改为默认目录，
#所有文件将默认储存到这个文件夹，点上栏文件夹符号即可查看该文件夹，注意
#要把单/改为双//
getwd()#查看当前目录是什么
install.packages('e1071')#安装一个包e1071
library('e1071')#打开安装包
x<-c(1,2,3,4,5)#写向量
y<-x^2
plot(x,y)#画图
?plot#查帮助，问电脑plot是做什么的
plot(cars)
plot(x,y,type='l',col='red')#让图是一条红色的线
x=c(1,2,3,4,5)
y=c(6,7,8,9,0)
mean(x)#求x的平均值
x<-c(1,2,3,4,5,2,4,3)
y<-c(x^2)
plot(x,y)
q<-c(x,x,7,8,5)#q里面装了x
mean(q)
q[5]#取q里的第五项，注意不安shift，是方括号，换成英文
q[3:5]#取第三项到第五项
A<-9
a<-8
#R语言的变量或者文件命名区分字母大小写，所有的字母和数字都可以用，
#后面还可以跟.或者下划线_，用.或者字母开头，但是若用.开头，
#后面紧接着的一位不能是数字，命名长度不限
wxr.x1999_R <- c(1,3,5,8)#注意c不可以大写
#写完一个代码或命令要用；隔开或者换行
#{}在执行循环体或函数语序体时会用
objects()#查看当前使用过的变量
rm(a,b)#删除某几个变量a，b
dev.new()#新开一个窗口，我们再执行画图，图像就在新窗口显示
assign("m",c(1,2,3))#把向量赋给m，也是写向量的一种方式
length(m)#向量m的长度
#向量运算，两个向量必须长度相当
n<-m+g
n<-m*g
n<-m/g
n<-m^g
sqrt(m)#对m开平方
exp(m)#e的m次方
log(m)#还可以取sin，cos
log(9,3)#log(x,base),后面是底，前面是真数
#复制rstudio中的图像到word，点击右栏polts，export，copy to clipboard，
#再黏贴到word
#control+l清空
#control+s 储存（笔记）
#control+加减号，可调整页面字体大小
data()#可显示R自带的数据库
data(CO2)#可调出数据库中的co2来看，在右栏中单击即可显示表格
#当上一行出现错误时，按键盘上箭头复制上一行，修改
#注意数字跟字母相乘时，也一定要有*号y=-1*x
#赋值符号前后有点空格，方便观看
#注意大小写，向量别忘了c
data.frame#数据框
data("iris")
iris<-iris
max(iris$Sepal.Length)#查看表中某一列或行的最大值
range(iris$Sepal.Width)#查看表中某一行或列的范围
prod(x,y) # 求乘积
var(x)# 求方差
x <- 1:30
x <- 2*1:15#产生向量1～15，且每一项都乘2
x <- seq(from=3,by=0.5,length=5)#seq 产生规则序列，by是公差
x <- seq(-10,10,by=0.1)#另一种写法,0.1的0可以省略
y=sin(x)
plot(x,y)
seq(2:8)#表示数列从1开始到8-2+1，即1 2 3 4 5 6 7
m <- 30
n <- (m-1):60
x <- c(-3,4,6.5,7,0)
y <- rep(x,times=5)#把x整体重复5遍
y <- rep(x,each=5)#把x每个数挨着重复5遍
x <- 4
y <- 5
if(x==4) z=5+6#双==表示判断谁等不等于谁,若if条件不满足，则找不到z的值
if(x<8) h=10-x
if(x<8&y>1)t=x+2*y#&“与”，表示前后两个条件同时成立
if(!(x<5))z=1#!"非“，表示否定后面的条件，x>=5,但输入的x=4，
#所以无法输出z
s <- 4
temp <- s<3#temp表示判断后面命题是否为真，再输入temp回车，系统会告诉你
#True，Flase
x <-5
temp <- x>3
if(temp)z=x+6#即如果temp出来是True，就能运算出z，F就算不出来
q <- 5
temp <- q>8
if(temp)p=11
x <-5
y <-6
if(x==6|y>3)z=11#|“或”，表示即前面或者后面只要有一个成立即可
setwd('//Users//wangxinran//Desktop//R语言')#先用setwd把目录换成默认
teens <- read.csv("snsdata.csv")#导入目录文件，把文件数据命名为teens
str(teens)#紧凑的表示对象的内部结构，即显示teens的内容
table(teens$gender)#统计表格中gender列每一种元素的个数
teens <- read.csv(file="snsdata.csv",header=T)
#header=T有表头（即第一行的变量名），header=F没表头，file是格式，没啥用
is.na(teens$age)#看age这个变量哪一列数据缺失
#$表示引用数据中的某一列
z <- c(1:3,NA)
is.na(z)#看z哪里缺失数据，缺失的地方是True，不缺失Flase
#缺失有两种，第一是数据直接是缺失的，第二种是数据包在运算过程中缺失，
#NaN 表示not a number，eg：0/0；Inf—Inf，当数据中出现这两个时，
#用is.nan（）或is.na（）都能判断数据在运算过程中是缺失的，
#但is.nan（）只能判断运算过程中的缺失
#而is.na（）两种缺失都能判断
x <- c(3,4,5,Inf-Inf,0/0)
is.na(x)
is.nan(x)
q <- c("x","y","z")#向量q可以是一个字符或字符串
q <- c('I want to eat','no')
labs <- paste(c("X","Y"),1:10)#将XY与数字1到10拟合（连接），赋值给lab
labs <- paste(c("X","Y"),1:10,sep='+')#sep表示拟合符号用加号+拟合
labs <- paste(c("X","Y","Z"),1:10,sep='&')#没有sep表示默认用空格
#作为中间的拟合符号
paste(1:12,c("st","nd","rd",rep("th",9)))#rep 重复9次
x <-c("hello","world","!")
nchar(x)#判断字符向量中，每一个元素字符串中字母的个数
DNA <- "AtGCtttACC"
tolower(DNA)#把所有的字母换成小写
toupper(DNA)#把所有的字母换成大写
chartr("Tt","Uu", DNA)#改变字母，把T变成U，t变成u
chartr("TtA", "Uua", DNA) # T换成U, t换成u, A换成a
X<-c("He say:\'hello', and then go")
X<-c("He say:'hello', and then go")
substr("abcdef",start=2,stop=4)#截取字符串string中的一部分sub，从第二个到第四个
#index vector 向量的索引，即在向量里取一部分，形式：向量[ ]
x<-c(-3:10,NA,NA)
x[is.na(x)] <- 0#把x里的NA改成0
y <- x[!is.na(x)]#x不是NA的地方赋给y
(x+1)[(!is.na(x)) & x>0] -> z#把x里不是NA并且大于0的数取出来，加1，最后赋给z
x[(!is.na(x)) & x>0]+1 -> z#与上面一样的
x[1:3]#看x第一到三个元素,[ ]里填想索引的子集
x<-c(-3:10,NA,NA,3:5,NA,1)
x[is.na(x)] <- 3
y<-c(-4:5)
y[y < 0] <- -y[y < 0]
y <- abs(y)#把y的每一项都变成正的
teens <- read.csv("snsdata.csv")
teens$gender <- ifelse(is.na(teens$gender),'I see',teens$gender)
#把teens里gender这一列里缺失的地方都换成I see，其余不变
#ifelse，如果（缺失，换成I see）否则（不变）
c("x","y")[rep(c(1,2,2,1), times=4)]#把字符“x","y"按1221的顺序重复思辨
x<-c(-3:10)
y <- x[-(1:5)]#除去x里的第1到5项
fruit<- c(5, 10, 1, 20)
names(fruit) <- c("orange", "banana", "apple", "peach")
fruit
lunch<- fruit[c("apple","orange")]#索引里面的apple，orange那两列
#想把上面的fruit表格存下来，首先把你想存的地址先换成当前目录setwd（）
#fruit表格有两种储存形式
write.table(fruit, file = "fruit.txt", row.names = F, quote = F)#txt形式
write.table(fruit, file = "fruit.txt")#有表头
install.packages('ggplot2')#作图包
library('ggplot2')#只需要安装一次，以后只要library用就行
pie(iris$Sepal.Length)#画饼状图
plot(iris$Sepal.Length)#散点图
gupiao <- read.csv('ss_stock.csv')
plot(gupiao)
plot(gupiao$Open,type='l')
#factor型变量,分类别或者分层次的，比如iris里的species
levels(iris$Species)#看species的每一层是什么
str(iris)#species分3类"setosa","versicolor",..我们可以分类画图
plot(iris$Sepal.Length,iris$Species)
library(MASS)#用require也可以加载一个包
a1 <- Cars93
plot(a1$Origin)
dev.new()
plot(a1$Origin,a1$Price)#箱线图，中间的黑线是中位数，可查资料看
plot(~iris$Sepal.Length+iris$Petal.Length)#画二者的关系图，注意～别忘了
#～表示取一部分数据，看两两之间有什么关系
plot(iris$Sepal.Length~iris$Sepal.Width+iris$Petal.Length)#～前是因变量纵坐标
#～后面是自变量，这一个命令画了两个图（sepallength分别跟width跟length的图），
#点箭头可以查看前面的图
#recursive递归结构，list里面再装一个list
x <- c(3:10,NA,NA)
mode(x) #看一个对象的属性
#object对象，可以是向量，也可以是list列，向量只能是简单形式eg：numeric、character
#多个模式不能交杂在一个vector里，但list里可以有多种形式
y <- as.character(x)#as.character，改变类型，把x从数字型转化成字符型as.numeric( )
#as.logical( )as.charactor( )as.matrix( )as.dataframe( )as.list( )
a <- numeric()#写一个空向量a
a[3] <- 7#填充a的第三个位置为7
length(a) <- 3 #a的长度就变为3了，用这种办法可以改变向量的长度
a[4:8] <- c(6,7,10,-2) #再把a的4～8位用向量c填充
length(a) <- 5# 截取，使a的长度变为5，只保留前五个元素
a <- a[-6:-8]#去掉a里的第六到底八个元素
x <- cbind(a=1:3,pi=pi)#列合并
x <- rbind(a=1:3,b=7:9)#按行合并
options(digits = 3)#改变小数的有效位数
attributes(x)#显示x的维度
#dim(x)看x行、列的维度，并写成一个向量
#结果出现“dim 3 2”:表示有三行两列
#dimnames[1]:看行的名字
#dimnames[2]:看列的名字
#dim(M)[1]:看数据框M行的维度（有几行）
#dim(M)[2]:看数据框M列的维度（有几列）
x <- 1:24
attr(x,'dim') <- c(3,8)#把x的维度变为三行八列
row.names(x) <- c('row1','row2','row3')#改变行的名字
colnames(x) <- c('col1','col2','col3','col4','col5','col6','col7','col8')
colnames(x) <- paste(c("col"),1:8,sep = "")#这样就不用一个一个输入col了
attr(x,"dimnames")<-list(paste("row",1:3,sep = ""))#另一种改变行名字的办法
class(x)#class看object的类型eg matrix，数组等，比mode范围更大
typeof#看是整形，范围最小最细致
as.integer(digits)#把小数变成整形
gl(2,5)#2代表两个水平，5代表每一个水平有5个，gl产生一个factor型，即因子型变量
#因子型，比如iris里的species就是因子型，有三层
a <- gl(2,5,label=c("Male","Female"))
state <- c("tas","sa","qld")
class(state)
statef <- factor(state)#把state换成因子变量
statef <- as.factor(state)#同上
class(statef)
is.factor(statef)
scores <- scan()#scan()在键盘上输入数据赋给scores，再看scores就有你输入的数据
#> scores <- scan()
#1: 68 89 98 100 67 28#输入你的数据
#7: #回车
# Read 6 items #系统自动出现
#> scores #查看scores
#[1] 68 89 98 100 67 28
cut(scores,breaks=c(0,70,80,90,100))#把数据按0～70，70～80····来分段，
#就可以看到每一个数据在哪一层，还可以给每一层用字符命名labels=
scoresF <-cut(scores,breaks=c(0,70,80,90,100),labels=c("差","中","良","优"))
table(scoresF)#统计每一层分别有多少个
class(scoresF)#scoresF是因子型
sexs <- sample(c('M','F'),length(scores),replace = T)#随机抽样，构造性别向量
#随机产生一些性别，按照scores的长度，自动给每个数据配上性别
#replace=T“有放回的随机抽样？Ture”
table(sexs,scoresF)#统计男女同学在不同层次的人数，做一个table，sex为行scoresF为列
tapply(scores,sexs,mean)#计算男女同学的平均成绩“成绩要按照性别求平均”
tapply(scores,sexs,sd)#计算男女同学成绩的标准差“成绩要按照性别求标准差”
#因子型变量不能直接比较大小，因为默认是无序，需要先变成有序型的因子变量ordered=T
scoresF1 <- factor(scoresF,ordered = T)
scoresF1[1]<scoresF1[3]
#矩阵就是一个二维的数组
x <- 1:24
dim(x) <- c(3,8)#产生矩阵的方法（1）只能按列产生，先产生第一列。。。
attr(x,'dim') <- c(3,8)#产生矩阵的方法（2）
#产生数组方法如下：
dim(x) <- c(2,6,2)#两行六列的矩阵两个（必须相乘等于元素个数）
x <- 1:36
dim(x) <- c(2,2,3,3)#四维数组，3*3个2*2维的矩阵
#数组：2*2的矩阵，3个为一组，有3组
#先产生一个3*3的大矩阵，每一个元素都是一个2*2的矩阵
dim(x) <- c(4,9)
x[c(1,3,4),c(1,3,5,7,9)]#取第1，3，4row与第1，3，5，7，9column交叉项
x[1:3,2:8]#矩阵索引的另一种，取1～3行，2～8列
x[,2:8]#","表示都要，行都要，列取2～8
x[5]#x[n]可以按矩阵生成的顺序一个一个地索引，先索引完第一列，再开始第二列
x <- 1:36
dim(x) <- c(2,2,3,3)
x[1,1,1,1]#数组索引
x[,,2,3]#","表示全要
x[3]#数组也可以按数组产生的顺序一个一个索引,顺序索引
a <- 1:4
b <- 2:5
cbind(a,b)#按列产生矩阵
rbind(a,b)#按行产生矩阵
x <- 1:36
matrix(x,nrow = 9,byrow=T)#按行产生矩阵，9行
matrix(x,nrow = 4,byrow=F)#按列产生矩阵，4行
x <- 1:36
array(x,c(4,9))
#产生矩阵的方法有？
#dim, attr,cbind,rbind,matrix,array
y <- 25:48
y.matrix <- matrix(y,nrow=4,byrow=F)
x <- 1:24
dim(x) <- c(4,6)
x+y.matrix
x*y.matrix#对应位置元素相乘,要求两个矩阵的行数列数都相同
L <- 1:24
H <- 25:48
L%*%H #两个向量内积：对应元素相乘相加，是一个数
L%o%H #两个向量外积：Anx1 x A1xn=Anxn,是一个向量，自动把前面那个变成一列的向量
outer(L,H,'*')#计算外积的另一种方式，*可以改成任意一种运算符号，甚至是一个函数
f <- function(x,y){cos(x)/sin(y)}
outer(L,H,'f')#L,H按照f做外积运算
#L里第一个元素，与H里第一个元素当成xy去做f运算，作为新矩阵的第一个元素
f <- function(x,y){min(x,y)}
outer(L,H,'f')
y <- 25:48
y.matrix <- matrix(y,nrow=4,byrow=F)
x <- 1:24
dim(x) <- c(6,4)
x%*%y.matrix#线代里的矩阵乘法
dim(y.matrix)
y.matrix <- t(y.matrix)#t()做矩阵转至
dim(y.matrix)
a <- 1:24
x <- array(a,c(2,3,4))
aperm(x)#数组转至 c（4，3，2）变成两个矩阵，以前的所有矩阵的第一行组成第一个矩阵
#以前所有矩阵的第二行挑出来组成第二个新的矩阵
A <- matrix(1:36,nrow = 6)
diag(A)#取对角线元素
sum(diag(A))#diagA所有元素相加
diag(3)#生成3x3阶单位矩阵
mean(A)#求A所有元素的均值
iris <- iris
mean(iris$Sepal.Length)#求Sepal.Length这一列的均值
colMeans(iris[,1:4])#同时求一到四列各自的均值
colMeans(iris[3:8,1:4])#同时对1～4列的3～8行求各自的均值
apply(iris[,1:4],1,mean)#对iris的所有行求均值，若中间换成2就是求列平均
#apply既可以对列求均值，也可以对行求均值
A <- 1:36
A <- matrix(A,nrow = 6,byrow = T)
B <- 37:72
B <- array(B,dim = c(6,6))
A*B#A B 对应元素相乘
A%*%B#A B做高代里的矩阵乘法
diag(A)#取出对角线元素
trace_A <- sum(diag(A))#trace_A一般用来命名矩阵的ji，就是对角线元素的和
diag(3)#生成3x3单位矩阵
mean(A)#所有元素求均值
colMeans(A)#每一列求均值
rowMeans(A)#每一行求均值
apply(A,2,mean)#每一列求均值,注意mean函数可以换，比如sum
apply(A,1,mean)#每一行求均值
#线性方程组A%*%x=b
#求解线性方程组solve(A,b)
x <- c(1,-2,3,-4,0,1,-1,1,1,3,0,1,0,-7,3,1)
A <- matrix(x,nrow=4,byrow = T)
b <- c(4,-3,1,-3)
x <- solve(A,b)#以行的形式展示结果
x <- solve(A)%*%b#以列的形式展示结果
#注意，只能求解方阵
#8e-16 8乘10的负16次方
A%*%x#可用于验证结果是否正确，若=b则正确
a <- c(1,-2,3,-4)
b <- c(0,1,-1,1)
c <- c(1,3,0,1)
d <- c(0,-7,3,1)
A <- rbind(a,b,c,d)#按行产生矩阵
b <- c(1,2,3,4)#注意不要写成列向量
aperm(b)#求数组b的转至
t(as.matrix(b))#把b改成矩阵，求矩阵转至
B[-4,]#去掉矩阵B中的某行或某列
#只有方阵才能求行列式，特征值，特征向量
B <- 37:72
B <- array(B,dim = c(6,6))
ev <- eigen(B)#求B的特征值特征向量
ev$values#只看eigen value
ev$vectors#只看eigen vector
typeof(ev)
ev[[1]]#只看list里第一项
ev[[2]]#只看list里第二项
#自己编个程序，求矩阵的正交变换
tapply(iris$Sepal.Length,iris$Species,mean)
tapply(a,b,mean)#按b求a的均值
A <- 1:36
A <- matrix(A,nrow = 6,byrow = T)
apply(A,1,sum)#对矩阵A按行求和
a <- sapply(1:3,function(x)x^2)#对A第1:3个元素按function变换，求出来a是个向量
#可以直接求和 求均值 sum(a)mean(a)
b <- lapply(1:3,function(x)x^2)#对A第1:3个元素按function变换，求出来是个list，
#list是一个杂框，想装什么都行
sapply(iris[,1:3],function(x)cos(x)/(sin(x)+exp(x)))
sapply(iris[,1:3],mean)
iris1 <- iris[1:4,1:4]
sapply(iris1,mean)#每一列求均值
b1 <- lapply(iris[,1:3],mean)
b1[[1]]#看b1的第一个位置
class(b1)
b2 <- unlist(b1)#把list解开，就可以求和求均值了
c <- mapply(mean,iris[,1:3])#先喂函数，再喂数据
sum(mapply(function(x)x^2,iris[,1]))#mapply出来的类型不一定
mapply(rep,times=1:4,x=4:1)#4重复1遍，3重复2遍。。。
Lst <- list(name="Wxr",husband="Ldm",lunch="中餐",food=c("beef","tomato"))
#构造一个list
Lst$name#查看其中一项
x <- c(1,3,4)
y <- matrix(1:24,nrow = 3,byrow=T)
z <- list(x,y,c("a","b"))#把xyc合成一个list
z[[1]]#查看list第一项
x <- c(1,3,4)
y <- matrix(1:24,nrow = 3,byrow=T)
z <- c("a","b","c","d")
L <- list(L1=x,L2=y,L3=z)
L$L1#查看第一项
Lst_iris <- list()#先定义一个空的list或者空的vector
Vector_iris <- c()
for(i in 1:50){
Lst_iris[[i]] <- apply(iris[i,1:4],1,mean)
Vector_iris[i] <- sd(iris[i,1:4])
}#对空变量的每一项参与for循环
Lst_iris
Vector_iris
Lst_iris[[1]]
c(Lst_iris,Vector_iris)#将两个list合并
data.frame#数据框 eg iris，比list方正，里面的向量必须长度一致，矩阵有相同行数
#数据框是一种矩阵形式的数据，数据框中的各列可以是不同类型的数据，每列是一个变量
#，每行是一个观测
ID <- c(1,2,3,4)
name <- c("A","B","C","D")
score <- c(60,70,80,90)
student1 <- data.frame(ID,name)#构造两个小数据框
student2 <- data.frame(ID,score)
total_student1 <- merge(student1,student2,by="ID")
#因为student1，2有相同的ID，就可以用merge按照ID合成一个大的数据框
ID <- c(1,2,3,4)
name <- c("A","B","C","D")
score <- c(60,70,80,90)
sex <- c("M","F","M","M")
student1 <- data.frame(ID,name)
student2 <- data.frame(score,sex)
total_student2 <- cbind(student1,student2)
#因为student1，2没有相同的变量，可以用cbind进行列合并成一个大的数据框
ID <- c(1,2,3,4)
name <- c("A","B","C","D")
student1 <- data.frame(ID,name)
ID <- c(5,6,7,8)
name <- c("E","F","G","H")
student2 <- data.frame(ID,name)
total_student3 <- rbind(student1,student2)
#因为student1，2没有相同的变量，可以用rbind进行列合并成一个大的数据框
total_student3
attach(iris)#把数据框引到界面里，有了此命令就可以直接用Sepal.Length
#可以attach多个数据在界面里
Sepal.Length
detach(iris)#把数据框从界面拿掉，后面就不能直接用Sepal.Length查找
Sepal.Length
mydataframe <- data.frame(
name=c("张三", "李四", "王五", "赵六", "丁一"),
sex=c("F","F","M","M", "M"),
age=c(16, 17, 18, 16, 19),
height=c(167.5, 156.3, 177.3, 167.5, 170.0),
weight=c(55.0, 60.0, 63.0, 53.0, 69.5) )
#可以在括号里写内容，直接列合并生成数据框
mydataframe[2:4,]#可以像矩阵一样索引
mydataframe[["name"]]#数据框独特的索引方式
mydataframe$name# $列索引
colnames(mydataframe) <- c("names","sexs","age","height","weight")#改列名字
row.names(mydataframe) <- c("第一行","第二行","第三行","第四行","第五行")
mydataframe1 <- t(mydataframe)#转至后就不是数据框了
mydataframe2 <- as.data.frame(mydataframe1)#强制转化成数据框
mydataframe2$第一行
mydataframe$H_W_ratio <- mydataframe$height/mydataframe$weight
#添加一列H_W_ratio，内容是身高与体重的比值
mydataframe$BIM <- mydataframe$weight/mydataframe$height^2
edit(mydataframe)#弹出一个表格，数据编辑器，可以直接在里面改（mac不行）
#但改完只能显示在屏幕里，不能保存在data.frame里
mydataframe <- edit(mydataframe)#给他赋值以后才能保存下来，赋值成新的mydataframe
fix(mydataframe)#跟edit一样，但是可以直接存到data.frame里
subset(mydataframe,name=="王五")#只调王五的信息出来
subset(mydataframe,name=="王五"|name=="赵六")#调王五跟赵六的信息
subset(iris,Species=="setosa")#注意是双等号，就是判断，是setosa的就调出来
library("MASS")#读一个以前下载的包
Cars93 <- Cars93#看cars93调出来
USA_pro <- subset(Cars93,Cars93$Origin=="USA")#单条件选取
USA_pro <- subset(Cars93,Cars93$Origin=="USA",select=c(Type,Price))
#挑选出满足条件的type跟price，只看他俩
pro <- subset(Cars93,select=c(Type,Price))#只看所有的type跟price
USA_pro_mode_Astro <- subset(Cars93,Cars93$Origin=="USA"&Cars93$Model=="Astro")
#多条件选取
pie(a,b)
a <- c(40,50,60,70)
b <- c("a","b","c","d")
pie(a,b)
pie(a,b,col = rainbow(4))
levels(Cars93$Manufacturer)#看manufacturer有几种，后面好挑选画图
#挑选"Acura" "Audi" "BMW" "Buick"
M <- subset(Cars93,Manufacturer=="Acura")
dim(M)#取M的维度，是个向量
dim(M)[1]#可以只看Acura行的维度，就知道种类是Acura的有几个了
length(which(Cars93$Type=="Compact"))#用这个也可以知道类型是Compact的有几个
a <- c(dim(subset(Cars93,Manufacturer=="Acura"))[1],
dim(subset(Cars93,Manufacturer=="Acudi"))[1],
dim(subset(Cars93,Manufacturer=="BMW"))[1],
dim(subset(Cars93,Manufacturer=="Buick"))[1])
#取"Acura","Audi","BMW","Buick"行的维度数来画图
table(Cars93$type)#也可以直接用table来代替上面的步骤，直接统计出每种类型的数量
b <- c("Acura","Audi","BMW","Buick")#给每一部分起名字
pie(a,b,col=c("skyblue","ligthgreen","red","yellow"))
pairs(iris)#跟plot(iris)一样#当两个数据相关性很大时，可以去掉一个只留一个
#scatterplot散点图
dev.new()
coplot(Cars93$Fuel.tank.capacity~Cars93$Price|Cars93$Origin)
#在origin的条件下，看油箱存油量跟价格的关系
#coplot（因变量，自变量|因子变量）按照因子变量来画因变量自变量的联合散点图
coplot(Cars93$Fuel.tank.capacity~Cars93$Price|Cars93$Origin+Cars93$DriveTrain)
#看油箱存油量跟价格在不同产地下的关系，油箱存油量跟价格在不同驱动下的关系
hist(iris$Sepal.Length)#画频数分布直方图
hist(iris$Sepal.Length,nclass=4)#nclass规定分成几组，就把前面那个图的条合并一下
#但有的时候nclass不一定管用，因为他是平均分配，把上面那个nclass改成5就不行
hist(iris$Sepal.Length,breaks=c(4.3,4.5,5.5,6.5,7.5,7.9))#自定义分界点、组数
hist(iris$Sepal.Lengt,probability = T)#纵坐标变成频数的占比
hist(iris$Sepal.Length,probability = F)#纵坐标是频数
hist(iris$Sepal.Length,probability = T,main = "Length")#main给图起个标题，可中文
par(mfrow = c(1,2))#把画图窗口分成一行两列的块（左右两块）
data <- c(rep(1,10), rep(2,5), rep(3,6))
data #rep(2,5)产生5个2
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), probability = T, main = "A")
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), probability = F, main = "B")
par(mfrow = c(1, 2))
data <- c(rep(1, 10), rep(2, 5), rep(3, 6))
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), axes = T, main = "axes = T")
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), axes = F, main = "axes = F")
#axes=F就是不要坐标轴
par(mfrow = c(1, 2))
data <- c(rep(1, 10), rep(2, 5), rep(3, 6))
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), col = "pink")#给颜色
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), col = rainbow(3))
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), col = rainbow(3),border=NA)
#去掉每一个框线
par(mfrow = c(1, 3))
data <- c(rep(1, 10), rep(2, 5), rep(3, 6))
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), density = 1, main = "density = 1")
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), density = 2, main = "density = 2")
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), density = 20, main = "density = 3")
#给条添点线，密度大就是添的线多，没具体值
par(mfrow = c(1, 3))
data <- c(rep(1, 10), rep(2, 5), rep(3, 6))
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), density=2,angle=45)
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), density=2,angle=90)
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), density=8,angle=30)#添的线的角度30度
hist(data, breaks = c(0.5, 1.5, 2.5, 3.5), density = 5,angle = 60,col = "pink")
barplot(GNP ~ Year, data = longley)#画柱状图
barplot(cbind(Employed, Unemployed) ~ Year, data = longley)#做堆叠条形图，即因变量可以多个
data("longley")#longley是软件自带的数据，跟iris一样
barplot(cbind(Employed, Unemployed,Population) ~ Year, data = longley)
barplot(cbind(Employed, Unemployed) ~ Year, data = longley,col=rainbow(2))
boxplot(iris$Sepal.Length~iris$Species)#把因子变量放后面，箱线图，超出横线是异常数据
boxplot(len ~ dose:supp, data = ToothGrowth,
boxwex = 0.5, col = c("orange", "yellow"),
main = "Guinea Pigs' Tooth Growth",
xlab = "Vitamin C dose mg", ylab = "tooth length",
sep = ":", lex.order = TRUE, ylim = c(0, 35), yaxs = "i")
ToothGrowth <- ToothGrowth
#len因变量（纵坐标），dose跟supp两个因子变量两两搭配做自变量
which(iris$Sepal.Length<5&iris$Species=='virginica')#查找满足条件的数据，会显示他所在的行
#用这个办法可以看箱线图的异常数据
#boxwex箱子的宽度 main给图起名字 xlab横坐标的名字 ylab纵坐标的名字
#sep 两个因子变量之间的连接符号 ylim纵坐标的范围 lex.order自变量的排列顺序
# yaxs='i'表示刻度线都在数据范围内部，紧贴着ylim的上下限
boxplot(len ~ dose:supp, data = ToothGrowth,
boxwex = 0.8, col = c("orange", "yellow","green"),
main = "boxplot",
xlab = "x axes", ylab = "y axes",
sep = "+", lex.order = FALSE, ylim = c(0, 50),yaxs = "r")
dotchart(cars$speed)#可以看到同样的值有几个
dotchart(cars$speed,xlim=c(10,20))#xlim规定横坐标的范围
length(which(Cars93$Type=="Compact"))
x <- seq(-10,10,0.1)
y <- x
z <- outer(x,y,function(a,b)a^2+b^2)#外积
class(z)
image(x,y,z)#画一个热力图,没有数据
contour(x,y,z)#画一个等高线图,因为z是x、y的平方和是个圆圈
z <- outer(x,y,function(a,b)a^2-b^2)
image(x,y,z)
contour(x,y,z)#画出来是个马鞍面的
persp(x,y,z)#画一个三维立体图，也是个马鞍面
#contour图是persp同一圈的高度向底面映射实现的（即xoy面的投影）
dev.new()
persp(x,y,z,col= heat.colors(30),theta = 45,phi=15,r=sqrt(3),d=1,axes = FALSE)
#theta 设置视角的方位角方向，phi为设置视角的余维度. r 观察点到图中心的距离.
#d 数值，用于增强或减弱透视效果. scale 设置在画图时是否要保持高度比例.
#expand 用于扩大或缩小z坐标.col 曲面表面的颜色. border 曲面边框的颜色.
#ltheta, lphi 如果设置了这个值，曲面的光源就是根据ltheta和lphi设置的角度来绘制.
#shade 计算曲面阴影的参数. axes 表示是否要画坐标轴.
#ticktype设置标记的类型.
require(grDevices) # the same to library(grDevices)
x <- seq(-10, 10, length= 30)
y <- x
f <- function(x, y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
z <- outer(x, y, f)
z[is.na(z)] <- 1
op <- par(bg = "pink")#bg添加背景颜色
persp(x, y, z, theta = 30, phi = 30,
expand = 0.5, col = rainbow(1000))
persp(x, y, z, theta = 45, phi = 20,
expand = 0.5, col = rainbow(1000),
r=180,
ltheta = 120,
shade = 0.75,
ticktype = "detailed",
xlab = "X", ylab = "Y", zlab = "Sinc( r )" ,
border=30,
axes = F)
x <- seq(-pi, pi, len=50)
y <- x
f <- outer(x, y, function(x, y) cos(y)/(1 + x^2))
contour(x,y,f)
persp(x,y,f,col=heat.colors(30))
contour(x,y,f)
contour(x,y,f,nlevels=3,col='red',add=TRUE)#nleves=3增加3条等高线
#add=TRUE 在上面画的那个图上增加一些细节,而不是画新的，一般用于两个高级做图之间
image(x,y,f)
contour(x,y,f,add=TRUE)
plot(1:5, 1:5, type = "n", xlim = c(0, 6), ylim = c(-1, 8))#先画了一个空框
symbols(x = 1:5, y = 1:5, circles = rep(1, 5), inches = FALSE, add = TRUE)
#add=T，就是把add前面的命令添加在上一个画的图里，上面两个命令相当于下面这一个
symbols(x = 1:5, y = 1:5, circles = rep(1, 5),inches = FALSE,
xlim = c(0, 6), ylim = c(-1, 8))
inches=T #英寸，画的圆大一点，
#x = 1:5, y = 1:5 目的是给圆配圆心(1,1)(2,2)...circles = rep(1, 5)画半径是1的圆五个
plot(1:5, 1:5, type = "n")#若没有type=‘n’，就会画成散点图，type=‘n’就是啥也不画
plot(1:5, 1:5)#x=1,2,3,4,5 y=1,2,3,4,5 搭配成坐标(1,1),(2,2)...画成散点图
contour(x,y,f,axes = F)#不要坐标轴
contour(x,y,f,axes = T)#要坐标轴
plot(1:100)
plot(1:100,log="x")#设置x轴成为对数轴
#相当于下面
x <- 1:100
y <- log(x)
plot(1:100,log="y")#设置y轴成为对数轴
plot(1:100,log="xy")
x <- seq(-10,10,0.5)
y=sin(x)
plot(x,y,type='p')#散点图
plot(x,y,type='l')#折线图
plot(x,y,type='b')#不穿过点的曲线图
plot(x,y,type='o')#穿过点的曲线图
plot(x,y,type='h')#连垂直线
plot(x,y,type='s')#连成阶梯形，一个点连另一个点的时候先横向走
plot(x,y,type='S')#连成阶梯形，一个点连另一个点的时候先竖向走
plot(x,y,type='o',col='red',font.main = 3,main="y=sinx",xlab="aaa",ylab="bbb",font.lab=2)
#font.main改标题的字体 font.lab改变坐标轴名称的字体 font.axis改变坐标轴的字体
#1：normal 2：斜体italic 3：粗体bold 4：粗斜体 bold italic
plot(x,y,type='o',col='red',font.axis=2)
plot(x,y,type='o',col='red',font.main = 3,main="y=sinx",sub="sin function")
#sub=''添加副标题
x <- 1:30
y <- sin(x)
plot(x,y)#先画一个散点图
lines(x,y,col="red")
points(x,y,type = 's')
#lines,points都是低级做图命令，它是在高级做图命令下(eg:plot)，往原来的图里添加一些别的信息
text(20,0.5,"lalala",col='blue')#在x=24，y=0.5的位置写一点文字
mtext('I want to China',col='green',side=4)
#在图形边框外写文字，side=1下面，side=2左面，side=3上面,side=4右面
setwd('//Users//wangxinran//Desktop//R语言')
HW <- read.table('HW.txt',header = T)
plot(HW)
x <- HW[,1]
y <- HW[,2]
plot(x,y)
lm.out <- lm(HW$Height~HW$Weight)#lm：找二者之间回归直线y=ax+b中的参数
lm.out#intercept对应b（截距），0.2719是a
abline(152.1136,0.2719,col='red' )#把回归直线直接画到上面的散点图里，注意别喂错数据位置
abline(lm.out$coefficients,col='blue')#第二种画法，直接取系数，就不会写错数了
abline(lm.out,col='green')#第三种画法
#做回归：跟数据，我们估计他们的关系是线性关系，用lm找到线性回归的系数，再用abline画出直线
#用poygon进行纯色填充
xx<-c(1:100,100:1)
yy<-rev(cumsum(xx))#cumsum是累计求和，rev是颠倒数字顺序
x <- c(1,2,3,4)
y1 <- cumsum(x)
y2 <- rev(y1)
plot(xx,yy,type='l')
plot(xx,yy,type='n',xlab="Time",ylab="Distence")#先画一个空白框，再用polygon给封闭位置填充
polygon(xx,yy,col="gray",border="red")#border边界颜色，col给封闭图形填充颜色
title("I am happy to do this")#给上面的图加标题
plot(c(1, 9), 1:2, type = "n")
polygon(1:9, c(2,1,2,1,NA,2,1,2,1), density = c(10, 20), angle = c(-45, 45))#填充线
#NA就是4～6中间的空白，就画成了两个封闭图形
#legend给图例
x <- seq(-pi, pi, len=65) #from –pi to pi, produce 65 data
plot(x, sin(x), type= "l",ylim = c(-1.2, 1.8), col=3,lty =2)#lty线型，lty=2是虚线
#col取3就代表绿色，若他换成green，legend里也要换成green，必须一一对应
points(x,cos(x),pch=3,col=4)#pch表示数据点位置用加号标出，不同的数值代表不同的加号
lines(x,tan(x),type="b",lty=1,pch=4,col=6)#type=b表示交替出现
title("legend(..., lty = c(2, -1, 1), pch = c(NA, 3, 4),
merge = TRUE)", cex.main = 1.1)#cex.main标题字体大小
legend(-1,1.9,c("sin","cos","tan"),col=c(3,4,6),text.col= "green4",
lty= c(2,-1,1),pch=c(NA,3,4), merge =TRUE, bg="gray90")
#-1，1.9 表示图例放的位置，text.col是图例里面字体的颜色，bg是图例背景颜色，merge=T居中对齐
#注意就是每一项的信息跟上面的函数信息一一对应，没有的位置就用NA
#不同数字代表不同类型eg不同颜色，不同线型
y <- function(x)log(x)+sqrt(x)+x^(1/3)
plot(y,1,1000,main=expression(y=log(x)+sqrt(x)+sqrt(x,3)),lwd=3,col="blue")
#用expression在标题中写数学表达式（1）
text(600,20,expression(paste(bgroup("(",atop(n,x),")"), p^x, q^{n-x})))
#用expression在图里面600，20的位置标注数学表达式（2）
#用paste把bgroup("(",atop(n,x),")")，p^x, q^{n-x}三个粘贴到一起
text(800,30,expression(paste(bgroup("{",atop(n,k),"}"), p^k, q^{n-k})))
title(expression(x%in%A))#体会不同数学表达式的写法x%in%A x属于A
#也可以不用main，用title加标题，配合expression写数学表达式（3）
par(family="Sarasa Gothic CL")#这个命令运行后就可以使用中文字体了
#交互式命令 locator，identify
x = rnorm(10)#产生十个随机数
plot(x)
locator(5,"o")
text(locator(1), "Outlier", adj=0)#把文字写到你点的地方
#当你画了一个散点图，你想把其中几个点连成折线，就用locator命令，数字表示你想点几个点
#这样你点到的那些点他就会给你连起来，也可以不是你画出来散点图中的点，图里所有的位置都可以
require(graphics)
hca<- hclust(dist(USArrests))
plot(hca)
(x <- identify(hca))#identify分层聚类，点上面的分支竖线，他就会把你点的那一类分支包括起来
#当你执行了交互式命令时，运行窗口有一个红色的圈写着stop，当你点完想点的，就点stop结束交互
x<-1:100
y<-sin(x)
plot(x,y,type = "l")
identify(40,1,"标一下") #你需要在坐标位置[40,1]那里点击一下，添加文字，点到别处则会出错
#par命令
#一、关于颜色
opar<-par(no.readonly=TRUE) #have a try to write par(no.readonly=TRUE)
#把默认的参数都存下来，存到opar里，最后再恢复它，以免画后面的图受到用par设置的参数的影响
x<-seq(1,10,length.out=100)-5#每个数据都要剪掉5
y<-c(log(x[x>0]),log(abs(x[x<=0])))
par(fg="red")#前景颜色
par(bg="yellow")#背景颜色
par(col='purple')#里面画的线的颜色
par(col.axis='green',col.lab="white")#坐标刻度值的颜色，坐标轴名字的颜色
plot(x,y, type="l")
title(main="BlackbgAndRedfg",col.main="pink",sub="Byprogram-dog.blogspot.com",
col.sub="blue")
par(opar)# par恢复默认参数，这样你后面画的图还是默认的设置（如默认颜色等）
#也可以不用opar，直接把昨天窗口扫空也行
par(fg="red",bg="yellow",col='purple',main="BlackbgAndRedfg",col.main="pink",
sub="Byprogram-dog.blogspot.com",col.sub="blue")#可以把上面的都合并到par里
#注意要把par放在高级做图命令plot之前，才能显示出par里的条件信息
#二、关于字体
#1】font 可以分别设置部分字体 粗体斜体粗斜体
opar<-par(no.readonly=TRUE)
x<-seq(-10,10,length.out=100)
y<-sin(x)
par(font.axis=1)# 1 normal
par(font.lab=2)# 2 bold
par(font.main=3)# 3 italic (type)
par(font.sub=4)# 4 bold Italic
plot(x,y,type='l')
title(main="fontstyle",xlab = "粗体", sub="Byprogram-dog.blogspot.com")
par(opar)
#2】family 统一设置全部字体 宋体，黑体，楷体 serif", "sans" and "mono"
dev.new()
opar<-par(no.readonly=TRUE)
par(mfrow=c(4,1))
x <- -10:10
y <- -(x^2)
par(family="mono") #mono字体
plot(x,y,type='l')
title(main="family mono style",sub="Byprogram-dog.blogspot.com")
par(family="") #默认字体default
plot(x,y,type='l')
title(main="family default style",sub="Byprogram-dog.blogspot.com")
par(family="serif") #serif字体
plot(x,y,type='l')
title(main="family serif style",sub="Byprogram-dog.blogspot.com")
par(family="sans") #sans字体
plot(x,y,type='l')
title(main="family sans style",sub="Byprogram-dog.blogspot.com")
par(opar)# Restoring the parameters of par
#三、关于字号
#1】ps直接设置磅值（字号大小）,只能设置全部字体
opar<-par(no.readonly=TRUE)
x <- seq(-10,10,length.out=100)
y <- log(x^2)
dev.new( )
par(mfrow=c(3,1))
par(ps=10)#10
plot(x,y,type='l')
title(main="fontsize: ps=10",sub="Byprogram-dog.blogspot.com")
par(ps=15)#15
plot(x,y,type='l')
title(main="fontsize: ps=15",sub="Byprogram-dog.blogspot.com")
par(ps=20)#20
plot(x,y,type='l')
title(main="fontsize:ps=20",sub="Byprogram-dog.blogspot.com")
par(opar)
#2】cex取normal情况下的倍数 par(cex.main=1.5)，可以只改部分字体
opar<-par(no.readonly=TRUE)
x<-seq(-10,10,length.out=100)
y<-sin(log(x^2))
dev.new()
par(mfrow=c(2,1))
par(cex.main=1)
plot(x,y,type='l')
title(main="fontsize:cex.main=1",sub="Byprogram-dog.blogspot.com")
par(cex.main=1.5)
plot(x,y,type='l')
title(main="fontsize:cex.main=1.5",sub="Byprogram-dog.blogspot.com")
par(opar)
#cex也可以取ps的倍数 par(ps=20,cex.main=0.5)
opar<-par(no.readonly=TRUE)
x<-seq(-10,10,length.out=100)
y<-sin(log(x^2))
dev.new()
par(mfrow=c(2,1))
par(ps=20,cex.main=0.5)
plot(x,y,type='l')
title(main="fontsize:cex.main=1",sub="Byprogram-dog.blogspot.com")
par(ps=20,cex.main=1.5)
plot(x,y,type='l')
title(main="fontsize:cex.main=1.5",sub="Byprogram-dog.blogspot.com")
par(opar)
#cex.axis for axis 坐标刻度值的字号
#cex.lab for labels坐标轴名字的字号
#cex.main for title主标题的字号
#cex.sub for sub-title副标题的字号
#四、关于线型 lty
#Style of line，1: full line；实线2: dashed line；虚线3: dotted line；
#点化线4: dot-dashed line；5: long dashed line
opar<-par(no.readonly=TRUE)
x<-seq(-10,10,length.out=100)
y<-sin(log(x^2))
par(lty=1)
plot(x,-y,type='l',col="red",ylim=c(-3,3))
par(lty=2)
lines(x,y,type='l',col="blue")
title(main="lty",sub="Byprogram-dog.blogspot.com")
par(opar)
#五、关于标识符,一个数据点给一个标识符
#pch 0～25种
x<-seq(-10,10,length.out=20)
y1<-0.1*x^2
y2<-0.2*x^2
y3<-0.4*x^2
y4<-0.8*x^2
y5<-1.6*x^2
y6<-3.2*x^2
par(pch=1)#pch1 圆圈
plot(x,y1,type='b',col="red",xlim=c(0,5))
par(pch=2)#pch2 三角
lines(x,y2,type='b',col="blue")
par(pch=3)#pch3 加号
lines(x,y3,type='b',col="green")
par(pch=4,lty=3)#pch4 乘号
lines(x,y4,type='b',col="red")
par(pch=0)#pch5 正方形
lines(x,y5,type='b',col="blue")
par(pch=6)#pch6 倒三角
lines(x,y6,type='b',col="green")
title(main="pch",sub="Byprogram-dog.blogspot.com")
par(opar)
plot(x,y1,type='b',col="red",xlim=c(0,5),pch=5)#也可以直接把pch放plot里
#六、线宽 #lwd 线宽 cex 标识符宽
opar<-par(no.readonly=TRUE)
x<-seq(1,10,length=20)
y<-1/x
dev.new()
par(mfrow=c(2,1))
plot(x,y,type="b",pch=2,cex=5,lty=3,lwd=1)
title(main="lwd=2 and cex=2",sub="Byprogram-dog.blogspot.com", cex.sub=0.5)
plot(x,y,type="b",pch=2,cex=1,lty=3,lwd=5)
title(main="lwd=1andcex=1",sub="Byprogram-dog.blogspot.com",cex.sub=0.5)
par(opar)
#
x<-seq(1,10,length.out = 20)
y<-1/x
dev.new()
par(pch=2,cex=2,lty=3,lwd=2,mfrow=c(2,1))#1
plot(x,y,type="b")
par(pch=2,cex=0.5,lty=3,lwd=1)#2
plot(x,y,type="b")
#把命令从plot里拿出来，要注意1.mfrow命令要放在第一个par命令的后面
#2.cex放在par里面就是改了全部的字体大小，而不是只改了标识符宽
#一定要明确par参量的意义，即在恢复默认设置之前，画图用到的参量都以par设置的为准
#若plot里也有对par里面设置出的参量做重新设置，那么以后设置的为准
#七、分界面
par(mfcol=c(3,2))#分作图界面，先按列走
par(mfrow=c(3,2))#分作图界面，先按行走
#建模1】最近邻居算法 K-NN
#根据离A最近的邻居的种类取类比推断出A的类别，重点在如何判断“最近”
#先把数据的指标画到坐标图中，再找到图中离待判断值比较近的点，作为序列样本
#再计算待判断值跟图中各点的距离，按大小排列好顺序标号，选取最近邻居个数k，
#看k里最多的种类是什么，用此来估计待测值的类别，注意k取的不同结果可能不同
#k的选取应该在两个极端值（1 跟 全部数据）之间，一般是图中离待判断值比较近的点的个数
tomato <- c(6,4)
grape <- c(8,5)
orange <- c(7,3)
sqrt(sum((tomato-grape)^2))
sqrt(sum((tomato-orange)^2))#计算两点之间的距离
#数据预处理
#1 数据的范围压缩（归一化）：把每一列数据等比例压缩至同一范围（0～1），不改变相对大小关系
#方法1，同除以每列的最大值
#方法2，Xnew=[X-min(X)]/[max(X)-min(X)]
#方法3，压缩至（-1，1）之间，用概率 Xnew=[X-Mean(X)]/StdDev(X)(标准差)
#不预处理，一些小的数据对距离的影响就会很小，量纲大的占优势，这样相对位置就不太准确
#2 把字符变量变成数值变量，便于距离计算 eg：male=1；female=0
iris <iris
sqrt(sum((iris[1,1:4]-iris[150,1:4])^2))#依据前四列数据，计算第1朵花跟150朵花的距离
#一、读数据
setwd('//Users//wangxinran//Desktop//R语言')
#判断一个病人是否得病，依据knn算法，将该病人的指标与数据库内的数据比较，由最近邻居得出结论
wbcd <- read.csv('wisc_bc_data.csv',stringsAsFactors = F)#读入数据，把因子变量转换成字符型
#二、查看并分析数据，做简单处理
str(wbcd)#看各类数据的类型分别是什么，diagnosis这一列就显示成chr 而不是 factor
wbcd <- wbcd[,-1]#去掉第一列，即id
table(wbcd$diagnosis)#看看样本里良性恶性分别是多少，以此来看样本数据是否合理
#尽量让各类样本在数据库里的比重均衡
wbcd$diagnosis <- factor(wbcd$diagnosis,levels=c("B","M"),labels="Benign","Malignant")
#想把第一列的字符B、M写成它完整的单词命名，要先改回成因子型，在用labels改名
summary(wbcd[c("radius_mean","area_mean","smoothness_mean")])
#查看数值型变量的一些情况，1st Qu 四分之一分为数 3st Qu 四分之三分为数
#三、归一化
normalize <- function(x){
return((x-min(x))/(max(x)-min(x)))
}#自定义一个归一化函数,对数据做预处理，x代表一列向量，return表示要该函数的返回值
wbcd_n <- as.data.frame(lapply(wbcd[,2:31],normalize))#对数据的2～31列归一化
summary(wbcd_n$radius_mean)#查看一下归一化是否成功
#四、取数据
wbcd_train <- wbcd_n[1:469,]#取一部分数据做训练样本，注意训练样本都是数值型的，也叫
#训练特征数据，所以应在归一化后的wbcd_n里取
wbcd_train_labels <- wbcd[1:469,1]#取训练标签（对应训练样本是良性还是恶性），注意
#要去wbcd里取，因为归一化的数据框里已经没有diagnosis这一列了
wbcd_test <- wbcd_n[470:569,]#取测试样本，也叫测试特征数据
wbcd_test_labels <- wbcd[470:569,1]#取测试标签
#从920行开始，如果怕弄错到底是从wbcd里取，还是从wbcd_n里取，可以先把标签跟归一化的数据
#合并成wbcd_n_1，再统一从取wbcd_n_1里取
wbcd_n_1 <- cbind(wbcd$diagnosis,wbcd_n)
wbcd_train_1 <- wbcd_n_1[1:469,2:31]
wbcd_train_1_labels <- wbcd_n_1[1:469,1]
wbcd_test_1 <- wbcd_n_1[470:569,2:31]
wbcd_test_1_labels <- wbcd_n_1[470:569,1]
#五、引入模型计算
library(class)#取class包是为了要用里面的knn模型
wbcd_test_pred<- knn(train = wbcd_train, test = wbcd_test,
cl = wbcd_train_labels, k = 21)
#直接用建模模型knn函数来计算测试样本的标签，看看跟我们选出来的真正的测试标签有多大差距
#若差距不大，说明我们取的k跟建立的knn模型比较合适，可以用来做良性恶性的判断
#一般先让k等于训练样本数开平方，注意knn里的各个元素别喂错了
#六、评价模型 1.table 2.CrossTable
#1.
table(wbcd_test_pred,wbcd_test_labels)#wbcd_test_pred模型得到的，wbcd_test_labels真实的
#评价模型，统计用模型测出来的良性恶性跟真实的良性恶性分别是多少，
#并且可以简单验证，统计出来主对角线上的结果是判断对了的个数，其他位置是判断错误的个数
#2.
install.packages("gmodels")#用这个包里的CrossTable来评价模型
library('gmodels')
CrossTable(wbcd_test_pred,wbcd_test_labels)#CrossTable叫混淆矩阵，把模型得到的跟真实的
#混在一起比较，横纵变量交叉点就是两者重合的，即正确的，其余是模型判断有误的
#表格最前面的Cell Contents会告诉你表格里每一个数据的意思
#Chi-square contribution卡方贡献率
CrossTable(wbcd_test_pred,wbcd_test_labels,chisq = F)
#七、模型的改进 1.更换归一化的方法 2.跟换k
#1.
wbcd_z<- as.data.frame(scale(wbcd[-1]))#scale压缩数值型数据，把每个数据减均值除以标准差
summary(wbcd_z$area_mean)
wbcd_train <- wbcd_z[1:469,]
wbcd_train_labels <- wbcd[1:469,1]
wbcd_test <- wbcd_z[470:569,]
wbcd_test_labels <- wbcd[470:569,1]
wbcd_test_pred<- knn(train = wbcd_train, test = wbcd_test,
cl = wbcd_train_labels, k = 21)
table(wbcd_test_labels,wbcd_test_pred)#发现错的更多了
#2
wbcd_test_pred<- knn(train = wbcd_train, test = wbcd_test,
cl = wbcd_train_labels, k = 15)
table(wbcd_test_labels,wbcd_test_pred)#测试不同的k看哪个好
#用knn实验iris数据
iris <- iris
#观察数据：由于iris数据的species是按顺序排列的，直接取前100个数据是不行的
normalize <- function(x){
return((x-min(x))/(max(x)-min(x)))
}
iris_n <- as.data.frame(lapply(iris[,1:4],normalize))
#第一种选数据的方法 打乱排列顺序,再随机选取
index <- sample(150,120)#随机产生1～150的行号,随机取里面的120个
iris_train <- iris_n[index,]
iris_train_labels <- iris[index,5]
iris_test <- iris_n[-index,]
iris_test_labels <- iris[-index,5]
library(class)
iris_test_pred<- knn(train = iris_train, test = iris_test,
cl = iris_train_labels, k = 11)
table(iris_test_labels,iris_test_pred)#产生的sample行号不同，结果不同
#k-折交叉验证 因为sample是随机产生的，所以每次测试样本的结果都不一样，为了判断模型是否
#有效，可以多测试几次不同的sample，取成功结果的比例的平均值
#第二种取数据的方法 每一类随机取一部分
table(iris$Species)#看每一类有几个
index1 <- sample(50,40)#随机产生1～50，50个行标随机取40个
index2 <- sample(51:100,40)#随机产生第51～100个行标随机取40个
index3 <- sample(101:150,40)#随机产生第101~150个行标随机取40个
iris_train <- iris_n[c(index1,index2,index3),]
iris_train_labels <- iris[c(index1,index2,index3),5]
iris_test <- iris_n[-c(index1,index2,index3),]
iris_test_labels <- iris[-c(index1,index2,index3),5]
library(class)
iris_test_pred<- knn(train = iris_train, test = iris_test,
cl = iris_train_labels, k = 11)
table(iris_test_labels,iris_test_pred)
#第三种取数据的方法 每一类顺序取前40个
iris_train <- iris_n[c(1:40,51:90,101:140),]
iris_train_labels <- iris[c(1:40,51:90,101:140),5]
iris_test <- iris_n[-c(1:40,51:90,101:140),]
iris_test_labels <- iris[-c(1:40,51:90,101:140),5]
library(class)
iris_test_pred<- knn(train = iris_train, test = iris_test,
cl = iris_train_labels, k = 11)
table(iris_test_labels,iris_test_pred)
#第四种取数据的方法 打乱排列顺序,再顺序选取
index<-sample(150)
iris<-iris[index,]
iris_n<-iris_n[index,]
iris_train<-iris_n[1:120,]
iris_train_labels<-iris[1:120,5]
iris_test<-iris_n[121:150,]
iris_test_labels<-iris[121:150,5]
iris_test_pred<- knn(train = iris_train, test = iris_test,
cl = iris_train_labels, k = 11)
table(iris_test_pred,iris_test_labels)
#发现：就算用样本数据去测试样本数据，结果也不一定全部成功
#第七章 读取及存储文件
#read.table一般用于读txt文件，其数据之间的分隔符号是空格形式，read.table读取空格分隔的数据
#read.csv一般用于读取csv文件，其数据之间的分隔符号是逗号形式，read.csv读取逗号分隔的数据
#当用read.table去读csv文件时，要用sep把分隔符改成逗号。read.csv则不用
#当用read.csv去读txt文件时，要用sep把分隔符改成空格。read.table则不用
#！！在读数据之前，要先看看数据之间的分隔符号，在读取时用sep把数据之间的格式修改
setwd('//Users//wangxinran//Desktop//R语言')
wine178 <- read.table(file='wine178.csv',header=T,sep=',')
wine178_1 <- read.csv(file='wine178.csv',header=T)
winequality_white <- read.csv(file='winequality-white.csv',sep=';')
winequality_white_1 <- read.table(file='winequality-white.csv',header=T,sep=';')
#read.table默认header=F，想要表头必须要加header=T，但read.csv默认header=T，可以不加
winequality_white_2 <- read.csv2(file='winequality-white.csv')
#read.csv2读取时默认数据之间的分隔符号是分号；，
HW <- read.table(file='HW.txt',header=T)
HW_1 <- read.csv(file='HW.txt',sep='',header=T)
usedcars <- read.table(file="usedcars.csv",sep=',',header=T)
usedcars_1 <- read.table("usedcars.csv")
str(usedcars)
mushrooms <- read.csv(file='mushrooms.csv')
mushrooms_1 <- read.table(file='mushrooms.csv',header=T,sep=',')
str(mushrooms)
#读取excel文件
install.packages('readxl')
library('readxl')
#安装'readxl'后用read_excel()命令可以读excel文件
#或者把.xls、.xlsx文件转化成.csv或.txt文件，再来读取，直接文件另存为就行了
temp_1 <- read_excel('temp_1.xlsx')
#也可以用如s下方法：File-》input dataset -》from excel-》把文件名跟位置拷贝过来
score <- sample(85:100,36,replace=T)#随机产生在85～100之间的数字36个
temp_1$score <- score#把score添加到temp_1表格中，新添一列
#readLines 读取纯文本
abstract <- readLines('abstract.txt')#读入文件，逐行读取文本
abstract
data<-readline()#从键盘中一行一行获取信息
#从键盘中输入想要的文字： I am a student
data
#用scan读数据
wine178_scan <- scan(file='wine178.csv',sep=',',skip=1)#skip=1向下一行再读，第一行不要
#scan()读取速度比较快，适用于大型数据，但第一行的字符型数据读不进去，且读进来不是数据框形式
#需要做如下改变，变成数据框并且加上表头第一行
wine178_Matrix <- matrix(wine178_scan,nrow=178,byrow=T)
wine178_dataframe <- as.data.frame(wine178_Matrix)
colnames(wine178_dataframe) <- c('Alcohol','Malic' ,'acid', 'Ash Alcalinity of ash' ,
'Total' ,'phenols', 'Flavanoids','Nonflavanoid','phenols'
,'Proanthocyanins', 'Color', 'intensitys' ,'Hue','wines')
usedcars_scan <- scan(file='usedcars.csv',what='character',sep=',',skip=1)
#如果数据是字符形式，要加what='character'，数据是numbeic不要紧
usedcars_Matrix <- matrix(usedcars_scan,nrow = 150,byrow = T)
usedcars_dataframe <- as.data.frame(usedcars_Matrix)
colnames(usedcars_dataframe) <- c('year','model','price','mileage','color','transmission')
#用write储存文件,将temp_1文件以新的文件名存储到指定目录下
write.table(temp_1,file="//Users//wangxinran//Desktop//R语言//temp_2.txt")
write.csv(temp_1,file="//Users//wangxinran//Desktop//R语言//temp_3.csv")
#也可以先将该目录改成当前目录，再储存时就不用指定了
setwd('//Users//wangxinran//Desktop//R语言')
write.table(temp_1,file="temp_2_2.txt")#write.table储存txt文件
write.csv(temp_1,file="temp_3_2.csv")#write.csv储存csv文件
write.csv(temp_1,file="temp_3_3.csv",row.names = F)#储存时不要行名字
write.table(temp_1,file="temp_2_3.txt",row.names = F,quote=F) #quote=F去掉双引号
iris <- iris
iris_sub <- iris[,1:4]
write.table(iris_sub,file='iris_sub.txt')
write.csv(iris_sub,file ='iris_sub.csv' )
#用cat存储一些新产生的数据到当前目录下
setwd('//Users//wangxinran//Desktop//R语言')
cat(1:10,file='temp11.txt')
cat(cbind(iris$Sepal.Length,iris$Petal.Length),file='iris_length.txt')
#用save存储.Rdata类型的文件，它是R语言下自己的文件格式，必须用R语言打开
setwd('//Users//wangxinran//Desktop//R语言')
save(iris_sub,file='iris_sub.Rdata')
save(list=ls(all=T),file="alldata.Rdata")#把environment里所有变量存盘
#给数据排序 sort rank order
S <- sample(40:55,8,replace=T)
sort(S)#把数据从小到大排序
sort(S,decreasing = T)#把数据从大到小排序
y <- c("wxr","bzx","skm","lxy")
sort(y)
#先按首字母排，如果首字母一样就比较第二个字母
order(S)#给出数据从小到大，在原始数据S中所在的位置
#给某个数据排序
library('readxl')
temp_1 <- read_excel('temp_1.xlsx')
score <- sample(85:100,36,replace=T)#随机产生在85～100之间的数字36个
temp_1$score <- score#把score添加到temp_1表格中，新添一列
order(temp_1$score)#显示第一个数据在原始数据中排第几，第二个数据在原始数据中排第几....
temp_2 <- temp_1[order(temp_1$score),]#给temp_1按成绩排序
S <- sample(40:55,8,replace=T)
S
rank(S)
#表示原始数据在排序后的序列里的位置，如果有相同数据排序，那么这个位置就取两个位置的平均值
#第一个数据在新数据中按从小到大顺序应该排第几,第二个...
#stack 自我补充，不作要求
#stack()函数：（堆栈的意思）
#长宽型数据的转换，长型：堆栈，数据间有不同的分类（如同属一类）；宽型：数据内容相对唯一
freshman <- c(12,23,24)
sophomores <- c(25,36,73)
juniors <- c(32,46,57)
data.frame(fr= freshman,so = sophomores,jun = juniors)
height <- stack(list(fresh =freshman,sopho = sophomores,juni = juniors))
height
tapply(height$values,height$ind,mean) #按照分类求均值
#练习，打开temp_1.xlsx文件，随机产生性别跟成绩，按成绩排序，并按性别求平均值
setwd('//Users//wangxinran//Desktop//R语言')
library('readxl')
temp_1 <- read_excel('temp_1.xlsx')
score <- sample(85:100,36,replace=T)
temp_1$score <- score
temp_1 <- temp_1[order(temp_1$score),]
sex <- sample(c("F","M"),36,replace=T)
temp_1$sex <- sex
tapply(temp_1$score,temp_1$sex,mean)
write.csv(temp_1,file="temp_1.csv")
#条件判断及循环
x<-3
y<-5
if(x<4)print("I want to go to school !")
if(x<1)print("I want to go to school !")
if(TRUE)print("I want to go to school !")
if(x<4|y>8)print("I want to go to school !")
if(x<4&y>8)print("I want to go to school !")
if(x>5)print('x>5')else print("x<=5")
if(x>5)print("I want to go to school !")else #注意else必须跟if在同一行
print("I stay at home")
if(x>4){
z = x+y;
print("z is :")
print(z)
}else
{
print("z is not exist")
}#条件满足的情况下，{}内所有的语句都会执行,条件不满足时执行else里所有的语句
#注意else必须紧跟在if结束的}位置
#if、else开始的{可以不紧跟在他们各自的后面，可以写在下一行
#如果if else在for循环里，else就可以不紧跟在if后面，可以写在下一行
teens <- read.csv("snsdata.csv")
teens <- ifelse(teens$age>15&teens$age<18,teens$age,NA)
#for循环
sum_x <- 0
x <- c(3,8,7,2,1,3,9,60,30,12)
for(i in x){
sum_x <- sum_x+i
}
sum_x
y <- 0
for(i in 1:dim(iris)[1]){
y <- y+iris$Sepal.Length[i]
}
y#求150个iris$Sepal.Length数据的和，dim(iris)[1]是iris的行数：150
y/dim(iris)[1]#求Sepal.Length的平均值
S <- 0#赋初值
for(n in 1:(10^5)){
S <- S+sqrt(3)/(2^n)
}
S
iris_na <- iris
for(i in 1:4){
iris_na[sample(1:nrow(iris),5,replace = F),i] <- NA
}
#在iris数据的第1～4列中，每一列里面都从150行里随机抽取5行(且不重复),把该位置数据变成NA
which(is.na(iris_na$Sepal.Length))#找到iris_na$Sepal.Length里NA所在的位置是第几行
for(i in 1:4){
iris_na[which(is.na(iris_na[,i])),i] <- mean(iris_na[,i],na.rm=T)
}
iris_na
# 把每一列里NA的位置，用这一列除了NA之外的数据所计算出来的平均值代替
mean(c(2,4,NA,6,8),na.rm=T)#na.rm=T去掉缺失值
#或者用以下方法，先求均值，再进行for循环
Mean_value <- sapply(iris_na[,1:4],mean,na.rm=T)
for(i in 1:4){
iris_na[which(is.na(iris_na[,i])),i] <- Mean_value[i]
}
iris_na
#while
S <- 0
n=1
while(n<=(10^5)){
S <- sqrt(3)/(2^n)
n=n+1
}
S
#repeat
f <- vector()
f[1] <- 1
f[2] <- 1
i=3
repeat{
f[i]=f[i-1]+f[i-2]
if((f[i])>=1000)break
i <- i+1
}
f
#switch 给它一个执行的位置比如1，它就执行第一个位置的操作
i=1
switch(i,mean(1:10),rnorm=10,print("I Love you"))#执行mean(1:10)求均值（第一项）
switch(i+1,mean(1:10),rnorm=10,print("I Love you"))#执行rnorm=10产生10个随机数（第二项）
switch(i+2,mean(1:10),rnorm=10,print("I Love you"))#执行print（第三行）
set.seed()#（）填初始种子点，随便填个自然数就行
#在执行rnorm或sample时先执行set.seed(自然数),就可以使每次产生的随机数都是一样的
#注意，每次都要运行一遍seed
set.seed(1)
rnorm(10)#rnorm产生服从标准正态分布的数
set.seed(113)
sample(85:100,10,replace = T)
#定义自己的函数
my_fun <- function(x,y) sin(x)+2*y
my_fun(3,8)
#做一个复杂计算，用大括号
twosam <- function(x1, x2) {
n1 <- length(x1); n2 <- length(x2)
xb1 <- mean(x1); xb2 <- mean(x2)
s1 <- var(x1); s2 <- var(x2)
s <- ((n1-1)*s1 + (n2-1)*s2)/(n1+n2-2)
tst<- (xb1 - xb2)/sqrt(s*(1/n1 + 1/n2))
tst #或者写成 return(tst), 还可以返回两个量，比如：return(list(tst, s))
}#return两个值的时候，必须把他们放到list里
x1<-c(1,2,3,3,4,5,6,7,8)
x2<-c(2,4,5,8,9,8)
a<-twosam(x1,x2)
a
#画分段函数图像
#method 1
My_f =function(x)
{
v1=2*x[x<=-1]+3
v2=((x+3)/(x+3))[(x>-1)&(x<=1)]#1要换成(x+3)/(x+3)，直接写1会报错
v3=x[x>1]^2
y=c(v1,v2,v3)
return(y)
}
#下面调用函数
x=c(-5:5)
y=My_f(x)
plot(x,y,type='l')
#method2
My_fun <- function(x){
if(x<=-1) y <- 2*x+3
else if((x>-1)&(x<=1)) y <- 1
else y <- x^2
return(y)
}
x <- seq(-3,3,0.1)
y <- vector()#y不能直接调用函数产生，需要用for循环把它的值放到一个vector里
for(i in 1:length(x)){
y[i] <- My_fun(x[i])
}#i指的是循环的次数，所以要注意(i in 1:length(x))而不是（i in x）
y
plot(x,y,type='l')
#定义一个二元运算
par(family="Sarasa Gothic CL")
"呵呵" <- function(x,y){cos(x)*sin(y)}
x <- c(1,2,3)
y <- c(4,5,6)
呵呵(x,y) #x"hehe"y不可用
"%!%" <- function(x, y) { cos(x)*sin(y)}
x<-c(1,2,3)
y<-c(4,5,6)
x%!%y #%!%(x,y)不可用
#函数的默认参数，如果不再赋值就用默认参数，后面如果有赋值，就用新的值
#缺省参数，把缺省参数"..."用大函数function里面的小函数比如summary来补上
my_func <- function(x, ...){
print(x)
summary(...)
}
my_func("This is iris data:", iris)
#函数里面的变量是局部变量，对全局没有影响
#线性回归
setwd('//Users//wangxinran//Desktop//R语言')
insurance.df <- read.table(file='Q1_Policies.txt',header=T)
#第一题，求回归系数 0.5134
attach(insurance.df)
#attch之后，在引用该数据中的某个变量就不用再写$了，全部已用完后用detach释放掉
insurance.lm <- lm(Policies~Quotes)#做回归，看回归系数
insurance.lm <- lm(Policies~Quotes,data=insurance.df)#当attach了不止一个数据时，要用data指明
insurance.lm
summary(insurance.lm)#可以看到更多的回归的结果
#Residual standard error 残差的标准误差（标准化残差）：3.451 拟合数据的自由度：58
#R-squared决定系数，越接近1说明模型拟合越好，一般一元回归看Multiple，二元回归看Adjusted
#*星号表示其对应量的显著程度
#第二题， 求回归系数的95%置信区间 [0.5134-2.001717*0.03308,0.5134+2.001717*0.03308]
qt(0.975,58)#求自由度为58的t分布的0.975分位数，=2.001717
#0.03308是标准误差Std. Error，可以从summary中得到
#第三题 用三种方法 P值、F分布、t分布 判断上述求得的回归系数在置信水平为95%时是否有效
#可以从summaty中得到t统计量跟F统计量 t value： 15.523 F-statistic：241
#也可以得到P值：p-value：< 2.2e-16,由于P<1-置信水平(0.05),所以我们算得的回归系数有效
#我们可以通过计算F(置信水平)(p,n-p-1) 跟上述得到的F比较，看该回归系数是否有效
qf(0.95,1,58)#求自由度是95%的F分布的(1,58)分位数
#算得qf(0.95,1,58)=4.006873<241,说明自变量对因变量有重要影响
#我们也可以通过计算t(1-置信水平) 跟上述得到的t比较，看该回归系数是否有效
qt(0.05,58)
#算得qt(0.05,58)=-1.671553，其绝对值< 15.523,说明自变量对因变量有重要影响
#第四题 用该回归去预测自变量是50时，因变量预测值的期望(E(y估计))的95%置信区间
predict(insurance.lm,data.frame(Quotes=50),se.fit=T,interval='confidence',level=0.95)
#predict(lm(y~x)，new，interval=“prediction”，level=0.95)
#——预测，lm(y~x)为之前通过训练数据拟合的回归方程；new为待预测的输入数据，其类型必须为数据框
#interval='confidence'表明想要y的期望的区间估计，interval='prediction' 表明想要y的区间估计
#level是置信水平 se.fit=T表示需要标准误差，=F则不需要
#结果中：fit表示自变量是50时因变量预测值(y估计)，lwr跟upr是你想求区间的上下限
#se.fit是E(y估计)公式中的最后一部分 df是拟合的自由度 residual.scale是残差的标准误差Sxx
#E(y估计)的置信区间=[y估计+t0.975(58)*se.fit,y估计-t0.975(58)*se.fit]
#第五题 用该回归去预测自变量是50时，因变量预测值(y估计)的95%置信区间
predict(insurance.lm,data.frame(Quotes=50),se.fit=T,interval='prediction',level=0.95)
#y估计的置信区间=[y估计+t0.975(58)*se.pred,y估计-t0.975(58)*se.pred]
#其中(se.pred)^2 = (se.fit)^2+(residual.scale)^2
#第六题 用该回归去预测自变量是40时，因变量预测值(y估计)的95%置信区间
#还可以把五、六题一起计算
predict(insurance.lm,data.frame(Quotes=c(50,40)),se.fit=T,interval='prediction',level=0.95)
plot(insurance.lm)#画图
#第一个图 残差大的数据在图中被标注了出来，残差点最好是随机分布的
residuals(insurance.lm)#计算残差
rstudent(insurance.lm)#计算删除学生残差
which(abs(rstudent(insurance.lm))>=3)#看删除学生残差中有没有大于3的，abs取绝对值
#第二个图 qq图，如果数据点基本在虚线附近（紧紧围绕直线），说明样本残差符合正态分布
#第三个图 标准化残差方根散点图，被标出来的点说明残差方跟较大，对回归影响不好，可以剔除
#第四个图 残差与杠图，横坐标是杠杆值，纵坐标是残差值，可以看两者的范围
#当杠杆值超过两倍或三倍的h杠，就说有高杠杆值点，但超过的这些点是否是异常点还要结合cook距离
cooks.distance(insurance.lm)#看各样本点的cook距离
hatvalues(insurance.lm)#看各样本点的杠杆值
plot(insurance.lm,which=c(1:6))#可以画出更多的图，比如cook距离
y=c(144,215,138,145,162,142,170,124,158,154,162,150,140,110,128,130,
135,114,116,124,136,142,120,120,160,158,144,130,125,175)
x=c(39,47,45,47,65,46,67,42,67,56,64,56,59,34,42,48,45,18,20
,19,36,50,39,21,44,53,63,29,25,69)
plot(y~x,type = 'p')
xueya_nianling.lm <- lm(y~x)
summary(xueya_nianling.lm)
plot(xueya_nianling.lm)
rstudent(xueya_nianling.lm)
which(abs(rstudent(xueya_nianling.lm))>=3)
cooks.distance(xueya_nianling.lm)
which(abs(cooks.distance(xueya_nianling.lm))>1)
h <-((1+1)/30)
2*h
3*h
x <- x[-2]
y <- y[-2]
xueya_nianling.lm <- lm(y~x)
summary(xueya_nianling.lm)
which(abs(rstudent(xueya_nianling.lm))>=3)
qt(0.975,27)
qf(0.95,1,27)
predict(xueya_nianling.lm ,data.frame(x=50),se.fit=T,interval='confidence',level=0.95)
predict(xueya_nianling.lm ,data.frame(x=50),se.fit=T,interval='prediction',level=0.95)
predict(xueya_nianling.lm ,data.frame(x=60),se.fit=T,interval='confidence',level=0.95)
#2
y=c(144,215,138,145,162,142,170,124,158,154,162,150,140,110,128,130,135,114,116,
124,136,142,120,120,160,158,144,130,125,175)
x1=c(39,47,45,47,65,46,67,42,67,56,64,56,59,34,42,48,45,18,20,19,36,50,39,21,44,
53,63,29,25,69)
x2=c(24.2,31.1,22.6,24.0,25.9,25.1,29.5,19.7,27.2,19.3,28.0,25.8,27.3,20.1,21.7,
22.2,27.4,18.8,22.6,21.5,25.0,26.2,23.5,20.3,27.1,28.6,28.3,22.0,25.3,27.4)
x3=c(0,1,0,1,1,0,1,0,1,0,1,0,0,0,0,1,0,0,0,0,0,1,0,0,1,1,0,1,0,1)
Pressure.lm <- lm(y~x1+x2+x3)
summary(Pressure.lm)
plot(Pressure.lm)
#从图上可以看出2号跟10号有些异常，现在用删除学生残差杆看一下
which(abs(rstudent(Pressure.lm))>=3)
#剔除2，10
y <- y[-c(2,10)]
x1 <- x1[-c(2,10)]
x2 <- x2[-c(2,10)]
x3 <- x3[-c(2,10)]
#重新做回归
Pressure.lm <- lm(y~x1+x2+x3)
summary(Pressure.lm)
plot(Pressure.lm)
#重新检测回归
which(abs(rstudent(Pressure.lm))>=3)
which(abs(cooks.distance(Pressure.lm))>1)
predict(Pressure.lm,data.frame(x1=43,x2=23.2,x3=1),se.fit=T,
interval='confidence',level=0.95)
predict(Pressure.lm,data.frame(x1=43,x2=23.2,x3=1),se.fit=T,
interval='prediction',level=0.95)
#正态假设检验，单个正态总体，方差未知且相等t检验
#一、右边检验
X<-c(159,280,101,212,224,379,179,264,222,362,168,250,149,260,485,170)
t.test(X,alternative="greater",mu=225)
#看p值，>0.05，不能拒绝原假设，接受H0
#计算t统计量qt(0.95,15)=1.75305，与上述得到的t=0.66852比较，0.66852<1.75305，接受H0
#alternative hypothesis:显示的是备择假设，由此我们可以知道原假设
#mean of x 上述数据x的均值
#二、左边检验
t.test(X,alternative="less",mu=225)
#p=0.743>0.05,接受H0
#计算qt(0.05,15),df表示自由度是15，0.66852>-1.75305，接受H0
#由于左边右边的假设结果一个是mu<=225,一个是mu>=225,所以做双边假设
#三、双边假设
t.test(X,alternative="two.sided",mu=225)
#qt(0.975,15)=2.13145,0.66852绝对值<2.13145,接受原假设，所以最后答案 mu=225
#综上，这个题三种检验都要做，因为前两种得到结果恰好相反，所以猜测只能取=号
#多个正态总体
#方法一
X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
t.test(X,Y,var.equal=TRUE, alternative="less")
#由p<0.05，拒绝原假设
#qt(0.05,18)=-1.734064， -4.2957<-1.734064,拒绝原假设
#方法二，letZ=X-Y，mu=0，Z跟mu比较，该换成单个正态总体
t.test(X-Y,mu=0, alternative="less")
#qt(0.05,9)=-1.833113 -4.2018<-1.833113 拒绝原假设
#p=0.00115<0.05 拒绝原假设
#神经网络模型
#蠓虫分类问题：现在有两种虫子Af、Apf，根据数据判断虫子属于哪一种
#每只虫子有两个数据，现在样本有两种虫子分别是9只、6只，根据样本建立模型并测试
#读入数据
MengChong<-read.table('MengChong.txt',sep=",",header = T)
#归一化
normalize <- function(x) {
return ((x - min(x)) / (max(x) - min(x)))
}
MengChong_n <- as.data.frame(lapply(MengChong[,1:2], normalize))
#取训练样本跟标签、测试样本跟标签
train_data<-as.matrix(MengChong_n);
test_data<-as.matrix(MengChong_n);
train_labels<-as.numeric(MengChong[,3])
test_labels<-train_labels
#newff构建神经网络
install.packages("AMORE")
library(AMORE)
net <- newff(n.neurons=c(2,8,2,1), learning.rate.global=1e-2, momentum.global=0.5,
error.criterium="LMS", Stao=NA, hidden.layer="tansig",
output.layer="purelin", method="ADAPTgdwm")
#神经网络有4层，依次是：‘2’输入层 2个圈、‘8’隐藏层 8个圈、‘2’隐藏层 2个圈、‘1’输出层 1个圈
#输入层的圈个数是看输入的每个样本有几个数据，输出层的圈个数是看输出的样本属于几个类别
#输入训练样本用该网络构建一个训练，并把该训练赋值为result
result <- train(net, train_data, train_labels, error.criterium="LMS",
report=TRUE, show.step=100, n.shows=5 )
#用测试样本测试上述训练，看用该训练得到的测试样本的结果是什么
y <- sim(result$net, test_data)
#比较测试样样本真实的结果跟用训练result产生的结果比较
y[which(y<1.5)] <- 1
y[which(y>=1.5)] <- 2
table(test_labels,y)
#计算该训练的正确率
n=length(test_labels)
Sum = 0
for(i in 1:n){
if(y[i]==test_labels[i]){
Sum =Sum+1
}
}
cat("正确率", (Sum/n)*100, "%")
#现在新来三只虫子的数据，想用训练判断他们的类别分别是什么
x<-rbind(c(1.24,1.80),c(1.28,1.84),c(1.40,2.04))
#用x储存这三只虫子的数据，一行是一只，第一列是antenna，第二列是wing
#对新来这三只做归一化，要用到原来训练数据各个属性值的最大和最小值
a<-x[,1]
b<-x[,2]
a1<-(a-min(MengChong[,1]))/(max(MengChong[,1])-min(MengChong[,1]))
b1<-(b-min(MengChong[,2]))/(max(MengChong[,2])-min(MengChong[,2]))
x_n<-cbind(a1,b1)
#输入归一化好的数据到训练中，用训练判断类别（即测试x_n）
y1 <- sim(result$net, x_n)
y1
#y1跟1接近，class就是1，跟2接近class就是2
#发现y1中的数据在1.5左右，跟1和2都不太接近，需要修改
#方法一：修改网络层中的第二层跟第三层的圈数，重新训练
net <- newff(n.neurons=c(2,8,1), learning.rate.global=0.5, momentum.global=0.5,
error.criterium="LMS", Stao=NA, hidden.layer="tansig",
output.layer="purelin", method="ADAPTgdwm")
#方法二 扩大样本数据
#方法三 给样本的每类数据加上不同的权重，关系大的一类权重大
#由于上述原数据相差不大，可以不用归一化，如下：
train_data<-as.matrix(MengChong[,1:2]);
test_data<-train_data;
train_labels<-c(1,1,1,1,1,1,1,1,1,2,2,2,2,2,2)
test_labels<-train_labels
library(AMORE)
net <- newff(n.neurons=c(2,8,2,1), learning.rate.global=1e-2, momentum.global=0.5,
error.criterium="LMS", Stao=NA, hidden.layer="tansig",
output.layer="purelin", method="ADAPTgdwm")
result <- train(net, train_data, train_labels, error.criterium="LMS",
report=TRUE, show.step=100, n.shows=5 )
y <- sim(result$net, test_data)
y[which(y<1.5)] <- 1
y[which(y>=1.5)] <- 2
table(test_labels,y)
n=length(test_labels)
Sum = 0
for(i in 1:n){
if(y[i]==test_labels[i]){
Sum =Sum+1
}
}
cat("正确率", (Sum/n)*100, "%")
x<-rbind(c(1.24,1.80),c(1.28,1.84),c(1.40,2.04))
y1 <- sim(result$net, x)
y1
#神经网络之拟合
#样本为每月的销售量，用来预测下一个月的销售量
sales<-read.table('sales.txt',header = F)
sales_n<-(sales-min(sales))/(max(sales)-min(sales))
sales_T<-sales_n[4:9,]
Train_data<-rbind(sales_n[1:3,],sales_n[2:4,],sales_n[3:5,],
sales_n[4:6,],sales_n[5:7,],sales_n[6:8,])#取训练样本的方法见图片
Train_labels<-sales_T
Test_data<-Train_data
Test_labels<-sales_T
library(AMORE)
set.seed(12345)
net <- newff(n.neurons=c(3,5,1), learning.rate.global=1e-2, momentum.global=0.5,
error.criterium="LMS", Stao=NA, hidden.layer="tansig",
output.layer="purelin", method="ADAPTgdwm")
result <- train(net, Train_data, Train_labels, error.criterium="LMS", report=TRUE,
show.step=100, n.shows=5 )
y <- sim(result$net, Test_data)
Y<-y*(max(sales)-min(sales))+min(sales)#返归一化，把压缩后的数据还原回去
real_Y<-sales[4:9,]#取出真实的值
plot(1:6,Y,xlim = c(1,6),ylim = c(1300,3000))
points(1:6,real_Y,col='blue',pch=15)#通过图形来比较真实值跟拟合值
?newff #可以看到激活函数的种类
#用另一个包来建立神经网络
Train_data1<-cbind(Train_data,Train_labels)
#install.packages("neuralnet")
library(neuralnet)
set.seed(12345) # to guarantee repeatable results
as.data.frame(Train_data1)->Train_data1
sales_model <- neuralnet(formula = Train_labels~ V1 + V2 + V3,data =Train_data1)
#关于neuralnet包，要求数据用dataframe的形式，
#要先as.data.frame(Train_data1)->Train_data1把Train_data1转换成数据框的形式
plot(sales_model)
model_results <- compute(sales_model, Test_data)
predicted_sales <- model_results$net.result
real_Y<-sales[4:9,]
y=predicted_sales
Y<-y*(max(sales)-min(sales))+min(sales)
plot(1:6,Y,xlim = c(1,6),ylim = c(1300,3000))
points(1:6,real_Y,col='blue',pch=15)
#sign()a符号函数，输入负数输出-1，输入正数输出1，输入0输出0
#计算例题output
#method 1
w1 <- 0.3
w2 <- 0.4
w3 <- 0.5
x1 <- 1.24
x2 <- 1.37
x3 <- 1.56
T <- 0.3
o <- sign(w1*x1+w2*x2+w3*x3-T)
#method 2
w0 <- T
x0 <- -1
W <- c(w0,w1,w2,w3)
X <- c(x0,x1,x2,x3)
o <- sign(W%*%X)
#iris_newff
iris<-iris
iris$Species<-factor(iris$Species,levels = c("setosa","versicolor","virginica"),
labels = c("1","2","3"))
iris$Species<-as.numeric(iris$Species)
#for (i in 1:5) {
# iris[,i] <- as.numeric(as.vector(iris)[,i])
#}
#normalize <- function(x) {
# return ((x - min(x)) / (max(x) - min(x)))
#}
#iris_n <- as.data.frame(lapply(iris[,1:4], normalize))
#train_data<-iris[c(1:40,51:90,101:140),1:4]
#test_data<-iris[c(41:50,91:100,141:150),1:4]
#train_labels<-iris[c(1:40,51:90,101:140),5]
#test_labels<-iris[c(41:50,91:100,141:150),5]
#另一种选训练测试样本的方法
train_index <- c(sample(1:50,40),sample(51:100,40),sample(101:150,40))#产生训练样本的行标
train_data <- iris[train_index,1:4]
test_data <- iris[-train_index,1:4]
train_labels <- iris[train_index,5]
test_labels <- iris[-train_index,5]
library(AMORE)
net <- newff(n.neurons=c(4,8,2,1), learning.rate.global=1e-2, momentum.global=0.5,
error.criterium="LMS", Stao=NA, hidden.layer="tansig",
output.layer="purelin", method="ADAPTgdwm")
#隐藏层8，2都是自己决定，主要注意4:iris一个数据有4个分量，1:最后判断种类，一个数据只对应一种
result <- train(net, train_data, train_labels, error.criterium="LMS", report=TRUE,
show.step=100, n.shows=5 )
y <- sim(result$net, test_data)
y[which(y<1.5)] <- 1
y[which(2<=y&y<3)] <- 2
y[which(y>=3)] <- 3
y<-as.numeric(y)
n=length(test_labels)
s<-0
for(i in 1:n){
if(y[i]==test_labels[i]){
s<-c(s+1)
}
}
cat("正确率", (s/30)*100,"%")
setwd('//Users//wangxinran//Desktop//R语言')
wine178_n <- read.csv('wine178.csv')
wine178_n$class<-as.numeric(wine178_n$class)
table(wine178_n[,14])
normalize <- function(x){
return((x-min(x))/(max(x)-min(x)))
}
wine178 <- as.data.frame(lapply(wine178_n[,-14],normalize))
train_index <- c(sample(1:59,50),sample(60:130,60),sample(130:178,40))#产生训练样本的行标
train_data <- wine178[train_index,1:13]
test_data <- wine178[-train_index,1:13]
train_labels <- wine178_n[train_index,14]
test_labels <- wine178_n[-train_index,14]
library(AMORE)
net <- newff(n.neurons=c(13,8,2,1), learning.rate.global=1e-2, momentum.global=0.5,
error.criterium="LMS", Stao=NA, hidden.layer="tansig",
output.layer="purelin", method="ADAPTgdwm")
result <- train(net, train_data, train_labels, error.criterium="LMS", report=TRUE,
show.step=100, n.shows=5 )
y <- sim(result$net, test_data)
y[which(y<1.5)] <- 1
y[which(2<=y&y<3)] <- 2
y[which(y>=3)] <- 3
y<-as.numeric(y)
n=length(test_labels)
s<-0
for(i in 1:n){
if(y[i]==test_labels[i]){
s<-c(s+1)
}
}
cat("正确率", (s/30)*100,"%")