首先安装几个R数据分析的包

funModeling:探索性数据分析(EDA)、数据准备和模型性能评估

tidyverse:数据科学集成包

Hmisc:高级计算函数、绘图功能等

matrixStats:数据分析

ggplot2:画图

载入需要的软件包

# library(funModeling) # library(tidyverse) # library(Hmisc) # library(matrixStats) library(fBasics) 计算偏度峰度 library(ggplot2) 

讯享网Loading required package: timeDate Loading required package: timeSeries 

读取文件

source<-read.table("exersice1_3.txt", header = TRUE, fileEncoding = "UTF8"); 

讯享网#数据概览 head(source) 

(1)求均值、方差、标准差、变异系数、偏度、峰度

 均值 n = dim(source)[1] # 行数 AA = source[, 2:4] ColMean <- apply(AA,2,mean) ColMean 

全国居民

1117

农村居民

747.6

城镇居民

2336.

讯享网 方差 ColVar <- apply(AA,2,var) ColVar 

全国居民

.

农村居民

.

城镇居民

.

 标准差 ColStd <- apply(AA,2,sd) ColStd 

全国居民

1015.

农村居民

632.2

城镇居民

2129.

讯享网 变异系数 ColCV <- ColStd/ColMean * 100 ColCV 

全国居民

90.8

农村居民

84.57

城镇居民

91.2

偏度公式：

$$g1=\frac{n}{(n-1)(n-2)}\frac{1}{s^3}\sum _{i=1}^n(x-\bar{x}{)}^3$$

峰度公式：
$$g2=\frac{n(n+1)}{(n-1)(n-2)(n-3)}\frac{1}{s^4}\sum _{i=1}^n(x_i-x{)}^4-\frac{3(n-1{)}^2}{(n-2)(n-3)}$$

 利用宏包fBasics计算偏度 skewness(AA) 

全国居民

0.8174

农村居民

0.7427

城镇居民

0.8679

讯享网 自定义函数计算公式 pianD<-function(xx=AA){ yy=sum((xx-mean(xx))^3)/sd(xx)^3*n/(n-1)/(n-2); return(yy) } apply(AA, 2, pianD) 

全国居民

1.018

农村居民

1.018

城镇居民

0.8573

 利用宏包fBasics计算峰度 kurtosis(AA) 

全国居民

-0.8152

农村居民

-0.53

城镇居民

-0.4999

讯享网 自定义函数计算峰度 fengD<-function(xx=AA){ aFeng=n*(n+1)/(n-1)/(n-2)/(n-3); bFeng=3*(n-1)^2/(n-2)/(n-3); yy=sum((xx-mean(xx))^4)/sd(xx)^4*aFeng-bFeng; return(yy) } apply(AA, 2, fengD) 

全国居民

-0.7379

农村居民

-0.318

城镇居民

-0.541

可以看出调用宏包计算结果与自定义函数计算结果不同，可能宏包不是无偏估计？

 画出核密度图 ggplot(source, aes(x=全国居民)) + geom_density(fill="#f4876a",colour="black",alpha=0.6,size=1) + geom_vline(data=source, aes(xintercept=mean(source[,2])), linetype="dashed", size=1) ggplot(source, aes(x=农村居民)) + geom_density(fill=9,colour="black",alpha=0.6,size=1) + geom_vline(data=source, aes(xintercept=mean(source[,3])), linetype="dashed", size=1) ggplot(source, aes(x=城镇居民)) + geom_density(fill=8,colour="black",alpha=0.6,size=1) + geom_vline(data=source, aes(xintercept=mean(source[,4])), linetype="dashed", size=1) 

讯享网

由图可以看出均属于左偏型的，即偏度大于零

讯享网 画出各列同方差正太分布于0-1正太分布的比较图 Norm = data.frame('全国居民'= rnorm(50000, mean = 0, sd = sd(source[,1])), norm = rnorm(50000, mean=0, sd=1)) ggplot(Norm) + geom_density(aes(x=全国居民), fill="#f4876a",colour="black",alpha=0.6,size=1) + geom_density(aes(x=norm), fill=9,colour="black",alpha=0.6,size=1) Norm = data.frame('农村居民'= rnorm(50000, mean = 0, sd = sd(source[,2])), norm = rnorm(50000, mean=0, sd=1)) ggplot(Norm, aes(x=农村居民)) + geom_density(fill="#f4876a",colour="black",alpha=0.6,size=1) + geom_density(aes(x=norm), fill=9,colour="black",alpha=0.6,size=1) Norm = data.frame('城镇居民'= rnorm(50000, mean = 0, sd = sd(source[,3])), norm = rnorm(50000, mean=0, sd=1)) ggplot(Norm, aes(x=城镇居民)) + geom_density(fill="#f4876a",colour="black",alpha=0.6,size=1) + geom_density(aes(x=norm), fill=9,colour="black",alpha=0.6,size=1) 

由图可以看出均属于细尾型的，即峰度小于零

(2)求中位数、上下四分位数、四分位极差、三均值

 分位数 summary(AA) 

讯享网 全国居民 农村居民 城镇居民 Min. : 184.0 Min. : 138.0 Min. : 405.0 1st Qu.: 321.8 1st Qu.: 255.2 1st Qu.: 617.8 Median : 727.5 Median : 530.5 Median :1499.5 Mean :1117.0 Mean : 747.9 Mean :2336.4 3rd Qu.:1642.2 3rd Qu.:1052.2 3rd Qu.:3675.0 Max. :3180.0 Max. :1973.0 Max. :6651.0 

 四分位极差 Q1 <- unname(quantile(AA[,1], probs=c(0.25, 0.5, 0.75))) Q2 <- unname(quantile(AA[,2], probs=c(0.25, 0.5, 0.75))) Q3 <- unname(quantile(AA[,3], probs=c(0.25, 0.5, 0.75))) R1 <- Q1[3] - Q1[1] R2 <- Q2[3] - Q2[1] R3 <- Q3[3] - Q3[1] R1; R2; R3; 

1320.5

797

3057.25

讯享网 三均值 M1 = 1/4*Q1[1] + 1/2*Q1[2] + 1/4*Q1[3] M2 = 1/4*Q2[1] + 1/2*Q2[2] + 1/4*Q2[3] M3 = 1/4*Q3[1] + 1/2*Q3[2] + 1/4*Q3[3] M1; M2; M3 

854.75

592.125

1822.9375

(3)作出直方图

 直方图 ggplot(source, aes(x=全国居民)) + geom_histogram(position= 'identity', stat="bin", binwidth=600) ggplot(source, aes(x=农村居民)) + geom_histogram(position= 'identity', stat="bin", binwidth=600) ggplot(source, aes(x=城镇居民)) + geom_histogram(position= 'identity', stat="bin", binwidth=800) 

(4)作出茎叶图

讯享网 茎叶图 stem(AA[,1],scale = 2,width=80,atom = 1e-08) stem(AA[,2],scale = 2,width=80,atom = 1e-08) stem(AA[,3],scale = 2,width=80,atom = 1e-08) sink() 

 The decimal point is 3 digit(s) to the right of the | 0 |  0 |  1 | 13 1 | 7 2 | 3 2 | 68 3 | 02 The decimal point is 2 digit(s) to the right of the | 0 | 468 2 | 025858 4 | 2157 6 | 22 8 | 6 10 | 2 12 | 14 | 3 16 | 7 18 | 807 The decimal point is 3 digit(s) to the right of the | 0 |  1 | 14679 2 | 4 3 | 09 4 | 9 5 | 48 6 | 27 

(5)找出异常值

讯享网计算上下截断点 except1 = list(Q1[1] - 1.5*R1, Q1[3] + 1.5*R1) except2 = list(Q2[1] - 1.5*R2, Q2[3] + 1.5*R2) except3 = list(Q3[1] - 1.5*R3, Q3[3] + 1.5*R3) except1; except2; except3 

1. -1659
2. 3623

1. -940.25
2. 2247.75

1. -3968.125
2. 8260.875

AA[, 1][AA[, 1] > except1[2] || AA[, 1] < except1[1]] 

讯享网AA[, 2][AA[, 2] > except2[2] || AA[, 2] < except2[1]] 

AA[, 3][AA[, 3] > except3[2] || AA[, 3] < except3[1]] 

可以看到都为空值，所以没有异常值

讯享网boxplot(AA[, 1])$out 

boxplot(AA[, 2])$out 

讯享网boxplot(AA[, 3])$out 

由箱图也可以看出没有异常点