r语言关于偏度的计算代码 r语言计算偏度峰度

首先安装几个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)

A data.frame: 6 × 4

	年份	全国居民	农村居民	城镇居民
	<int>	<int>	<int>	<int>
1	1978	184	138	405
2	1979	207	158	434
3	1980	236	178	496
4	1981	262	199	562
5	1982	284	221	576
6	1983	311	246	603

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

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

1117 农村居民 747.863636363636 城镇居民 2336.40909090909

## 方差
ColVar <- apply(AA,2,var)
ColVar

1031680.28571429 农村居民 399673.837662338 城镇居民 4536136.44372294

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

1015.71663652531 农村居民 632.197625479832 城镇居民 2129.82075389525

## 变异系数
ColCV <- ColStd/ColMean * 100
ColCV

90.932554747118 农村居民 84.5338100076357 城镇居民 91.157869663422

偏度公式：

$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^n _{i=1} (x_i -x )^4- \frac{3(n-1)^2}{(n-2)(n-3)}$

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

0.889330205728174 农村居民 0.878668808067427 城镇居民 0.842138275058679

## 自定义函数计算公式
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.02484718945818 农村居民 1.01256119786818 城镇居民 0.970464107448573

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

-0.831996482068152 农村居民 -0.8278408069353 城镇居民 -0.915094641514999

## 自定义函数计算峰度
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.457241207817379 农村居民 -0.45144409308318 城镇居民 -0.57316209891541

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

## 画出核密度图
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 | 22233344
  0 | 567889
  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 | 445666789
  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
3. -940.25
4. 2247.75
5. -3968.125
6. 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

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