模拟产生统计专业同学的名单(学号区分),记录数学分析,线性代数,概率统计三科成绩,然后进行一些统计分析
先简单介绍几个函数:
正态分布函数rnorm()
泊送分布函数rpois()
指数分布函数rexp()
Gamma分数函数rgamma()
均匀分布函数runif()
二项分布函数rbinom()
几何分布函数rgeom()
> num<-seq(10378001,10378100)
> num
[1] 10378001 10378002 10378003 10378004 10378005 10378006 10378007 10378008
[9] 10378009 10378010 10378011 10378012 10378013 10378014 10378015 10378016
[17] 10378017 10378018 10378019 10378020 10378021 10378022 10378023 10378024
[25] 10378025 10378026 10378027 10378028 10378029 10378030 10378031 10378032
[33] 10378033 10378034 10378035 10378036 10378037 10378038 10378039 10378040
[41] 10378041 10378042 10378043 10378044 10378045 10378046 10378047 10378048
[49] 10378049 10378050 10378051 10378052 10378053 10378054 10378055 10378056
[57] 10378057 10378058 10378059 10378060 10378061 10378062 10378063 10378064
[65] 10378065 10378066 10378067 10378068 10378069 10378070 10378071 10378072
[73] 10378073 10378074 10378075 10378076 10378077 10378078 10378079 10378080
[81] 10378081 10378082 10378083 10378084 10378085 10378086 10378087 10378088
[89] 10378089 10378090 10378091 10378092 10378093 10378094 10378095 10378096
[97] 10378097 10378098 10378099 10378100
模拟数学分析成绩
> x1<-round(runif(100,min=80,max=100))#80到100分之间均匀分布,然后四舍五入
> x1
[1] 83 95 81 100 87 85 94 89 91 91 83 85 99 98 82 81 80 96
[19] 82 95 84 85 85 84 98 87 91 99 87 84 91 84 84 81 96 95
[37] 97 94 88 89 82 99 86 93 100 93 83 90 99 92 94 88 94 81
[55] 90 81 86 83 93 93 100 88 87 94 85 85 94 97 92 87 86 89
[73] 90 87 80 97 83 94 87 81 83 84 97 91 92 98 93 83 80 97
[91] 91 82 91 89 86 85 90 90 94 87
模拟线性代数成绩
> x2<-round(rnorm(100,mean=80,sd=7))#产生100个正态分布,平均值是80,标准差是7,然后四舍五入
> x2
[1] 74 89 87 79 83 75 77 79 72 96 86 86 79 79 75 72 73 80
[19] 85 87 73 77 75 82 79 81 84 78 68 73 90 79 79 87 78 71
[37] 70 84 75 97 84 74 85 78 78 83 70 75 78 75 81 82 83 75
[55] 73 81 82 75 73 81 78 100 77 77 83 79 80 82 80 76 79 78
[73] 79 73 78 74 77 69 78 85 77 74 80 77 84 81 79 86 66 80
[91] 93 75 87 79 73 78 77 86 78 86
模拟概率统计成绩
> x3<-round(rnorm(100,mean=83,sd=18))#产生100个正态分布,平均值是83,标准差是18,然后四舍五入
> x3
[1] 68 88 113 89 98 82 101 74 76 115 82 63 114 73 88 92 103 65
[19] 77 89 79 73 101 85 97 100 108 52 50 78 90 83 80 120 107 92
[37] 75 73 119 47 85 94 49 74 126 79 95 90 76 120 74 71 119 57
[55] 89 83 78 113 73 53 101 83 78 91 88 77 87 78 92 80 106 106
[73] 67 65 64 98 59 98 56 69 51 100 69 67 72 79 83 89 93 67
[91] 116 100 71 86 65 59 84 79 53 71
> x3[which(x3>100)]=100#对产生的值中超过100的改写成100
> x3
[1] 68 88 100 89 98 82 100 74 76 100 82 63 100 73 88 92 100 65
[19] 77 89 79 73 100 85 97 100 100 52 50 78 90 83 80 100 100 92
[37] 75 73 100 47 85 94 49 74 100 79 95 90 76 100 74 71 100 57
[55] 89 83 78 100 73 53 100 83 78 91 88 77 87 78 92 80 100 100
[73] 67 65 64 98 59 98 56 69 51 100 69 67 72 79 83 89 93 67
[91] 100 100 71 86 65 59 84 79 53 71
合成数据框
> x<-data.frame(num,x1,x2,x3)
> x
num x1 x2 x3
1 10378001 83 74 68
2 10378002 95 89 88
3 10378003 81 87 100
4 10378004 100 79 89
5 10378005 87 83 98
6 10378006 85 75 82
7 10378007 94 77 100
8 10378008 89 79 74
9 10378009 91 72 76
10 10378010 91 96 100
11 10378011 83 86 82
12 10378012 85 86 63
13 10378013 99 79 100
14 10378014 98 79 73
15 10378015 82 75 88
16 10378016 81 72 92
17 10378017 80 73 100
18 10378018 96 80 65
19 10378019 82 85 77
20 10378020 95 87 89
21 10378021 84 73 79
22 10378022 85 77 73
23 10378023 85 75 100
24 10378024 84 82 85
25 10378025 98 79 97
...
保存到硬盘
> write.table(x,file="d:\\mark.txt",col.names=F,row.names=F,sep=" ")#将数据保存到d:\mark.txt中,不保存行名和列名,数据以空格分割开来
对列求平均值
> colMeans(x)[c("x1","x2","x3")]
x1 x2 x3
89.26 79.37 81.46
apply()函数说明,举例apply(x,2,mean)中的第一个参数表示传入的数据,第二个参数表示是对行还是列进行运算,1表示行,2表示列,第三个参数表示所求的形式,mean求均值,max求最大值
所以求成绩的平均值,最大值,最小值可以这样(这里的缺点是把学号也求进去了):
> apply(x,2,mean)
num x1 x2 x3
10378050.50 89.26 79.37 81.46
> apply(x,2,max)
num x1 x2 x3
10378100 100 100 100
> apply(x,2,min)
num x1 x2 x3
10378001 80 66 47
>
当然去掉学号我们可以这样做:
> apply(x[c("x1","x2","x3")],2,mean)
x1 x2 x3
89.26 79.37 81.46
> apply(x[c("x1","x2","x3")],1,sum)
[1] 225 272 268 268 268 242 271 242 239 287 251 234 278 250 245 245 253 241
[19] 244 271 236 235 260 251 274 268 275 229 205 235 271 246 243 268 274 258
[37] 242 251 263 233 251 267 220 245 278 255 248 255 253 267 249 241 277 213
[55] 252 245 246 258 239 227 278 271 242 262 256 241 261 257 264 243 265 267
[73] 236 225 222 269 219 261 221 235 211 258 246 235 248 258 255 258 239 244
[91] 284 257 249 254 224 222 251 255 225 244
> which.max(apply(x[c("x1","x2","x3")],1,sum))
[1] 10
> x$num[which.max(apply(x[c("x1","x2","x3")],1,sum))]
[1] 10378010
绘制直方图函数hist()
> hist(x$x1)
图形如下:
> plot(x$x1,x$x2)
图形如下:
列联函数table(),柱状图绘制函数barplot()
table()计算各分数有多少人
> table(x$x1)
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
3 6 4 7 6 7 4 8 3 4 5 7 3 5 8 3 2 5 3
99 100
4 3
barplot()绘制柱状图
> barplot(table(x$x1))
图形如下:
饼图绘制函数pie()
> pie(table(x$x1))
图形如下:
箱子的上下横线为样本的25%和75%分位数
箱子中间的横线为样本的中位数
上下延伸的直线称为尾线,尾线的尽头为最高值和最低值
异常值
> boxplot(x$x1,x$x2,x$x3)
图形如下:
水平放置的箱尾图:
> boxplot(x$x1,x$x2,x$x3,horizOntal=T)
图形如下:
> boxplot(x[2:4],col=c("red","green","blue"),notch=T)
图形如下:
每个观测单位的数值表示为一个图形每个图的每个角表示一个标量,字符串类型会标注在图的下方
角线的长度表达值的大小
> stars(x[c("x1","x2","x3")])
图形如下:
> stars(x[c("x1","x2","x3")],full=T,draw.segment=T)
> stars(x[c("x1","x2","x3")],full=F,draw.segment=T)
用五官的宽度和高度来描绘数值
人对脸谱高度敏感和强记忆
适合较少样本的情况
脸谱图不在R的基本包中,需要额外的安装aplpack包
> library(aplpack)
载入需要的程辑包:tcltk
> faces(x[c("x1","x2","x3")])
effect of variables:
modified item Var
"height of face " "x1"
"width of face " "x2"
"structure of face" "x3"
"height of mouth " "x1"
"width of mouth " "x2"
"smiling " "x3"
"height of eyes " "x1"
"width of eyes " "x2"
"height of hair " "x3"
"width of hair " "x1"
"style of hair " "x2"
"height of nose " "x3"
"width of nose " "x1"
"width of ear " "x2"
"height of ear " "x3"
其他的脸谱图
安装TeachingDemos包
> library(TeachingDemos)
载入程辑包:‘TeachingDemos’
The following objects are masked from ‘package:aplpack’:
faces, slider
> faces2(x)
用数字和字符来表示
> stem(x$x1)
The decimal point is at the |
80 | 000000000
82 | 00000000000
84 | 0000000000000
86 | 000000000000
88 | 0000000
90 | 000000000000
92 | 00000000
94 | 00000000000
96 | 0000000
98 | 0000000
100 | 000
可用于判断是否正态分布
直线的斜率是标准差,截距是均值
点的散布越接近直线,则越接近正态分布
> qqnorm(x1)
> qqline(x1)
散点图的进一步设置
> plot(x$x1,x$x2,main="数学分析与线性代数成绩的关系",xlab="数学分析",ylab="线性代数",xlim=c(0,100),ylim=c(0,100),xaxs="i",yaxs="i",col="red",pch=19)
图形如下:
连线图
> a<-c(2,3,4,5,6)
> b<-c(4,7,8,9,12)
> plot(a,b,type="l")
图形如下:
密度图
> plot(density(rnorm(1000)))
R语言有部分内置的数据集,使用函数data()可以查看
> data()
如下:
利用内置的mtcars数据集绘制
> heatmap(as.matrix(mtcars),Rowv=NA,Colv=NA,col=heat.colors(256),scale="column",margins=c(2,8),main="Car characteristics by Model")
图形如下:
Sepal花萼
Petal花瓣
Species种属
> iris
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
7 4.6 3.4 1.4 0.3 setosa
8 5.0 3.4 1.5 0.2 setosa
9 4.4 2.9 1.4 0.2 setosa
10 4.9 3.1 1.5 0.1 setosa
11 5.4 3.7 1.5 0.2 setosa
12 4.8 3.4 1.6 0.2 setosa
13 4.8 3.0 1.4 0.1 setosa
14 4.3 3.0 1.1 0.1 setosa
15 5.8 4.0 1.2 0.2 setosa
....
用来克服散点图中数据点重叠问题
在有重叠的地方用一朵“向日葵花”的花瓣数目来表示重叠数据的个数
> sunflowerplot(iris[,3:4],col="gold",seg.col="gold")
遍历样本中全部的变量配对画出二元图
直观地了解所有变量之间的关系
> pairs(iris[,1:4])
用plot也可以实现相同的效果
> plot(iris[,1:4],main="Relationships between characteristics of iris flowers",pch=19,col="blue",cex=0.9)
图形如下:
利用par()在同一个device输出多个散点图
par命令博大精深,用于设置绘图参数,help(par)
> par(mfrow=c(3,1))
> plot(x1,x2)
> plot(x2,x3)
> plot(x3,x1)
图形如下:
有哪些颜色?colors()
> colors()
[1] "white" "aliceblue" "antiquewhite"
[4] "antiquewhite1" "antiquewhite2" "antiquewhite3"
[7] "antiquewhite4" "aquamarine" "aquamarine1"
[10] "aquamarine2" "aquamarine3" "aquamarine4"
[13] "azure" "azure1" "azure2"
[16] "azure3" "azure4" "beige"
[19] "bisque" "bisque1" "bisque2"
[22] "bisque3" "bisque4" "black"
[25] "blanchedalmond" "blue" "blue1"
[28] "blue2" "blue3" "blue4"
[31] "blueviolet" "brown" "brown1"
[34] "brown2" "brown3" "brown4"
[37] "burlywood" "burlywood1" "burlywood2"
[40] "burlywood3" "burlywood4" "cadetblue"
[43] "cadetblue1" "cadetblue2" "cadetblue3"
....
dev.new()生成一个新的图形框
dev.list()查看有多少个图形框
dev.cur()查看当前的图形框是哪个
dev.next()查看下一个图形框
安装scatterplot3d包
> library(scatterplot3d)
> scatterplot3d(x[2:4])
三维作图
> x<-y<-seq(-2*pi,2*pi,pi/15)
> f<-function(s,y)sin(x)*sin(y)
> z<-outer(x,y,f)
> contour(x,y,z,col="blue")
> persp(x,y,z,theta=30,phi=30,expand=0.7,col="lightblue")
安装maps包
> library(maps)
> map("state",interior=FALSE)
图形如下:
> map("state",boundary=FALSE,col="red",add=TRUE)
> map("world",fill=TRUE,col=heat.colors(10))
先下载安装maps包和geosphere包并加载
加载包:
> library(maps)
> library(geosphere)
画出美国地图:
> map("state")
图形如下:
通过设置坐标范围使焦点集中在美国周边,并且设置一些有关颜色
> xlim<-c(-171.738281,-56.601563)
> ylim<-c(12.039321,71.856229)
> map("world",col="#f2f2f2",fill=TRUE,bg="white",lwd=0.05,xlim=xlim,ylim=ylim)
图形如下:
画一条弧线连线,表示社交关系
> lat_ca<-39.164141
> lon_ca<--121.64062
> lat_me<-45.21300
> lon_me<--68.906250
> inter<-gcIntermediate(c(lon_ca,lat_ca),c(lon_me,lat_me),n=50,addStartEnd=TRUE)
> lines(inter)
图形如下:
装载数据
> airports<-read.csv("http://datasets.flowingdata.com/tuts/maparcs/airports.csv",header=TRUE)
> flights<-read.csv("http://datasets.flowingdata.com/tuts/maparcs/flights.csv",header=TRUE,as.is=TRUE)
画出多重联系
> map("world",col="#f2f2f2",fill=TRUE,bg="white",lwd=0.05,xlim=xlim,ylim=ylim)
> fsub<-flights[flights$airline=="AA",]
> for(j in 1:length(fsub$airline)){
+ air1<-airports[airports$iata==fsub[j,]$airport1,]
+ air2<-airports[airports$iata==fsub[j,]$airport2,]
+ inter<-gcIntermediate(c(air1[1,]$long,air1[1,]$lat),c(air2[1,]$long,air2[1,]$lat),n=100,addStartEnd=TRUE)
+ lines(inter,col="black",lwd=0.8)
+ }
图形如下: