• Archive by category "R语言实践"

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投票支持率
总体 4000 53253.55 37.5%
1次抽样 400 54170.98 39.5%
500次抽样 400*500 53235.58 37.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

打赏作者

R语言中的配置管理yaml

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

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

前言

我们做程序开发,经常需要编写配置文件。早期时,我们多用properites文件,之后用XML,现在很多的应用都用JSON,直到YAML的出场。YAML是专门用来写配置文件的语言,YAML的首要设计目的是为了方便人们读写,而JSON的首要设计目标是简单性和通用性。因此,YAML可以看作是JSON的自然超集,可以提高人们的可读性和更完整的信息模型,每个JSON文件也是一个有效的YAML文件。

目录

  1. YAML是什么
  2. YAML语法介绍
  3. R语言yaml包使用
  4. R语言config包使用
  5. 案例:ETL过程的配置管理

1. YAML是什么

YAML 是 YAML Ain’t Markup Language的缩写,是一种数据序列化语言YAML主页。YAML强调以数据为中心,设计目标包括易于阅读,为不同编程语言提供方便的数据交换,适合描述数据结构,标准接口支持通用工具,支持一次性处理,丰富的表达能力与可扩展性,易于使用。

YAML支持3种数据结构:

  • 键值表,键值对的集合,包括映射,哈希,字典。
  • 序列,为一组排列的值,包括数组,列表。
  • 常量,为单个的不可再分隔的值,包括字符串,布尔值,整数,浮点数,Null,时间,日期

2. YAML语法介绍

由于YAML是JSON的自然超集,所以我们每个YAML语法段,都可以用JSON进行表示。本文介绍 YAML 的语法,你通过在线 JS-YAML 解析器进行验证。

1. YAML文件使用 Unicode 编码作为字符标准编码,默认为UTF-8。
2. 使用“#”来表示注释内容。
3. 使用空格作为嵌套缩进工具,通常建议使用两个空格缩进,不建议使用 tab。
4. YAML文件后缀为 .yml,如:abc.yml 。
5. YAML文件可以由一或多个文档组成,文档间使用—(三个横线)在每文档开始作为分隔符。同时,文档也可以使用…(三个点号)作为结束符。
6. 键值表:使用 :(冒号)空格表示单个键值对,使用”{}”表示一个键值表,”? ” 问号+空格表示复杂的键。当键是一个列表或键值表时,就需要使用本符号来标记。


date: 2015-02-01                                    # 单个键值对
items: {no: 1234, descript: cpu, price: 800.00}     # 多个键值对

JSON输出


{ date: Sun Feb 01 2015 08:00:00 GMT+0800 (中国标准时间),
  items: { no: 1234, descript: 'cpu', price: 800 } }

7. 序列表示:使用 -(横线)单个空格表示单个列表项, 使用”[]”表示一组数据, 组合表示。每个结构都可以嵌套组成复杂的表示结构。

7.1 序列


--- # 序列
- blue                                              # 列表
- red
- green

JSON输出


[ 'blue', 'red', 'green' ]

7.2 嵌套结构


-
  - [blue, red, green]        
  - [Age, Bag]

JSON输出


[ [ [ 'blue', 'red', 'green' ], [ 'Age', 'Bag' ] ] ]

7.3 复合结构对象


--- # 复合结构
languages:
 - English
 - 中文
 - 日本語
websites:
 blog: fens.me
 YAML: yaml.org 
 R: www.r-project.org

JSON输出


{ languages: [ 'English', '中文', '日本語' ],
  websites: { blog: 'fens.me', YAML: 'yaml.org', R: 'www.r-project.org' } }

8. 文本块:使用 “|” 和文本内容缩进表示的块:保留块中已有的回车换行,相当于段落块。使用 “>” 和文本内容缩进表示的块:将块中回车替换为空格,最终连接成一行。使用定界符“”(双引号)、‘’(单引号)或回车表示的块:最终表示成一行。

8.1 段落: 每行保留回车


lines: |          
  北京市朝阳区
  ABC小区#292

JSON输出

{ lines: '北京市朝阳区\nABC小区#292\n' }

8.2 一行:最后一个回车


lines: > 
  北京市朝阳区
  ABC小区#292

JSON输出

{ lines: '北京市朝阳区 ABC小区#292\n' }

8.3 连续一行:无回车


lines:    
  “北京市朝阳区
  ABC小区#292”

JSON输出

{ lines: '“北京市朝阳区 ABC小区#292”' }

9. 数据类型:整数,浮点数,字符串,NULL,日期,布尔,时间。


--- # 数据类型
integer: 12345     # 整数标准形式
octal: 0o34        # 八进制表示,第二个是字母 o
hex: 0xFF          # 十六进制表示
 
float: 1.23e+3     # 浮点数
fixed: 13.67       # 固定小数
minmin: -.inf      # 表示负无穷
notNumber: .NaN    # 无效数字
 
null:              # 空值
boolean: [true, false] # 布尔值
string: '12345'    # 字符串
 
date: 2015-08-23   # 日期
datetime: 2015-08-23T02:02:00.1z  # 日期时间
iso8601: 2015-08-23t21:59:43.10-05:00  # iso8601 日期格式
spaced: 2015-08-23 21:59:43.10 -5      # ?

JSON输出


{ integer: 12345,
  octal: '0o34',
  hex: 255,
  float: 1230,
  fixed: 13.67,
  minmin: -Infinity,
  notNumber: NaN,
  null: null,
  boolean: [ true, false ],
  string: '12345',
  date: Sun Aug 23 2015 08:00:00 GMT+0800 (中国标准时间),
  datetime: '2015-08-23T02:02:00.1z',
  iso8601: Mon Aug 24 2015 10:59:43 GMT+0800 (中国标准时间),
  spaced: Mon Aug 24 2015 10:59:43 GMT+0800 (中国标准时间) }

10. 特殊字符:!(叹号),!!(双叹号)


--- # 特殊字符
isString: !!str 2015-08-23     # 强调是字符串不是日期数据
picture: !!binary |            # Base64  图片
  R0lGODlhDAAMAIQAAP//9/X
  17unp5WZmZgAAAOfn515eXv
  Pz7Y6OjuDg4J+fn5OTk6enp
  56enmleECcgggoBADs=

JSON输出


{ isString: '2015-08-23',
  picture: 
   [ 71,
     73,
     70,
     56,
     57,
     97,
     12,
     0,
     12,
     0,
     132,
     0,
     0,
     255,
     255,
     247,
     245,
     245,
     238,
     233,
     233,
     229,
     102,
     102,
     102,
     0,
     0,
     0,
     231,
     231,
     231,
     94,
     94,
     94,
     243,
     243,
     237,
     142,
     142,
     142,
     224,
     224,
     224,
     159,
     159,
     159,
     147,
     147,
     147,
     167,
     167,
     167,
     158,
     158,
     158,
     105,
     94,
     16,
     39,
     32,
     130,
     10,
     1,
     0,
     59 ] }

11. 下面是内置类型
!!int # 整数类型
!!float # 浮点类型
!!bool # 布尔类型
!!str # 字符串类型
!!binary # 也是字符串类型
!!timestamp # 日期时间类型
!!null # 空值
!!set # 集合
!!omap, !!pairs # 键值列表或对象列表
!!seq # 序列,也是列表
!!map # 键值表

11.1 omap


--- !!omap
- Mark: 65
- Sammy: 63
- Key: 58

JSON输出

[ { Mark: 65 }, { Sammy: 63 }, { Key: 58 } ]

11.2 set


--- !!set           # 注意,“?”表示键为列表,在这里列表为 null
? Mark
? Sammy
? Key

JSON输出

{ Mark: null, Sammy: null, Key: null }

12. 锚点与引用:使用 “&” 定义数据锚点(即要复制的数据),使用 “*” 引用上述锚点数据(即数据的复制目的地)。


hr:
   - Mark McGwire
   # Following node labeled SS
   - &SS 定义要复制的数据
   
rbi:
   - *SS # Subsequent occurrence   这里是数据复制目标
   - Ken Griffey

JSON输出


{ hr: [ 'Mark McGwire', '定义要复制的数据' ],
  rbi: [ '定义要复制的数据', 'Ken Griffey' ] }

YAML 语法结构还是挺简单的,主要是3个数据结构类型的组合使用。我们理解了核心语法后,接下来就可以使用R语言来调用YAML了,进行在R语言中的配置管理。

3. R语言yaml包使用

R语言中,有多个包都可以进行YAML的数据解析,我来介绍2个包一个是yaml包,另一个是config包,下面我们分别介绍一下。2个包都能解析YAML文件,但设计目标是不同的。

  • yaml包,是专门用来解析YAML语法的包,核心功能就是解析YAML文件到R对象,再把R对象生成YAML文件。
  • config包,是专门用来进行配置管理的。用于项目配置管理的场景,当开发,测试和生产环境同一个参数需要配置不同的值,使用config包使用起来比yaml包更简单。

开发所使用的系统环境

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

先让我来学习一下yaml包,yaml包安装和使用,都很容易。首先,我们先进行安装。


> install.packages("yaml")
> library(yaml)

查看yaml包,发现一共就5个函数

  • read_yaml,读取一个YAML文件
  • write_yaml,写入到一个YAML文件
  • as.yaml,转换R对象为YAML字符串
  • yaml.load,转换YAML字符串到R对象
  • yaml.load_file,读文件转换YAML到R对象,看了一个代码和read_yaml()函数没有什么区别

接下来,我们新建一个abc.yml文件,用于描述一个订单。


invoice: 31223
date   : 2020-01-23
bill-to: &id001
    given  : 小明
    address:
        lines: |
            北京市朝阳区
            ABC小区 #292
        city    : 北京
ship-to: *id001
product:
    - sku         : BL394D
      quantity    : 4
      price       : 450.00
    - sku         : BL4438H
      quantity    : 1
      price       : 2392.00
tax  : 251.42
total: 4443.52

使用yaml包,来加载这个文件。

> abc<-read_yaml("abc.yml",fileEncoding = "UTF-8")

查看R语言结构,就是映射为R语言中的list结构。

在R语言中,修改结果。


# 修改日期
abc$date<- '2019-01-01'

# 去掉第二件商品
abc$product[2]<-NA

输出保存为新的文件abc2.yml

write_yaml(abc,"abc2.yml",fileEncoding = "UTF-8")

打开文件abc2.yml


invoice: 31223
date: '2019-01-01'
bill-to:
  given: 小明
  address:
    lines: |
      北京市朝阳区
      ABC小区 #292
    city: 北京
ship-to:
  given: 小明
  address:
    lines: |
      北京市朝阳区
      ABC小区 #292
    city: 北京
product:
- sku: BL394D
  quantity: 4
  price: 450.0
- .na
tax: 251.42
total: 4443.52

我们可以进行JSON转换,再通过JSON的格式进行观察。在R语言中JSON和R的详细介绍,请参考文章R和JSON的傻瓜式编程


# 加载JSON库
> library(RJSONIO)
> abcj<-toJSON(abc)

# 查看JSON
> cat(abcj)
{
 "invoice": 31223,
"date": "2019-01-01",
"bill-to": {
 "given": "小明",
"address": {
 "lines": "北京市朝阳区\nABC小区 #292\n",
"city": "北京" 
} 
},
"ship-to": {
 "given": "小明",
"address": {
 "lines": "北京市朝阳区\nABC小区 #292\n",
"city": "北京" 
} 
},
"product": [
 {
 "sku": "BL394D",
"quantity": 4,
"price":      450 
},
null 
],
"tax":   251.42,
"total":  4443.52 
}

4. R语言config包使用

config包,是RStudio公司开发的,是专门用来进行配置管理的。用于项目配置管理的场景,区别开发,测试和生产环境可能需要不同的值,使用起来比yaml更简单。

首先是安装。


> install.packages("config")
> library(config)

使用config包,我们新一个datasource.yml。用来描述,在软件项目开发时,开发环境(default),测试环境(test),生产环境(production)对于datafile这个参数,会用于不同的文件。


# 开发环境
default:
  time: 5
  datafile: "data-dev.csv"

# 测试环境
test:
  time: 30
  datafile: "data-test.csv"

# 生产环境
production:
  inherits: test
  datafile: "data-prod.csv"

使用config包读取datasource.yml。


# 默认情况会加载default配置
> dev <- config::get(file = "datasource.yml")
> dev
$time
[1] 5

$datafile
[1] "data-dev.csv"

attr(,"config")
[1] "default"
attr(,"file")
[1] "c:\\work\\R\\config\\datasource.yml"

当我们在生产环境时,可以在R的运行时环境配置好全局的环境变量,这个时候就可以直接取对prod的配置信息了。


# 设置全局变量
> Sys.setenv(R_CONFIG_ACTIVE = "production")

# 判断是否production被激活
> config::is_active("production")
[1] TRUE

# 获得配置
> prod <- config::get(file = "datasource.yml")
> prod
$time
[1] 30
$inherits
[1] "test"
$datafile
[1] "data-prod.csv"
attr(,"config")
[1] "production"
attr(,"file")
[1] "c:\\work\\R\\config\\datasource.yml"

所以,我们看到config包,提供了对于yaml的在软件管理的场景的使用方法,很简单,更方便。

5. 案例:ETL过程的配置管理

现在,我们已经完全了解了yaml的语法,以及程序如何调用的使用方法了,那么接下来就是真正地发挥YAML的作用了,把YAML的配置管理用于具体的业务逻辑中。

我设计了一个真实的场景,做数据项目的时候,经常需要对数据进行各种变形处理,就是常说的ETL过程。先加载数据,然后数据处理,最后输出数据,进行可视化。

数据的连续处理过程,在这个思路上其实也出现了一大批软件,通过界面拖拽就实现的功能操作。从软件设计的角度来说,数据连续处理的过程,是非常适合用配置管理的。以配置的方法简化了软件开发的复杂度,又能让用户自己参与其中,用户体验那叫一个爽。

我们就可以用YAML格式,来做为不同软件的核心数据交换结构。基于这个文字型的数据结构,用户可以自己编写,因为易读易懂。同时,基于这个文字型的数据结构,程序可以解析关键的配置信息,这样就打通人和程序的对话格式。把原本看起来很复杂的过程,进行了抽象,便利地实现,用户可参与的配置任务。

我设计一个task.yml的文件,用于描述一个数据处理的过程。


version: v1.0              # 定义版本

resource:                   # 数据源
  - database: &mysql1
      dbname: yaml
      host: localhost
      port: 3306
      username: yaml
      password: yaml
  - file: &file1 
      path: C:/work/demo/test.csv   
      file_encoding: 'utf-8'

data_input:                 # 数据读取
  method: mysql             # 数据读取方法,格式为 method: 方法名
  dataref: *mysql1
  param:                    
    table_name: iris_table
    columns: Sepal.Length

data_analysis:     
  method: mean              # 分析方法,格式为 method: 方法名
  param:                    # 参数
    na_rm: FALSE            # 具体参数格式 参数名: 参数值

data_output:
  method: csv               # 数据输出方法,格式为 method: 方法名
  dataref: *file1

配置解释:

  • version:定义配置文件版本
  • resource:定义数据源,数据库,文件等
  • data_input:定义数据数据输入,引用数据源
  • data_analysis: 定义数据处理方法,计算平均值
  • data_output:定义数据数据输出,引用数据源

运行R语言的代码,就可以执行上面的配置过程,下面的代码是只是一个代码结构,涉及到很多具体的函数调用细节,等下次在详细介绍。


library(magrittr)
task<-read_yaml("task.yml",fileEncoding = "UTF-8")
task

run<-function(task){
  data_input<-function(obj){}          #数据加载
  data_analysis<-function(dat,obj){}   #数据处理
  data_output<-function(dat,obj){}     #数据输出
  
  data_input(task$data_input)  %>% 
    data_analysis(a,task$data_analysis) %>% 
    data_output(task$data_output)
}

# 运行函数
run(task)

本文对YAML进行比较完整的介绍,配置管理工作是软件开发中非常重要的,现在工具越来越便利,把程序员可以从大量复杂的场景剥离出来,只专注于代码本身,可以极大地提高开发效率,并保证程序的质量。把YAML试试用起来吧。

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

打赏作者

2019中国R大会上海 : 面向数据的思维模式和R语言的数据项目开发

跨界知识聚会系列文章,“知识是用来分享和传承的”,各种会议、论坛、沙龙都是分享知识的绝佳场所。我也有幸作为演讲嘉宾参加了一些国内的大型会议,向大家展示我所做的一些成果。从听众到演讲感觉是不一样的,把知识分享出来,你才能收获更多。

关于作者

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://fens.me
  • email: bsspirit@gmail.com

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

前言

本次R语言(上每)大会由上海华东师范大学与统计之都合办,本次会议规模比上一年少了很多,只有2天4个场次,包括主会场,工业应用专场,深度计算专场和R语言数据思维专场。

本次大会,回归数据本质,以落地为主,没有商业广告,更多的是工程实践的思路和方法。结合了学术界、产业界、工业界等,跨领域跨专业的交流,通过碰撞形成了新的火花。很多同仁们都在自己的领域进行创新和尝试,确实是一场不错的盛会。

目录

  1. 我分享的主题:面向数据的思维模式和R语言的数据项目开发
  2. 会议体验和照片分享

1. 我分享的主题:面向数据的思维模式和R语言的数据项目开发

本次大会我被安排在第二天的R语言与数据思维专场,以数据思维为切入点,深入探讨数据思维和toB的数据分析项目实施过程。我的分享的PPT下载,本次活动的官方参会指南

会议日程安排:

我分享主题:面向数据的思维模式和R语言的数据项目开发

目前很多公司机构已经完成了数据的原始积累,如何让沉睡的数据发挥价值,是急需要功课的难关。

我主要为分三个部分进行介绍:

  • 面向数据的思维模式
  • 如何开展一个数据项目
  • R语言项目案例

数据项目和软件项目、互联网项目都有非常大的不同,不确定性、跨学科知识点、工程落地,都是影响数据项目成功与失败的重要因素。掌握数据思维,科学的方法论,专业的团队,便利的工具,才能让数据项目走向成功。

数据分析师每天都有大量的数据需要处理,我们会根据业务的要求做各种复杂的报表,包括了 聚合、分组、排序、筛选、转置、差分、填充、移动、清洗、回归、分布检验、高数计算 等等。有时为了计算一个业务指标,你的SQL怎么写都不会少于10行时。用R语言可以高效地、优雅地解决数据分析中的问题,统计建模、数据挖掘、原型开发,可以节约大量的时间,让我们专业于建模!

2. 会议体验和照片分享

本次大会我体会到的一些关键字:数据思维,工业数据,深度学习,创新和新想法。

2.1 会议体验证和总结

终于最近加班实在是太严重了,没有提前到会场,错过了与很多老朋友的见面机会。第二天的所场次,我都完整的听完了,只是为了赶火车错过了最后一场。

对我来说,最精彩的部分,确实都集中在了第二天的下午场,在我之前的2位嘉宾的分享(黄天元,任坤),确实是精品和干货,而且我们都很享受着R语言编程的过程。我很感慨啊,正是由于我们这些R语言爱好者的兴趣,在一直推动着R语言的落地,非常惊喜!!

我的分享。

能容纳380人的会场。

2.2 相关照片

黄天元:复旦大学在读PHD,R语言数据操纵之美。

任坤,VSCode vs RStudio

张丹,青萌数海CTO,面向数据的思维模型

陈堰平,微软,机器学习模型的可解释性

刘心广,PHD,西门子,汽车车体焊接质量智能预测

我 和 黄天元。

张翔,车轮互联,产品经理眼中灵活管理响应延迟的计算系统

张先秩,澎湃科技CEO,如何做设备端只能化

汤银才,华东师范大学教授,统计学习与机器学习的比较。

年度最后一场,上海R语言大会。最精彩的下午场,都是干货!各嘉宾都在自己的领域为R语言做出贡献,R会越来越强大。感谢华东师范大学和统计之都,所有主办方的小伙们,辛苦了!明年再见。

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

打赏作者

线性和非线性的最小二乘回归

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

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

关于作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

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

前言

很多经典的机器学习模型都是在最小二乘法的基础上发展起来的,理解最小二乘法的原理,对于我们掌握机器学习的高级算法是非常重要的,尤其是回归模型中。最小二乘法是主要是用来做函数拟合,或者求函数极值等问题。

本文以回归模型为例,介绍在最小二乘法在线性和非线性模型中的应用。

目录

  1. 最小二乘法介绍
  2. 最小二乘求解
  3. 线性最小二乘回归
  4. 非线性最小二乘回归

1. 最小二乘法介绍

最小二乘法(Least Squares Method,简记为LSE),源于天文学和测地学上的应用需要,是勒让德( A. M. Legendre)于1805年在其著作《计算慧星轨道的新方法》中提出的,它的主要思想是通过计算理论值与观测值之差的平方和最小,求解未知的参数。

最小二乘法的公式为:

公式解读:

  • H, 目标值,观测值与理论值之差的平方和
  • H’, 最小化的H值
  • y, 观测值
  • yi, 理论值
  • argmin, 最小化函数

我们先了把术语进行统一定义,观测值通常说的是我们的样本数据,理论值就是假设拟合函数产生的预测数据,目标函数是判断一个模型的好不好的标准,进行最小二乘法建立目标函数,让目标函数最小化。

以上图为线性拟合为例,观测值为y1,y2,…,y7,理论值为y1′,y2′,….,y7’,是线上的点并与观测值的x相同,红色线就是我们要做的拟合直线,我们可以画出无数条这样的直线,判断标准是让这条线上的点的理论值与观测值的误差平方和最小,即最小二乘法。

SSE=sum((y1-y1')^2+(y2-y2')^2+(y3-y3')^2+(y4-y4')^2+(y5-y5')^2+(y6-y6')^2+(y7-y7')^2)

利用最小二乘法,可以对数据进行线性和非线性拟合,使误差平方和(SSE)或残差平方和最小。如果观测到的误差近似正态分布时,这种方法是非常有效的。我们在线性回归模型中,要进行残差的正态分布检验,就是基于这个目的的,可以参考文章R语言解读一元线性回归模型

2. 最小二乘求解

在线性方程的求解或是数据曲线拟合中,利用最小二乘法求得的解则被称为最小二乘解。最小二乘法(又称最小平方法)是一种数学优化技术,它通过最小化误差的平方和寻找数据的最佳函数匹配。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

定义数据集x和y,两个变量。

> x<-c(6.19,2.51,7.29,7.01,5.7,2.66,3.98,2.5,9.1,4.2)
> y<-c(5.25,2.83,6.41,6.71,5.1,4.23,5.05,1.98,10.5,6.3)

可视化,画出x,y变量的散点图。

> plot(x,y)

对数据进行最小二乘求解,找到最小二乘解。


> fit<-lsfit(x,y);fit
$coefficients
Intercept         X
0.8310557 0.9004584

$residuals
[1] -1.1548933 -0.2612063 -0.9853975 -0.4332692 -0.8636686  1.0037250  0.6351198 -1.1022017
[9]  1.4747728  1.6870190

$intercept
[1] TRUE

$qr
$qt
[1] -17.1901414   6.2044421  -0.7047339  -0.1530724  -0.5856560   1.2766692   0.9102648
[8]  -0.8295242   1.7584540   1.9625308

$qr
Intercept           X
[1,] -3.1622777 -16.1718880
[2,]  0.3162278   6.8903149
[3,]  0.3162278  -0.2782874
[4,]  0.3162278  -0.2376506
[5,]  0.3162278  -0.0475287
[6,]  0.3162278   0.3936703
[7,]  0.3162278   0.2020970
[8,]  0.3162278   0.4168913
[9,]  0.3162278  -0.5409749
[10,]  0.3162278   0.1701682

$qraux
[1] 1.316228 1.415440

$rank
[1] 2

$pivot
[1] 1 2

$tol
[1] 1e-07

attr(,"class")
[1] "qr"

结果解读:

  • coefficients, 系数的最小二乘估计,包括截距Intercept和自变量X
  • residuals, 残差
  • intercept, 有截距
  • qt, 矩阵QR分解。QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵。关于矩阵的详细操作,请参考文章R语言中的矩阵计算

我们看到lsfit()函数使用矩阵的QR分解的方法,计算得出了最小二乘估计的截距Intercept和自变量系数。

目标:y=a*x+b,用最小二乘法估计出参数a和b。

根据最小二乘法的矩阵计算公式,我直接跳过了推到过程,得出结论。

公式解释:

  • X,为自变量矩阵
  • Y,为因变量矩阵
  • B,为参数,需要求解
  • XT,X的转置矩阵
  • *, 矩阵乘法
  • -1, 计算逆矩阵

手动计算过程:


# 因变量矩阵
> y1<-matrix(y,ncol=1)  

# 自变量矩阵,增加截距列
> x1<-matrix(c(x,rep(1,length(x))),ncol=2)

# 求解最小二乘解
> solve(t(x1)%*%x1) %*% t(x1) %*% y1
          [,1]
[1,] 0.9004584
[2,] 0.8310557

对应到一元线性回归方程中,y = 0.9004584 * x + 0.8310557,与lsfit()的结果一致。

3. 线性最小二乘回归

接下来,我们先进行线性回归分析,以上面的数据集作为样本,只包含x和y。以普通最小二乘法进行线性回归,计算出预测值并和观测值进行比较,计算预测值和观测值的误差平方和,使得误差函数最小。线性回归模型的详细介绍,可以参考文章R语言解读一元线性回归模型, R语言解读多元线性回归模型

建立线性回归模型


# 建立回归模型
> line<-lm(y~x)

# 查看回归模型
> summary(line)
Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1549 -0.9550 -0.3472  0.9116  1.6870 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8311     0.9440   0.880 0.404338    
x             0.9005     0.1698   5.302 0.000726 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.17 on 8 degrees of freedom
Multiple R-squared:  0.7785,	Adjusted R-squared:  0.7508 
F-statistic: 28.12 on 1 and 8 DF,  p-value: 0.0007263

# 残差平方和
> sum(line$residuals^2)
[1] 10.95334

通过查看模型的结果数据,我们可以发现通过T检验自变量x是非常显著,但截距Intercept不显著。通过F检验判断出整个模型的自变量方差是非常显著。通过R^2的相关系数检验可以判断自变量和因变量是不是高度相关的,残差平方和(residual sum-of-squares)为10.95334,效果不好。

接下来,使用模型进行预测,生成预测值。


# 设置x的取值区间
> new <- data.frame(x = seq(min(x),max(x),len = 100))

# 预测结果
> pred<-predict(line,newdata=new)

进行可视化,把观测值与预测值进行展示,以散点表现观测值,以线表示预测值。


> plot(x,y)
> lines(new$x,pred,col="red",lwd=2)

从图中我们可以很容易看出,拟合的效果并不好,与R^2检验不通过是一致的。那么,如果线性模型表现不好,我们可以试试非线性的模型,还是用最小二乘法进行非线性的回归。

4. 非线性最小二乘回归

由于散点的分布并不是线性的,我们再尝试用非线性的函数,进行最小二乘回归的拟合,包括一元二次函数,一元三次函数,指数函数,倒数函数。在R语言中,我们可以使用nls()函数,进行非线性最小二乘回归。

4.1 一元二次函数非线性拟合
定义公式 y ~ a*x+b*x^2+c ,其中x, x^2是一次函数和二次函数,a,b,c分别是参数,我们需要给定参数的初始值。


# 建立回归模型
> nline1 <- nls(y ~ a*x+b*x^2+c, start = list(a = 1,b = 1,c=2))

# 查看模型结果
> nline1
Nonlinear regression model
  model: y ~ a * x + b * x^2 + c
   data: parent.frame()
       a        b        c 
0.008987 0.082188 2.850389 
 residual sum-of-squares: 9.758

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.024e-06

# 分析模型
> summary(nline1)

Formula: y ~ a * x + b * x^2 + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a 0.008987   0.977890   0.009    0.993
b 0.082188   0.088760   0.926    0.385
c 2.850389   2.379763   1.198    0.270

Residual standard error: 1.181 on 7 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.024e-06

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c都不显著,残差平方和(residual sum-of-squares)为9.758,数值过大,效果不好。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


# 模型预测
> npred1<-predict(nline1,newdata=new)

# 可视化
> plot(x,y)
> lines(new$x,npred1,col="blue",lwd=2)

从图中我们可以看出,拟合的是一条曲线,效果不太好。接下来,我们再试试一元三次函数。

4.2 一元三次函数非线性拟合

定义公式 y ~ a*x+b*x^2+c*x^3+d ,其中x, x^2, x^3是一次函数,二次函数和三次函数,a,b,c,d分别是参数,我们需要给定参数的初始值。

使用一元三次函数,进行非线性拟合。


> nline2 <- nls(y ~ a*x+b*x^2+c*x^3+d, start = list(a = 1,b = 1,c=2,d=1))

# 查看模型结果
> nline2
Nonlinear regression model
  model: y ~ a * x + b * x^2 + c * x^3 + d
   data: parent.frame()
      a       b       c       d 
 10.011  -1.822   0.110 -12.516 
 residual sum-of-squares: 3.453

Number of iterations to convergence: 2 
Achieved convergence tolerance: 3.983e-08

# 分析模型
> summary(nline2)

Formula: y ~ a * x + b * x^2 + c * x^3 + d

Parameters:
   Estimate Std. Error t value Pr(>|t|)  
a  10.01051    3.08630   3.244   0.0176 *
b  -1.82152    0.57797  -3.152   0.0198 *
c   0.10998    0.03323   3.310   0.0162 *
d -12.51586    4.88778  -2.561   0.0429 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7586 on 6 degrees of freedom

Number of iterations to convergence: 2 
Achieved convergence tolerance: 3.983e-08

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c,d在0.01的置信区间显著,残差平方和(residual sum-of-squares)为3.453,数值较小,效果应该算是还可以的,能够表达拟合的特征。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred2<-predict(nline2,newdata=new)
> plot(x,y)
> lines(new$x,npred2,col="green",lwd=2)

从图中我们可以看出,拟合的曲线基本穿过观测值的区域,直观上说明效果还不错,与数据检验给出的结论是一致的。

4.3 指数函数非线性拟合

定义公式 y ~ a*exp(b*x)+c ,其中exp(x)是指数函数,a,b,c分别是参数,我们需要给定参数的初始值。


> nline3 <- nls(y ~ a*exp(b*x)+c, start = list(a = 1,b = 1,c=2))

> nline3
Nonlinear regression model
  model: y ~ a * exp(b * x) + c
   data: parent.frame()
     a      b      c 
0.6720 0.2718 2.2087 
 residual sum-of-squares: 8.966

Number of iterations to convergence: 17 
Achieved convergence tolerance: 7.136e-06

> summary(nline3)

Formula: y ~ a * exp(b * x) + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a   0.6720     1.2127   0.554    0.597
b   0.2718     0.1764   1.541    0.167
c   2.2087     2.2447   0.984    0.358

Residual standard error: 1.132 on 7 degrees of freedom

Number of iterations to convergence: 17 
Achieved convergence tolerance: 7.136e-06

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c都不显著,残差平方和(residual sum-of-squares)为8.966,数值过大,效果不好。


> npred3<-predict(nline3,newdata=new)
> plot(x,y)
> lines(new$x,npred3,col="yellow",lwd=2)

从图中我们可以看出,指数函数的拟合曲线,和二次函数表现差不多,效果不太好。

4.4 倒数函数非线性拟合

定义公式 y ~ y ~ b/x+c ,其中1/x是倒数函数,b,c分别是参数,我们需要给定参数的初始值。

使用倒数函数,进行拟合。

> nline4 <- nls(y ~ b/x+c, start = list(b = 1,c=2))

# 查看模型结果
> nline4
Nonlinear regression model
  model: y ~ b/x + c
   data: parent.frame()
    b     c 
-17.0   9.5 
 residual sum-of-squares: 15.74

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.259e-07

> summary(nline4)

Formula: y ~ b/x + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
b  -17.003      4.108  -4.139  0.00326 ** 
c    9.500      1.077   8.817 2.15e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.403 on 8 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.259e-07

通过查看模型的结果数据,我们可以发现通过T检验相关系数b,c非常显著,残差平方和(residual sum-of-squares)为15.74,数值较大,效果不好。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred4<-predict(nline4,newdata=new)
> plot(x,y)
> lines(new$x,npred4,col="gray",lwd=2)


从图中我们可以看出,倒数函数的拟合曲线,确实也效果不太好,而且最后一个点有个想离群点,如果没有最后一个点这个效果应该是更好的。

4.5 模型优化
根据上面的几种模型验证的结果,我手动进行了优化,发现高次函数进行拟合效果会更好。定义公式 y ~ a*x^3+b*x^5+c*x^7+e*d^9+e*x^11 ,其中高次函数定义为11次,并去掉了截距。


> nline5 <- nls(y ~ a*x^3+b*x^5+c*x^7+d*x^9+e*x^11,
+               start = list(a=1,b=2,c=1,d=1,e=1))

# 查看模型结果
> nline5
Nonlinear regression model
  model: y ~ a * x^3 + b * x^5 + c * x^7 + d * x^9 + e * x^11
   data: parent.frame()
         a          b          c          d          e 
 2.959e-01 -2.064e-02  5.837e-04 -7.324e-06  3.367e-08 
 residual sum-of-squares: 2.583

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.343e-07

# 模型分析
> summary(nline5)

Formula: y ~ a * x^3 + b * x^5 + c * x^7 + d * x^9 + e * x^11

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
a  2.959e-01  5.132e-02   5.767   0.0022 **
b -2.064e-02  6.176e-03  -3.341   0.0205 * 
c  5.837e-04  2.468e-04   2.365   0.0644 . 
d -7.324e-06  3.940e-06  -1.859   0.1222   
e  3.367e-08  2.139e-08   1.574   0.1763   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7187 on 5 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.343e-07

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c是显著,残差的误差平方和(residual sum-of-squares)为2.583,数值是最小的,而且比三次函数增长了2上系数,效果应该算不错的,能够表达拟合的特征。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred5<-predict(nline5,newdata=new)
> plot(x,y)
> lines(new$x,npred5,col="purple3",lwd=2)

从图中我们可以看出,拟合的曲线穿过观测值的区域,并且非常贴近部分数据点,直观上说明效果很好,与数据检验给出的结论是一致的。

4.6 总结

我们把上面的1种线性函数和5种非线性的拟合函数的模型评价指标,整理到一个表格中。其中,三次函数和高次函数的拟合效果是不错的,而其他的几种拟合效果,都是不太好的。判断好坏的标准,就是残差的误差平方和,同时参数的T检验显著,才是比较好的拟合模型。

拟合函数 残差平方和 参数T检验显著性 解释性 颜色
一次函数 10.95334 不显著 红色
二次函数 9.758 不显著 蓝色
三次函数 3.453 显著 不好 绿色
指数函数 8.966 不显著 不好 黄色
倒数函数 15.74 显著 灰色
高次函数 2.583 显著 无法解释 紫色

最后,合并所有拟合曲线,画到同一张图上,进行对比。绿色的三次函数拟合曲线,和紫色的高次函数拟合曲线效果最好。

但是,高次函数也带来了的负面问题,一是过拟合,二是业务的可解释性非常差。如果x和y分别带入到业务场景,比如x是货物的重量和y是货物的单价,那么x的11次方,完全是无法理解的,等同于一个黑盒的模型。所以,真正落地项目在使用的时候,要根据业务场景的要求,在准确度,可解释性,复杂度等标准之下,进行权衡。

本文利用最小二乘法,用R语言进行了线性回归和非线性回归的实验。验证了最小二乘法在处理回归问题上是非常有效的,除了回归问题,最小二乘法还能解决分类的问题。最小二乘法作为机器学习的底层算法,是非常值得我们学习和掌握的,明白原理就能掌握一类的高级模型。

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

打赏作者