• Posts tagged "抽样"

Blog Archives

抽样三步曲之二:R语言进行抽样估计

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-sampling-estimation

前言

在国家大选前夕,媒体进行民意调查,通过对不同地区不同职业的400人进行调查,其中的160人支持A,根据调查结果,媒体给出结论,A的民众支持度为40%。这个结论是否合理呢?

本文将继续上一篇文章抽样三步曲之一:R语言进行随机抽样sampling,将介绍如何使用抽样结果,进行抽样估计,推断出科学的结论。抽样三步曲系列文章:用R语言进行随机抽样用R语言进行抽样估计和用R语言进行抽样校准。

目录

  1. 抽样估计介绍
  2. 点估计
  3. 抽样误差
  4. 区间估计

1. 抽样估计介绍

抽样估计(Sampling estimation)又称为抽样推断,是在抽样调查的基础上所进行的数据推测,即用样本数量特征来估计和推算总体数量特征。抽样估计,是对总体进行描述的另一种重要方法,具有花费小、适用性强、科学性高等特点。抽样估计,只能得到对总体特征近似的测试,因此抽样估计还必须同时考察所得结果的可能范围和可靠程度。

在现实场景中,我们对总体数据进行抽样得到的样本,分析样本的统计特征形成样本特征,然后利用样本特征对总体特征做一个估计或者推理。比如,采用样本均值对总体均值进行估计,采用样本方差对总体的方差进行估计,采用样本中某个特征的元素对总体中具有这个特征的元素进行估计等。

在统计学中,通常把总体特征也称为总体参数,总体参数是总体的某种特征值,是所研究的总体某些方面的数量表现。总体特征比如有总体均值,总体总值,总体比例,总体比率,总体方差等。统计估计,就是用对应的样本特征来估计出总体的特征。

抽样估计,使用大数定律做为用样本估计总体的基础理论,当样本量增大时,随机影响会相互抵消,样本均值与总体均值接近。使用中心极限定理进行区间估计的基础理论,只要方差有限,在样本足够多时,估计量的分布将趋于正态分布。

抽样估计,主要有两种方法点估计和区间估计。点估计,是通过样本估计量的某一次估计值来推断总体参数的可能取值。区间估计,是根据样本估计量以一定的可靠程序推断总体参数在所在区间的范围。

2. 点估计

点估计,又称定值估计,就是用实际样本指标数值作为总体参数的估计值。点估计的方法简单,一般不考虑抽样误差和可靠程度,它适用于对推断准确程度与可靠程度要求不高的情况。

比如:支持率调查,从总体选民4000人中抽样400个人调查,其中有160人投票支持A,以此推断A在总体选民的支持率为40%(160/400)。这就是一个典型的点估计,用某个实际的抽样结果来推断总体的结果,但准确度可能不高。

投票支持率就是一个点估计的统计量,一个好的点估计统计量,必须满足具有3个性质:无偏性,一致性,有效性。

  • 无偏性。即要求所有可能样本指标的平均数(样本指标的数学期望)与被估计的总体参数之间没有偏差。虽然每一次的样本指标值和总体指标值之间都可能有误差,但在多次反复的估计中,所有抽样指标值的平均数应该等于所估计的总体指标值本身,即用样本指标去估计总体参数,平均说来是没有偏误的。
  • 一致性。用统计量估计总体参数要求当样本的单位数充分大时,抽样指标也充分地靠近总体指标。就是说,随着样本单位数n的无限增加,统计量和未知的总体参数之差的绝对值小于任意小的数,它的概率也趋近于1,即实际上是几乎肯定的。
  • 有效性。以统计量估计总体参数时,优良估计量的方差应该比其他估计量的方差小。例如用样本平均数或总体某一变量值来估计总体平均数,虽然两者都是无偏的,而且在每一次估计中,两种估计量和总体平均数都可能有离差,但样本平均数更靠近于总体平均数的周围,平均说来其离差比较小。所以对比说来,抽样平均数是更为有效的估计量。

显然上例中的,基于抽样选民投票率推断出40%支持率的结论,是不够科学合理的。我们需要设计更多的统计量,来对总体特征进行估计。

接下来,我们设计一个场景,来完成抽样估计的实验过程。假设总体选民4000人,分布在全国的20个地区,随机抽样400人(无重复)进行调查。计算抽样特征:样本平均工资年薪,样本工资方差,样本投票支持率。1次抽样的结果,可能不能说明问题,我们进行500次抽样,把每次抽样结果汇总,可以形成抽样特征的抽样分布。以样本特征的平均值,来估计总体特征,得出结论。

操作步骤如下:

  1. 生成总体数据
  2. 抽样获得样本,并计算样本特征
  3. 多次抽样
  4. 抽样分布
  5. 获得结论

1. 生成总体数据
生成总体选民4000人数据。

# N总体4000人
> N<-4000

# 工资,使用均匀分布和正态分布进行叠加
> set.seed(0)
> salary1<-runif(2000,30000,80000)  # 均匀分布
> set.seed(0)
> salary2<-rnorm(2000,51814,3000)   # 正态分布
> salary<-c(salary1,salary2)

# 20个地区
> set.seed(1)                 #设置随机种子
> area<-sample(1:20,N,TRUE)

# 投票,1为投票给A,0为不投票给A
> set.seed(2)                 #设置随机种子
> vote<-sample(1:N,1500)

# 总体数据集
> df<-data.frame(salary,area)
> df$area<-as.factor(df$area)
> df$vote<-0
> df$vote[vote]<-1
> df$vote<-as.factor(df$vote)

# 查看数据前10条
> head(df,10)
     salary area vote
1  74834.86    4    0
2  43275.43    7    0
3  48606.19    1    1
4  58642.67    2    0
5  75410.39   11    0
6  40084.10   14    0
7  74919.48   18    0
8  77233.76   19    1
9  63039.89    1    0
10 61455.70   10    1

# 查看统计概率
> summary(df)
     salary           area      vote    
 Min.   :30030   13     : 226   0:2500  
 1st Qu.:48302   20     : 217   1:1500  
 Median :52015   9      : 215           
 Mean   :53254   14     : 215           
 3rd Qu.:56654   8      : 211           
 Max.   :79997   6      : 206           
                 (Other):2710   

总体数据集,包括3列:

  • salary,每个选民的工资
  • area,每个选民所在地区,1到20
  • vote,对A的投票

2. 抽样获得样本,并计算样本特征
进行随机抽样,选出400人。

# 加载工具包
> library(sampling)
> library(plyr)
> library(magrittr)

# 抽样400人
> s<-400
> set.seed(30)             # 随机种子
> s1<-srswor(s,N)      # 简单随机抽样

# 样本数据集 
> sample_df<-getdata(df,s1)

# 样本工资均值
> mean(sample_df$salary)
[1] 53936.74

# 样本工资标准差
> sd(sample_df$salary)
[1] 10930.51

# 样本工资方差
> var(sample_df$salary)
[1] 119476082

# 样本投票比例
> table(sample_df$vote)
  0   1 
242 158

# 地区分布
> table(sample_df$area)
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
12 21 24 17 22 17 20 25 19 17 20 15 24 24 18 22 20 23 26 14 

对选样本选民资进行正态分布检验。


> shapiro.test(sample_df$salary)
	Shapiro-Wilk normality test
data:  sample_df$salary
W = 0.94848, p-value = 1.389e-10

# 画出qq图
> qqnorm(sample_df$salary)
> qqline(sample_df$salary,col="red")

由于p-value=0.1389 > 0.05,不能拒绝正在分布的假设,所以,工资符合正态分布。但从图中看,已经偏离很大了。这样我们就完成一次抽样的过程了。
3. 多次抽样
1次抽样不能说明什么问题,我们进行多次抽样,进行500次随机抽样,可能得到不同的结果,形成不同的抽样分布。


# 定义抽样数据框
> sample500<-data.frame(time=1:500,mean=0,var=0,vote=0)

# 进行500次随机抽样
> for(i in 1:500){
+   s1<-srswor(s,N)
+   sample_df<-getdata(df,s1)
+   sample500$mean[i]<-mean(sample_df$salary)
+   sample500$var[i]<-var(sample_df$salary)
+   sample500$vote[i]<-length(which(sample_df$vote==1))/s
+ }

# 查看前6条数据
> head(sample500)
  time     mean       var   vote
1    1 53457.14 118965345 0.4000
2    2 53138.77  98917316 0.4150
3    3 52852.20 113358728 0.3925
4    4 53655.00 104086558 0.3525
5    5 52530.54 107175129 0.3625
6    6 53491.33 119491300 0.3700

这时,我们完成了500次的抽样,并把抽样的结果进行了整理。
4. 抽样分布
为了观察方例,以柱状图进行可视化,查看每个统计量的抽样分布。


# 图片输出布局,一行三列
> par(mfrow=c(1,3))  

# 输出样本均值,样本方差,样本投票率的柱状图
hist(sample500$mean,10,ylim=c(0,200),main='mean')
hist(sample500$var,10,ylim=c(0,200),main='var')
hist(sample500$vote,10,ylim=c(0,200),main='vote')

从图中可以观察到样本均值、样本方差和投票支持率,都看起来都呈现正态分布的形状。根据抽样分布性质,当总体分布已知且为正态分布或接近正态分布时,无论样本容量大小,样本均值都为正态分布。当总体分布未知,根据中心极限定理,样本容量足够大时趋近于正态分布。经验上验证,当样本容量大于30时,样本均值接近正态分布。

接下来,我们分别进行正态分布的检查。

# 均值的正态分布检查
> shapiro.test(sample500$mean)
	Shapiro-Wilk normality test
data:  sample500$mean
W = 0.99745, p-value = 0.6438

# 方差的正态分布检查
> shapiro.test(sample500$var)
	Shapiro-Wilk normality test
data:  sample500$var
W = 0.99503, p-value = 0.1084

# 投票率的正态分布检查
> shapiro.test(sample500$vote)
	Shapiro-Wilk normality test
data:  sample500$vote
W = 0.99465, p-value = 0.07877

左一图,500次抽样的工资的均值,为正态分布,而且p-value远大于0.05,比之前1次的抽样效果好很多。

 5. 获得结论

我们检验的结果与性质相符合,由于样本特征都满足正态分布,所以我们以这3个样本特征的均值,做为总体特征的点估计值,说服力就会更强一些。

> colMeans(sample500[,c("mean","var","vote")])
        mean          var         vote 
5.323558e+04 1.136433e+08 3.751550e-01 

计算500次抽样的样本特征的均值,民选平均年薪为53235.58元,对A的投票支持率为37.52%。

3. 抽样误差

在抽样调查中,通常调查结果和实际结果之间会有一定的偏差,根据形成原因不同,一般会分抽样误差和非抽样误差。

非抽样误差可以产生于抽样调查的各个阶段,包括抽样调查及抽样设计、数据采集及数据的处理与分析阶段。非抽样误差的产生不是由于抽样的随机性,而是由于其他多种原因引起的,主观因素偏多。

抽样误差,表示总体的点估计与总参数之间的差异,是由样本统计量推断总体参数时的误差,是不可避免的。抽样误差,会随着样本量增大而减小,受总体变异大小影响,抽样方法不同也会产生误差。

在上面的场景中,用总体数据的实际值与点估计值进行对比。

数据集数量salary均值vote投票支持率
总体400053253.5537.5%
1次抽样40054170.9839.5%
500次抽样400*50053235.5837.52%

我们观察到,1次抽样得到的点估计值与总体参数之间是存在较大的误差,这个误差被称为抽样误差。而进行500次的抽样后的估计值与总体参数之间,误差明显缩小了许多。

抽样分布的目的,就是为了确保我们在使用点估计值进行总体估计的时候能够清楚的知道误差的范围到底有多大,该如何去调整抽样的大小或采取相应的校正以使得其可以更加准确的近似总体的参数。

对抽样效果的评估时,当样本容量相同,抽样的样本误差越小,则效果越好。我们可以使用样本均值的标准差统计量,描述样本的误差,对抽样效果进行评价。样本均值的标准差可用来度量均值与总体均值的距离,即为可能的误差,又被称为均值标准误。


N为总体数量,s为样本数量

均值标准误,有限总体:
sd_sample = sqrt((N-s)/(N-1)) * (sd_total / sqrt(s))
其中,fpc = sqrt((N-s)/(N-1))为有限总体修正因子

均值标准误,无限总体:
sd_sample = sd_total / sqrt(s)

如果总体很大样本很小,当s/N<=0.05时,样本不足5%时,有限总体修正因子趋近于1.

计算场景中均值标准误。


> N<- 4000
> s<- 400

# 在有限个总体的情况下,修正因子 
> fpc<-sqrt((N-s)/(N-1)) ;fpc
[1] 0.9488019

# 均值标准误
> fpc * sd(df$salary) / sqrt(sd(sample500$mean))
[1] 447.9483

抽样误差,为样本统计量的估计值与要测度的总体参数值之间的绝对差距。抽样分布能够用来提供抽样误差大小的可能(概率)。


# 总体均值为红色线
> total_mean<-mean(df$salary)
> par(mfrow=c(1,4))
> hist(sample500$mean[1:30],breaks=20,xlim=c(51000,55000),ylim=c(0,80),main="30")
> abline(v = total_mean,col="red")
> hist(sample500$mean[1:100],breaks=20,xlim=c(51000,55000),ylim=c(0,80),main="100")
> abline(v =total_mean,col="red")
> hist(sample500$mean[1:200],breaks=20,xlim=c(51000,55000),ylim=c(0,80),main="200")
> abline(v =total_mean,col="red")
> hist(sample500$mean[1:500],breaks=20,xlim=c(51000,55000),ylim=c(0,80),main="500")
> abline(v =total_mean,col="red")

在上面场景中,抽样的次数越多,样本的均值越集中于实际的均值。如果媒体认为被抽样的选中的选民,他们的工资与总体的工资均值(53253.55)误差在200元内是合理的抽样,那么我们可以把问题转化为,计算标准正态分布工资在200元以内的可能性是多少,即抽样落在[53253.55-200, 53253.55+200]区间的概率。


# 标准正态分布概率
> ((53253.55+200)-(53253.55-200))/447.9483
[1] 0.8929602

# 计算面积
> pnorm(0.8929602, 0, 1)
[1] 0.8140608

在上述的简单随机抽样中,有81.41%的概率,使得样本均值与总体实际均值误差不超过200元。所以,上面的场景中,媒体有81.41%的概率会抽中他们目标的选民,1次抽样会造成偏差,当经过500次的抽样后,取投票率的平均值是科学合理的。

4. 区间估计

由于点估计始终会存在一定的估计误差,对于未知参数点估计值只是一个近似值。那么我们在进行估计时,不要给出一个具体值,而是给出一个区间值,可能是更合适,也是更可信的,就能解决估计不准的问题。

举个例子,去欧洲旅游10日自由行,你做预算时,可能给不出准确的花费18000元还20000元,你通常会说一个大概的范围,比如15000元到30000元之间,也可能是18000元到18500元之间。如果给出的范围太大,虽然可信度高,但没有什么参考意义;如果给出的范围大小,准确度高了,但是可信度又变低了。找到一个合适的估计范围,可以置信区间来解决。

区间估计,是根据样本指标、抽样误差和概率保证程度去推断总体参数的可能范围。在统计学中,通常用置信区间及样本出现的概率来估计总体参数,以一定的概率保证总体参数包含在估计区间内。区间估计要完成两个方面的估计:1是根据样本指标和抽样平均误差估计总体指标的可能范围,2是估计推断总体指标真实值在这个范围的可靠程度。

置信区间涉及两个问题,一个是置信水平,另一个是如何建立置信区间。所谓置信水平就是给出一个区间的信心,这个信心以概率来表示,绝大多数情况下取0.95,表示你对所估计的总体参数有95%的信心落在你所给的区间内。通常置信水平以1-alpha表示,alphe称为显著性水平。

继续上面欧洲旅游例子:我们从网上找到去过自游行去过欧洲旅游的人,找到了30个人为样本,样本的平均花费是为19377元,最低的5000元,最高的65000元,标准差为14750。我们以这个样本的95%的置信区间为标准计算,推算本次旅游的预算。

首先,假设总体的均值为20000元,方差为5000,生成样本数据集30个样本。


> total_m<-20000              # 总体均值
> total_sd<-5000              # 总体标准差
> s<-30                 # 样本数据
> min<-5000             # 样本最小值 
> max<-65000            # 样本最大值

# 生成样本数据
> set.seed(10)
> price<-rnorm(s-2,total_m,total_sd)          # 正态分布

# 查看样本数据
> price<-c(min,max,price);price               # 拼接最大值,最小值
 [1]  5000.000 65000.000 20093.731 19078.737 13143.347
 [6] 17004.161 21472.726 21948.972 13959.619 18181.620
[11] 11866.637 18717.608 25508.898 23778.908 18808.832
[16] 24937.224 23706.951 20446.736 15225.281 19024.248
[21] 24627.606 22414.893 17018.447  9073.566 16625.670
[26]  9404.694 13674.010 18131.692 16562.223 15639.206

# 样本均值
> sample_m<-mean(price);sample_m
[1] 19335.87

# 样本标准差
> sample_sd<-sd(price);sample_sd
[1] 9946.058

画出样本数据,x轴为样本点序号,y轴为旅游的花费

# 排序
> price<-price[order(price)]

# 画出散点图
> plot(price,type='both')

接下来,我们会分几种情况,进行区间估计。

  • 大样本(>=30),已知总体方差,对总体均值进行区间估计。
  • 大样本(>=30),未知总体方差,对总体均值进行区间估计。
  • 小样本(<30),总体符合正态分布时,已知总体方差。
  • 小样本(<30),总体符合正态分布时,未知总体方差。

公式解释:

  • 样本均值: x¯
  • 样本数量: n
  • 样本方差: s
  • 总体方差: σ
  • Z分布:z
  • T分布:  t
  • 显著性水平alpha: α

从总体均值的估计程序中,其实就可以为2种情况,是否已知总体方差,已知总体方差就用总体方差,未知总体方差就样本方差。对于是否是正态分布的假设,可以根据中心极限定理的进行人为的预判。对于样本的数量,我们抽保证获得30个以上,在目前的大数据时代,应该没啥难度。

# 进行区间估计
# 参数为:样本,标准差,置信区间
> conf.int=function(x,sigma,conf.level=0.95) {
+   mean<-mean(x)
+   n<-length(x)
+   alpha<-1-conf.level
+   z=qnorm(1-alpha/2,mean=0,sd=1,lower.tail = T)
+   c(mean-sigma*z/sqrt(n),mean+sigma*z/sqrt(n))
+ }

当已知总体标准差时,对总体均值进行区间估计。

> conf.int(price,total_sd,0.95)
[1] 17546.68 21125.07

当未知总体标准差时,用样本方标准差代替,对总体均值进行区间估计。

> conf.int(price,sample_sd,0.95)
[1] 15776.79 22894.96

上面的过程,其实就是Z分布的检验过程,我们也可以直接使用Z分布的检验函数。

# 安装BSDA包
> install.packages('BSDA')
> library(BSDA)

# 当已知总体标准差时,对总体均值进行区间估计。
> z.test(price,sigma.x=total_sd)
	One-sample z-Test
data:  price
z = 21.181, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 17546.68 21125.07
sample estimates:
mean of x 
 19335.87 

# 当未知总体标准差时,用样本方标准差代替,对总体均值进行区间估计。
> z.test(price,sigma.x=sd(price))
	One-sample z-Test
data:  price
z = 10.648, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 15776.79 22894.96
sample estimates:
mean of x 
 19335.87 

我们用z分布检查,得到的计算结果,与手动计算的结果是一致的。

最后,如果我们实现无法获得大于30个样本,而总体分布未知,我们可以使用单总体的T检验进行区间估计。T检验的详细使用说明,请参考文章用R语言解读统计检验-T检验

> t.test(price)
	One Sample t-test
data:  price
t = 10.648, df = 29, p-value = 1.558e-11
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 15621.96 23049.79
sample estimates:
mean of x 
 19335.87 

好的抽样方法,能让我们的样本有更好的代表性。那么好的抽样估计,可以让调查结论更科学、更合理,真正是起到事半功倍的效果。不要再随意拍脑袋做决定了,用统计的专业方法来改进工作吧。如果你想了解如何进行抽样,请参考上一篇文章抽样三步曲之一:R语言进行随机抽样sampling

转载请注明出处:
http://blog.fens.me/r-sampling-estimation

打赏作者

抽样三步曲之一:R语言进行随机抽样sampling

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-sampling

前言

抽样是统计上的基本问题,在实现中生活有非常多的应用场景。美国大选时,媒体进行民意调查,通过对不同地区不同职业的少量人进行询问,就推算出谁的支持率更多。我们进到机场时,会被要求防爆检查,由于人力资源有限,并不是所有人都被检查,而是有一定几率的抽检,来防范重大的安全风险。做机器学习模型的时候,需要训练集和测试集,样本的好坏直接会影响到模型的训练效果。

这些问题其实都可以用统计学的抽样方法来解释和解决。抽样三步曲系列文章:用R语言进行随机抽样用R语言进行抽样估计和用R语言进行抽样校准。

目录

  1. 抽样的基本概念
  2. 简单随机抽样
  3. 分层抽样
  4. 整群抽样
  5. 系统抽样
  6. PPS抽样和Brewer抽样

1. 抽样的基本概念

抽样,是一种可以基于样本的统计信息来获取总体信息,而无需分别调查每个个体信息的方法,其基本要求是要保证所抽取的样本对总体具有充分的代表性。

我们把包含所研究的全部个体的集合看作总体,构成总体的每一个元素作为个体,从总体中抽取一部分的个体所组成的集合叫做样本,样本中的个体数目叫做样本数量。

1.1 为什么需要抽样?

抽样,是为了从样本中得出有关人群的结论,以便通过直接观察群体的一部分(样本)来确定该人群的特征。

  • 比起观察群体中的每个个体,观察样本所需的时间更少,所以观察样本效率更高效。
  • 比起分析整个群体行为,分析样本更简单更方便,所以分析样本难度更小。

抽样会大量的帮助我们节约工作时间,减少工作量,降低分析难度。

1.2 抽样步骤
我们开始进行抽样时,可以按照标准抽样步骤来进行,从而保证抽样过程的合理性。通常来说,会按照下面5个步骤进行抽样。

抽样步骤解释:

  • 第1步:界定总体,在具体抽样前,首先对从总抽取样本的总体范围与界限作明确的界定。比如进行民意调查的抽样,需要考虑18岁以上且有资格的人群。
  • 第2步:制定抽样框,收集总体中全部抽样名单,构成抽样样本的个体或人群的列表。
  • 第3步:决定抽样方法,可以使用概率抽样或非概率抽样,概率抽样方法是从客观上设计抽样规则,每个投票人都具有同样的价值,非概率抽样方法更多是主观上的抽样。
  • 第4步:设置样本数量,样本中要采集的个人或物品的数量,要足够对这一人群做出精准的推断。样本量越大,对这一人群的推断就越准确。
  • 第5步:实际抽取样本,一旦确定了目标人群、抽样框架、抽样技术和样本数量,就可以开始收集数据,进行抽样了。
  • 第6步:评估样本质量,对样本的质量、代表性、偏差等等进行初步的检验和衡量,其目的是防止由于样本的偏差过大而导致的失误

1.3 抽样方法
抽样方法,主要分为概率抽样和非概率抽样。

  • 概率抽样:群体中的每个人都有被选择的平等机会,概率抽样提供了一个真正代表群体的样本,客观性强。
  • 非概率抽样:群体中的每个人都没有被选择的平等机会。因此,可能出现非代表性样本,这种样本无法产生概括性的结果,主观性强。

本文的后面,将重点介绍概率抽样的各种方法。

1.4 有放回和无放回
我们在抽样的时候,需要考虑2种情况,一种是有放回的情况,另一种是无放回的情况。

  • 有放回,指从总体集合中抽出一个元素后,再把这个元素放回总体集合,第二次可以继续抽到这个元素,样本不是唯一的。
  • 无放回,指从总体集合中抽出一个元素后,不把这个元素放回总体集合,第二次可以不能抽到这个元素,样本都是唯一的。

我们假设总体数量为N,样本数量为s,从N个物体中抽取s个样本。

2. 简单随机抽样

简单随机抽样是指,从总体数据中任意抽取指定数量的数据作为样本,每个可能被抽取的样本概率相等。抽取方法分为:放回抽样,不放回抽样。

开发所使用的系统环境

  • Win10 64bit
  • R: 3.6.1 x86_64-w64-mingw32/x64 b4bit

在R语言中,我们可以用基础包的sample()函数进行随机抽样,使用自带的鸢尾花数据集iris。

使用sample()函数,总体数据为N=150个,抽取出样本s=15个。


# 打印数据集iris
> head(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

# 定义总体数和样本数
> N<-nrow(iris) # 总体
> s<-15         # 样本

> # 简单随机抽样,默认为无放回
> s1<-sample(N,s);s1
 [1]  27  33  16  86  43  36  79  98  10  44  82   4  19 145  75

# 打印抽样结果
> iris[s1,]
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
27           5.0         3.4          1.6         0.4     setosa
33           5.2         4.1          1.5         0.1     setosa
16           5.7         4.4          1.5         0.4     setosa
86           6.0         3.4          4.5         1.6 versicolor
43           4.4         3.2          1.3         0.2     setosa
36           5.0         3.2          1.2         0.2     setosa
79           6.0         2.9          4.5         1.5 versicolor
98           6.2         2.9          4.3         1.3 versicolor
10           4.9         3.1          1.5         0.1     setosa
44           5.0         3.5          1.6         0.6     setosa
82           5.5         2.4          3.7         1.0 versicolor
4            4.6         3.1          1.5         0.2     setosa
19           5.7         3.8          1.7         0.3     setosa
145          6.7         3.3          5.7         2.5  virginica
75           6.4         2.9          4.3         1.3 versicolor

sample()函数,默认是不放回的抽样,保证每个结果是唯一的。我们在某些场景下,也会用到放回抽样,在使用sample()函数时,增加一个参数。通过放回抽样的方法,抽出来的结果就会有重复。


> s1a<-sample(N,s,TRUE);s1a # 有放回
 [1]  57  27  57 100  51  93  81 137  46  96   7  15 114  24  34

在R语言中专门有一个包用来进行随机抽样,sampling包。后面,我们的抽样方法,都会基于sampling包进行介绍。sampling包也提供了一个函数srswor()函数,可以进行简单随机抽样。


# sampling包安装
> install.packages("sampling")
> library(sampling)

# 使用srswor()函数,进行随机抽样,无放回抽样
> s2<-srswor(s,N)

# 打印被选中的索引
> s2
  [1] 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 [49] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [97] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1
[145] 0 0 0 0 0 0

# 打印抽样结果
> getdata(iris,s2)
    ID_unit Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1         1          5.1         3.5          1.4         0.2     setosa
8         8          5.0         3.4          1.5         0.2     setosa
10       10          4.9         3.1          1.5         0.1     setosa
20       20          5.1         3.8          1.5         0.3     setosa
22       22          5.1         3.7          1.5         0.4     setosa
28       28          5.2         3.5          1.5         0.2     setosa
48       48          4.6         3.2          1.4         0.2     setosa
54       54          5.5         2.3          4.0         1.3 versicolor
78       78          6.7         3.0          5.0         1.7 versicolor
81       81          5.5         2.4          3.8         1.1 versicolor
82       82          5.5         2.4          3.7         1.0 versicolor
117     117          6.5         3.0          5.5         1.8  virginica
127     127          6.2         2.8          4.8         1.8  virginica
136     136          7.7         3.0          6.1         2.3  virginica
144     144          6.8         3.2          5.9         2.3  virginica

ID_unit列,是新增加的一列,把索引序号直接加在结果数据集了,方便我们肉眼识别随机数是否均匀。

有放回的抽样,可以使用srswr()函数,当有重复选择时,数据会出现大于1的情况,比如第46个值被抽出来2次。


# 有放回抽样
> s2a<-srswr(s,N);s2a 
  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1
 [49] 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0
 [97] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
[145] 0 0 0 0 0 0

# 查看被抽样数据的数量
> table(s2a)
s2a
  0   1   2 
136  13   1 

简单随机抽样,是最简单,也是最常被大家所使用的。

3. 分层抽样

分层抽样是将抽样单位按某种特征或者某种规则划分为不同的层,然后从不同的层中规定的比例,独立、随机的抽取样本,最后将个层的样本结合起来,对总体的目标量估计。分层抽样的优点是样本的代表性比较好,抽样误差比较小。举例说明:对药品进行抽查合格率,从生产环境,运输环节,销售环境进行分层抽样,每个环节按一定的比例抽取药品进行检查。

下面我们使用strata()函数,来模拟分层抽样的操作。取iris鸢尾花的数据集,以Species种属作为分层标准,把数据分成3个层,按照1:3:2的比例进行抽样。


# 判断每个类型的数量
> table(iris$Species)
    setosa versicolor  virginica 
        50         50         50 

# 设置抽样比例为1:3:2
> size <-c(1,3,2)

# 进行分层抽样
> s4=strata(iris,c("Species"),size=size, method="srswor");s4
       Species ID_unit Prob Stratum
47      setosa      47 0.02       1
75  versicolor      75 0.06       2
83  versicolor      83 0.06       2
84  versicolor      84 0.06       2
112  virginica     112 0.04       3
132  virginica     132 0.04       3

# 查看抽样的数据集
> getdata(iris,s4)
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species ID_unit Prob Stratum
47           5.1         3.8          1.6         0.2     setosa      47 0.02       1
75           6.4         2.9          4.3         1.3 versicolor      75 0.06       2
83           5.8         2.7          3.9         1.2 versicolor      83 0.06       2
84           6.0         2.7          5.1         1.6 versicolor      84 0.06       2
112          6.4         2.7          5.3         1.9  virginica     112 0.04       3
132          7.9         3.8          6.4         2.0  virginica     132 0.04       3

结果的数据集增加了3列,ID_unit为索引序号,Prob为被抽样的概率,Stratum为分层的层号。这样就完成了分层抽样。

4. 整群抽样

整群抽样又称聚类抽样,是将总体中各单位归并成若干个互不交叉、互不重复的集合,然后以群为抽样单位抽取样本的一种抽样方式。

整群抽样与分层抽样是容易被混淆的,分层抽样是将抽样单位按某种特征或某种规则划分为不同的层,然后从不同的层中独立,随机地抽取样本。整群抽样则是将总体中若干个单位合并为群,抽样时直接抽取群,然后对选中群中的所有单位全部实施调查。

区别点:分层抽样要求各层之间的差异很大,层内个体或单元差异小,而整群抽样要求群与群之间的差异比较小,群内个体或单元差异大。分层抽样的样本是从每个层内抽取若干单元或个体构成,而整群抽样则是要么整群抽取,要么整群不被抽取。

下面我们使用cluster()函数,来模拟整群抽样的操作。取iris鸢尾花的数据集,以Sepal.Length萼片长度作为分群标准,把数据分成3个群进行抽样。

首先,查看Sepal.Length列数据的分布情况


# 加载工具包
> library(magrittr)
> library(plyr)

# 查看Sepal.Length列数据的分布情况
> Sepal_Length<-table(iris$Sepal.Length) %>% ldply
> names(Sepal_Length)<-c("Sepal.Length","times")
> Sepal_Length
   Sepal.Length times
1           4.3     1
2           4.4     3
3           4.5     1
4           4.6     4
5           4.7     2
6           4.8     5
7           4.9     6
8             5    10
9           5.1     9
10          5.2     4
11          5.3     1
12          5.4     6
13          5.5     7
14          5.6     6
15          5.7     8
16          5.8     7
17          5.9     3
18            6     6
19          6.1     6
20          6.2     4
21          6.3     9
22          6.4     7
23          6.5     5
24          6.6     2
25          6.7     8
26          6.8     3
27          6.9     4
28            7     1
29          7.1     1
30          7.2     3
31          7.3     1
32          7.4     1
33          7.6     1
34          7.7     4
35          7.9     1

以Sepal.Length列进行分群


> s5 <- cluster(iris,clustername="Sepal.Length",size=3,method="srswor",description=TRUE);s5
Number of selected clusters: 3 
Number of units in the population and number of selected units: 150 8 
  Sepal.Length ID_unit       Prob
1          4.5      42 0.08571429
2          6.8     113 0.08571429
3          6.8      77 0.08571429
4          6.8     144 0.08571429
5          6.9     121 0.08571429
6          6.9     140 0.08571429
7          6.9     142 0.08571429
8          6.9      53 0.08571429

# 打印抽样结果
> getdata(iris,s5)
    Sepal.Width Petal.Length Petal.Width    Species Sepal.Length ID_unit       Prob
42          2.3          1.3         0.3     setosa          4.5      42 0.08571429
113         3.0          5.5         2.1  virginica          6.8     113 0.08571429
77          2.8          4.8         1.4 versicolor          6.8      77 0.08571429
144         3.2          5.9         2.3  virginica          6.8     144 0.08571429
121         3.2          5.7         2.3  virginica          6.9     121 0.08571429
140         3.1          5.4         2.1  virginica          6.9     140 0.08571429
142         3.1          5.1         2.3  virginica          6.9     142 0.08571429
53          3.1          4.9         1.5 versicolor          6.9      53 0.08571429

分群的时候,我们定义了抽样3个群,随机抽出了Sepal.Length分别等于4.5,6.8,6.9的值,其他的取值就排除了本次的抽样了。

5. 系统抽样

系统抽样,先将总体的全部单元按照一定顺序排列,采用简单随机抽样抽取第一个样本单元(或称为随机起点),再顺序抽取其余的样本单元,以固定的采样间隔进行选择,这类抽样方法被称为系统抽样。

我们可以使用UPsystematic()函数进行系统抽样的实现。我们以Sepal.Length列为抽样判断的列,先计算出每个总体单元的入样概率,为每个值被选为样本的概率。


# 取样本15个
> s <- 15

# 计算Sepal.Length的值出现概率,入样概率
> pik <-inclusionprobabilities(iris$Sepal.Length,s);pik
> pik
  [1] 0.08727895 0.08385625 0.08043354 0.07872219 0.08556760 0.09241301 0.07872219 0.08556760
  [9] 0.07529949 0.08385625 0.09241301 0.08214489 0.08214489 0.07358813 0.09925841 0.09754706
 [17] 0.09241301 0.08727895 0.09754706 0.08727895 0.09241301 0.08727895 0.07872219 0.08727895
 [25] 0.08214489 0.08556760 0.08556760 0.08899030 0.08899030 0.08043354 0.08214489 0.09241301
 [33] 0.08899030 0.09412436 0.08385625 0.08556760 0.09412436 0.08385625 0.07529949 0.08727895
 [41] 0.08556760 0.07701084 0.07529949 0.08556760 0.08727895 0.08214489 0.08727895 0.07872219
 [49] 0.09070165 0.08556760 0.11979464 0.10952653 0.11808329 0.09412436 0.11123788 0.09754706
 [57] 0.10781517 0.08385625 0.11294923 0.08899030 0.08556760 0.10096977 0.10268112 0.10439247
 [65] 0.09583571 0.11466058 0.09583571 0.09925841 0.10610382 0.09583571 0.10096977 0.10439247
 [73] 0.10781517 0.10439247 0.10952653 0.11294923 0.11637193 0.11466058 0.10268112 0.09754706
 [81] 0.09412436 0.09412436 0.09925841 0.10268112 0.09241301 0.10268112 0.11466058 0.10781517
 [89] 0.09583571 0.09412436 0.09412436 0.10439247 0.09925841 0.08556760 0.09583571 0.09754706
 [97] 0.09754706 0.10610382 0.08727895 0.09754706 0.10781517 0.09925841 0.12150599 0.10781517
[105] 0.11123788 0.13006275 0.08385625 0.12492869 0.11466058 0.12321734 0.11123788 0.10952653
[113] 0.11637193 0.09754706 0.09925841 0.10952653 0.11123788 0.13177410 0.13177410 0.10268112
[121] 0.11808329 0.09583571 0.13177410 0.10781517 0.11466058 0.12321734 0.10610382 0.10439247
[129] 0.10952653 0.12321734 0.12664005 0.13519681 0.10952653 0.10781517 0.10439247 0.13177410
[137] 0.10781517 0.10952653 0.10268112 0.11808329 0.11466058 0.11808329 0.09925841 0.11637193
[145] 0.11466058 0.11466058 0.10781517 0.11123788 0.10610382 0.10096977

# 入参概率合计为样本数
> sum(pik)
[1] 15

使用系统抽样的方法,进行抽样。


# 打印抽样的索引序号
> s6<-UPsystematic(pik);s6
  [1] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
 [49] 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
 [97] 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0
[145] 0 0 0 0 0 0

# 打印抽样的结果
> g6<-getdata(iris,s6);g6
    ID_unit Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
2         2          4.9         3.0          1.4         0.2     setosa
15       15          5.8         4.0          1.2         0.2     setosa
26       26          5.0         3.0          1.6         0.2     setosa
37       37          5.5         3.5          1.3         0.2     setosa
49       49          5.3         3.7          1.5         0.2     setosa
59       59          6.6         2.9          4.6         1.3 versicolor
69       69          6.2         2.2          4.5         1.5 versicolor
78       78          6.7         3.0          5.0         1.7 versicolor
88       88          6.3         2.3          4.4         1.3 versicolor
98       98          6.2         2.9          4.3         1.3 versicolor
108     108          7.3         2.9          6.3         1.8  virginica
117     117          6.5         3.0          5.5         1.8  virginica
125     125          6.7         3.3          5.7         2.1  virginica
134     134          6.3         2.8          5.1         1.5  virginica
143     143          5.8         2.7          5.1         1.9  virginica

系统抽样比简单随机抽样更加方便,但可能存在某种潜在模式,会导致抽样结果是有偏差的。

6. PPS抽样和Brewer抽样

PPS抽样,是按概率比例抽样,在多阶段抽样中,尤其是二阶段抽样中,初级抽样单位被抽中的机率取决于其初级抽样单位的规模大小,初级抽样单位规模越大,被抽中的机会就越大,初级抽样单位规模越小,被抽中的机率就越小。就是将总体按一种准确的标准划分出容量不等的具有相同标志的单位在总体中不同比率分配的样本量进行的抽样。

PPS抽样,总体中含量大的部分被抽中的概率也大,可以提高样本的代表性。PPS抽样的主要优点是使用了辅助信息,减少抽样误差;主要缺点是对辅助信息要求较高,方差的估计较复杂等。PPS抽样为有放回抽样,Brewer抽样为无放回抽样。

我们先进行PPS抽样,使用UPmultinomial()函数,PPS抽样是有放回的,抽样的索引号大于1的值时为,表示被抽到多次。


# 取15个样本,定义样本的入样概率
> s <-15
> pik <-inclusionprobabilities(iris$Sepal.Length,s)

# PPS抽样,打印索引
> s7 <- UPmultinomial(pik);s7
  [1] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
 [49] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0
 [97] 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
[145] 0 0 0 0 0 0

# 查看抽样结果
> getdata(iris,s7)
     Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
11            5.4         3.7          1.5         0.2     setosa
46            4.8         3.0          1.4         0.3     setosa
52            6.4         3.2          4.5         1.5 versicolor
63            6.0         2.2          4.0         1.0 versicolor
72            6.1         2.8          4.0         1.3 versicolor
72.1          6.1         2.8          4.0         1.3 versicolor
74            6.1         2.8          4.7         1.2 versicolor
88            6.3         2.3          4.4         1.3 versicolor
95            5.6         2.7          4.2         1.3 versicolor
104           6.3         2.9          5.6         1.8  virginica
114           5.7         2.5          5.0         2.0  virginica
115           5.8         2.8          5.1         2.4  virginica
118           7.7         3.8          6.7         2.2  virginica
126           7.2         3.2          6.0         1.8  virginica
142           6.9         3.1          5.1         2.3  virginica

Brewer抽样使用UPbrewer()函数,Brewer抽样是无放回的,样本不会被抽到多次。


> # Brewer抽样
> s8 <- UPbrewer(pik);s8
  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
 [49] 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
 [97] 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[145] 0 0 0 0 0 0
> getdata(iris,s8)
    ID_unit Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
19       19          5.7         3.8          1.7         0.3     setosa
22       22          5.1         3.7          1.5         0.4     setosa
31       31          4.8         3.1          1.6         0.2     setosa
41       41          5.0         3.5          1.3         0.3     setosa
49       49          5.3         3.7          1.5         0.2     setosa
59       59          6.6         2.9          4.6         1.3 versicolor
63       63          6.0         2.2          4.0         1.0 versicolor
75       75          6.4         2.9          4.3         1.3 versicolor
89       89          5.6         3.0          4.1         1.3 versicolor
104     104          6.3         2.9          5.6         1.8  virginica
105     105          6.5         3.0          5.8         2.2  virginica
116     116          6.4         3.2          5.3         2.3  virginica
119     119          7.7         2.6          6.9         2.3  virginica
124     124          6.3         2.7          4.9         1.8  virginica
128     128          6.1         3.0          4.9         1.8  virginica

本文对于抽样技术的介绍,算是起了一个头,借助R语言可以非常方便的上手实现。抽样的好坏和实际工作量是有直接的关系,用科学的抽样方法,让抽样更准确。下一篇文章,我将继续介绍抽样估计,科学地得出抽样结论,请参考抽样三步曲之二:用R语言进行抽样估计

转载请注明出处:
http://blog.fens.me/r-sampling

打赏作者