分层随机抽样中九个非常实用的R函数

Nine Very Practical R Functions in Stratified Random Sampling

摘要: 分层抽样技术是在实际工作中应用得非常广泛的抽样技术之一。但在文献中，还没有方便地可以用于在分层随机抽样中仅给定基本的样本数据时就能解决总体均值和总体比例的点估计和区间估计问题，计算总体均值时样本量的确定及分配问题，计算总体比例时样本量的确定及分配问题，事后分层抽样下总体均值和总体比例的点估计和区间估计等问题通用的R函数。本文自编了九个通用的R函数：Compute_Y_bar_st()、Compute_Y_bar_prop_from_y_bar_h_s_h_st()、Compute_Y_bar_srs_pst()、Compute_P_st()、Compute_P_from_a_h_st()、Compute_P_srs_pst()、Compute_nh_given_n_Y_bar_st()、Compute_n_nh_Y_bar_st()及Compute_n_nh_P_st()，它们将会为需要使用分层抽样技术以提高估计精度进行实际问题分析的使用者提供极大的方便。

Abstract: Stratified sampling technique is one of the sampling techniques widely used in practical work. But in the literature, there are no convenient generic R functions to solve the problem of point estimation and interval estimation of population mean and population proportion, the problem of total sample size and each layer sample size when calculating population mean, the problem of total sample size and each layer sample size when calculating population proportion, and the problem of point estimation and interval estimation of population mean and population proportion in post-stratification sampling, in stratified random sampling when only basic sample data are given. We compile nine generic R functions: Compute_Y_bar_st(), Compute_Y_bar_prop_from_y_bar_h_s_h_st(), Compute_Y_bar_srs_pst(), Compute_P_st(), Compute_P_from_a_h_st(), Compute_P_srs_pst(), Compute_nh_given_n_Y_bar_st(), Compute_n_nh_Y_bar_st(), and Compute_n_nh_P_st(), which will provide great convenience for users who need to use stratified sampling technology to improve the estimation accuracy for practical problem analysis.

1. 引言

分层抽样技术 [1] - [11] 是在实际工作中应用得非常广泛的抽样技术 [12] - [17] 之一。同时R软件 [18] [19] 作为统计学的常用编程工具，它具有完全免费、简洁高效、运行方便等优点。本文选取R软件对分层随机抽样中的总体均值和总体比例的点估计和区间估计问题，计算总体均值时样本量的确定及分配问题，计算总体比例时样本量的确定及分配问题，事后分层抽样下总体均值和总体比例的点估计和区间估计等问题进行了程序实现。针对分层随机抽样，本文自编了九个非常实用的R函数。我们在对这九个R函数进行输入变量及输出变量的解释后给出了相应实际问题的R程序实现。这些内容构成了本文第一作者毕业论文的核心内容 [20]。我们相信，这九个R函数将会为需要使用分层抽样技术以提高估计精度进行实际问题分析的使用者提供极大的方便。

2. 分层随机抽样中九个非常实用的R函数及应用举例

我们推荐分层随机抽样中九个非常实用的R函数。R函数1~3用于解决给定各种样本信息时分层随机抽样下总体均值的点估计和区间估计问题。R函数4~6用于解决给定各种样本信息时分层随机抽样下总体比例的点估计和区间估计问题。R函数7~9用于解决给定各种样本信息时分层随机抽样下总体均值和总体比例的样本量的确定及分配问题。

R函数1：Compute_Y_bar_st()

对于分层随机抽样，给定样本单位的观察值组成的矩阵y_matrix等信息，得到计算总体均值的分层随机抽样的点估计和区间估计的R函数(程序) Compute_Y_bar_st()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中(下载链接：https://pan.baidu.com/s/1y0UvE24vfVm8dTVAUYnskg，提取码：1234)。

下面我们举一个例子来说明该R函数的使用方法。

例1 ( [16] 中例4.1)为调查某地区住户的平均家庭成员数，将该地区分成城市和乡村2层，每层按简单随机抽样抽取10户，调查所获得的数据见表1所示。请估计该地区住户的平均家庭成员数及其95%的置信区间。

Table 1. Household membership survey data

表1. 家庭成员数调查数据

解：对于分层随机抽样，由理论公式，可以计算：

$t=Z_{\alpha /2}\approx 1.959964,N={\displaystyle \underset{h=1}{\overset{2}{\sum }}{N}_h}=750$

${\bar{y}}_{st}={\displaystyle \underset{h=1}{\overset{2}{\sum }}{W}_h{\bar{y}}_h}\approx \text{3}\text{.733333},v\left({\bar{y}}_{st}\right)={\displaystyle \underset{h=1}{\overset{2}{\sum }}{W}_h^2\left(1-f_h\right)\frac{s_h^2}{n_h}}\approx \text{0}\text{.06708148}$

$se\left({\bar{y}}_{st}\right)=\sqrt{v\left({\bar{y}}_{st}\right)}\approx \text{0}\text{.2590009},\Delta =t\cdot se\left({\bar{y}}_{st}\right)\approx \text{0}\text{.5076325},\gamma =\frac{\Delta }{{\bar{y}}_{st}}\approx \text{0}\text{.135973}$

$L_{\bar{Y}}={\bar{y}}_{st}-\Delta \approx \text{3}\text{.225701,}{U}_{\bar{Y}}={\bar{y}}_{st}+\Delta \approx \text{4}\text{.240966}$

代入数据，调用R函数Compute_Y_bar_st()进行计算，详细的R程序输入及输出结果请见附录A.1。

因此，估计得该地区住户的平均家庭成员数为3.733人，抽样标准误为0.259人，平均家庭成员数的95%置信区间为(3.226，4.241)人。

R函数2：Compute_Y_bar_prop_from_y_bar_h_s_h_st()

对于按比例分配的分层随机抽样，给定各层的样本均值y_bar_h和各层的样本标准差s_h等信息，得到计算总体均值的分层随机抽样的点估计和区间估计的R函数(程序) Compute_Y_bar_prop_from_y_bar_ h_s_h_st()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中。

下面我们举一个例子来说明该R函数的使用方法。

例2 ( [16] 中练习4.5)某开发区利用电话调查对区内居民消费冷冻食品情况进行调查，将电话号码(6位数字)的前2位作为一部分，后4位作为一部分，前2位代表局号，局号及每个局号中拥有的电话数可以找到，按局号分层，按每个局号(去掉商户后)拥有的电话数比例分配样本量(各层抽样比可以忽略)。调查后各层样本户购买冷冻食品支出的中间结果见表2所示。试估计该开发区居民户购买冷冻食品的平均支出及估计的95%置信区间。

Table 2. Sample households purchase frozen food expenditure

表2. 样本户购买冷冻食品支出

解：对于按比例分配的分层随机抽样，由理论公式，可以计算：

$t=Z_{\alpha /2}\approx \text{1}\text{.959964},n={\displaystyle \underset{h=1}{\overset{10}{\sum }}{n}_h}=\text{200,}{\bar{y}}_{st}={\displaystyle \underset{h=1}{\overset{10}{\sum }}{W}_h{\bar{y}}_h}\approx 75.792$

$v\left({\bar{y}}_{st}\right)=\frac{1-f}{n}{\displaystyle \underset{h=1}{\overset{10}{\sum }}{W}_h{s}_h^2}\approx 59.46035,se\left({\bar{y}}_{st}\right)=\sqrt{v\left({\bar{y}}_{st}\right)}\approx 7.711054$

$L_{\bar{Y}}={\bar{y}}_{st}-t\cdot se\left({\bar{y}}_{st}\right)\approx \text{60}\text{.67861,}{U}_{\bar{Y}}={\bar{y}}_{st}+t\cdot se\left({\bar{y}}_{st}\right)\approx \text{90}\text{.90539}$

代入数据，调用R函数Compute_Y_bar_prop_from_y_bar_h_s_h_st()进行计算，详细的R程序输入及输出结果请见附录A.2。

故该开发区居民户购买冷冻食品的平均支出为75.792元，标准误差为7.711元，其95%置信区间为(60.679，90.905)元。

R函数3：Compute_Y_bar_srs_pst()

对于事后分层抽样，给定各层的样本均值y_bar_h和各层的样本标准差s_h等信息，得到计算总体均值的事后分层抽样和简单随机抽样的点估计和区间估计的R函数(程序) Compute_Y_bar_srs_pst()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中。

下面我们举一个例子来说明该R函数的使用方法。

例3 ( [16] 中例4.6)某高校欲了解在校学生用于课外进修(如各种考证辅导班、外语辅导班等)的开支，在全校8000名学生中抽出了一个200人的简单随机样本。根据学生科的统计，本科生人数为全校学生的70%，调查最近一个学期课外进修支出的结果见表3所示，试估计全校学生用于课外进修的平均开支。

Table 3. Results of the survey on expenditure on after-school education of current students

表3. 在校学生课外进修开支调查结果

解：由理论公式，可以计算：

$t=Z_{\alpha /2}\approx \text{1}\text{.959964},n={\displaystyle \underset{h=1}{\overset{2}{\sum }}{n}_h}=\text{200},f=\frac{n}{N}=0.025,\bar{y}=\frac{1}{n}{\displaystyle \underset{h=1}{\overset{2}{\sum }}{n}_h{\bar{y}}_h}=283.8$

$s^2={\displaystyle \underset{h=1}{\overset{2}{\sum }}\frac{n_h-1}{n-1}{s}_h^2}+{\displaystyle \underset{h=1}{\overset{2}{\sum }}\frac{n_h}{n-1}{\left({\bar{y}}_h-\bar{y}\right)}^2}\approx 86772.05,s=\sqrt{s^2}\approx 294.571$

采用事后分层估计，则

$v\left({\bar{y}}_{pst}\right)=\frac{1-f}{n}{\displaystyle \underset{h=1}{\overset{2}{\sum }}{W}_h{s}_h^2}+\frac{1}{n^2}{\displaystyle \underset{h=1}{\overset{2}{\sum }}\left(1-W_h\right)s_h^2}\approx 381.8343$

${\bar{y}}_{pst}={\displaystyle \underset{h=1}{\overset{2}{\sum }}{W}_h{\bar{y}}_h}=276.2,se\left({\bar{y}}_{pst}\right)=\sqrt{v\left({\bar{y}}_{pst}\right)}\approx 19.54058$

$L_{{\bar{Y}}_{pst}}={\bar{y}}_{pst}-t\cdot se\left({\bar{y}}_{pst}\right)\approx \text{237}\text{.9012,}{U}_{{\bar{Y}}_{pst}}={\bar{y}}_{pst}+t\cdot se\left({\bar{y}}_{pst}\right)\approx \text{314}\text{.4988}$

采用简单随机估计，则

$v\left(\bar{y}\right)=\frac{1-f}{n}{s}^2\approx \text{423}\text{.0137,}se\left(\bar{y}\right)=\sqrt{v\left(\bar{y}\right)}\approx \text{20}\text{.5673}$

$L_{{\bar{Y}}_{srs}}=\bar{y}-t\cdot se\left(\bar{y}\right)\approx \text{243}\text{.4888,}{U}_{{\bar{Y}}_{srs}}=\bar{y}+t\cdot se\left(\bar{y}\right)\approx \text{324}\text{.1112}$

代入数据，调用R函数Compute_Y_bar_srs_pst()进行计算，详细的R程序输入及输出结果请见附录A.3。

因此，采用事后分层估计的平均开支为276.2元，其抽样标准误为19.541元，总体均值的事后分层抽样的95%置信区间为(237.901，314.499)元。而采用简单随机估计的平均开支为283.8元，其抽样标准误为20.567元，总体均值的简单随机抽样的95%置信区间为(243.489，324.111)元。事后分层估计的抽样标准误比简单随机估计的抽样标准误要小。

R函数4：Compute_P_st()

对于分层随机抽样，给定样本单位的观察值组成的矩阵y_matrix等信息，得到计算总体比例的分层随机抽样的点估计和区间估计的R函数(程序) Compute_P_st()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中。

下面我们举一个例子来说明该R函数的使用方法。

例4 ( [16] 中例4.2)对某地区的居民拥有家庭电脑的情况进行调查，以居民户为抽样单位，根据收人水平将居民户划分为4层，每层按简单随机抽样抽取10户，调查获得数据见表4所示。估计该地区居民拥有家庭电脑的比例及抽样标准误。

Table 4. The sample households have a home computer

表4. 样本户拥有家庭电脑情况

解：对于分层随机抽样，由理论公式，可以计算：

$t=Z_{\alpha /2}\approx \text{1}\text{.959964},N={\displaystyle \underset{h=1}{\overset{\text{4}}{\sum }}{N}_h}=\text{2850,}{p}_{st}={\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h{p}_h}=0.2$

$v\left(p_{st}\right)={\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h^{\text{2}}v\left(p_h\right)}\approx \text{0}\text{.004998324},se\left(p_{st}\right)=\sqrt{v\left(p_{st}\right)}\approx \text{0}\text{.07069883}$

$\Delta =t\cdot se\left(p_{st}\right)\approx \text{0}\text{.1385672},\gamma =\frac{\Delta }{p_{st}}\approx \text{0}\text{.6928358}$

$L_P=p_{st}-\Delta \approx \text{0}\text{.06143284,}{U}_P=p_{st}+\Delta \approx \text{0}\text{.3385672}$

代入数据，调用R函数Compute_P_st()进行计算，详细的R程序输入及输出结果请见附录A.4。

因此，估计得该地区居民拥有家庭电脑的比例为0.2，抽样标准误为0.071，总体比例的分层估计的95%置信区间为(0.061，0.339)。

R函数5：Compute_P_from_a_h_st()

对于分层随机抽样，给定各层样本中具有所考虑特征的单位数a_h等信息，得到计算总体比例的分层随机抽样的点估计和区间估计的R函数(程序)Compute_P_from_a_h_st()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中。

下面我们举一个例子来说明该R函数的使用方法。

例5 ( [16] 中练习4.4)随着经济发展，某市居民正在悄悄改变过年的习惯，虽然仍有大多数居民除夕夜在家吃年夜饭、看电视节目，但也有些家庭到饭店吃年夜饭，或逛夜市，或利用过年的假期到外地旅游。为研究这种现象，某研究机构以市中心165万居民户作为研究对象，将居民户按6个行政区分层，在每个行政区随机抽出30户进行了调查(各层抽样比可以忽略)。每个行政区的情况以及在家吃年夜饭、看电视节目的居民户比例见表5。试估计该市居民在家吃年夜饭的比例，并给出抽样标准误。

Table 5. Proportion of households and number of households

表5. 居民户比例及在家居民户数

解：对于分层随机抽样，由理论公式，可以计算：

$t=Z_{\alpha /2}\approx \text{1}\text{.959964},p_{st}={\displaystyle \underset{h=1}{\overset{6}{\sum }}{W}_h{p}_h}\approx 0.924$

$v\left(p_{st}\right)={\displaystyle \underset{h=1}{\overset{6}{\sum }}{W}_h^2\left(1-f_h\right)\frac{p_h{q}_h}{n_h-1}}\approx \text{0}\text{.0003969808},se\left(p_{st}\right)=\sqrt{v\left(p_{st}\right)}\approx \text{0}\text{.01992438}$

$L_P=p_{st}-t\cdot se\left(p_{st}\right)\approx \text{0}\text{.8849489,}{U}_P=p_{st}+t\cdot se\left(p_{st}\right)\approx \text{0}\text{.9630511}$

$p_h=\frac{a_h}{n_h},即p_1=0.900,p_2\approx 0.933,p_3=0.900,p_4\approx 0.867,p_5\approx 0.933,p_6\approx 0.967$

$q_h=1-p_h,即q_1=0.100,q_2\approx 0.067,q_3=0.100,q_4\approx 0.133,q_5\approx 0.067,q_6\approx 0.033$

$W_1{p}_1\approx 0.162,W_2{p}_2\approx 0.196,W_3{p}_3\approx 0.126,W_4{p}_4\approx 0.078,W_5{p}_5\approx 0.149,W_6{p}_6\approx 0.213$

代入数据，调用R函数Compute_P_from_a_h_st()进行计算，详细的R程序输入及输出结果请见附录A.5。

故该市居民在家吃年夜饭的比例为0.924，抽样标准误为0.020，95%置信区间为(0.885，0.963)。

R函数6：Compute_P_srs_pst()

对于事后分层抽样，给定各层样本中具有所考虑特征的单位数a_h等信息，得到计算总体比例的简单随机抽样和事后分层抽样的点估计和区间估计的R函数(程序) Compute_P_srs_pst()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中。

下面我们举一个例子来说明该R函数的使用方法。

例6 ( [16] 中练习4.9)某公司进行财务审计，需要对原始凭证进行审核，该公司先后有两名出纳，由A出纳登记的原始凭证占70%，B出纳登记的原始凭证占30%。审计人员从原始凭证中随机抽出100份，结果发现，由A、B出纳登记的原始凭证分别为43份和57份，差错分别为1份和2份。

1) 用简单随机抽样的公式估计登记原始凭证的差错率，并估计抽样标准误；

2) 用事后分层的公式估计登记原始凭证的差错率，并估计抽样标准误(有限总体校正系数 $1-f\approx 1$ )。

解：1) 对于简单随机抽样，由理论公式，可以计算：

$t=Z_{\alpha /2}\approx \text{1}\text{.959964},a={\displaystyle \underset{h=1}{\overset{2}{\sum }}{a}_h}=3,n={\displaystyle \underset{h=1}{\overset{2}{\sum }}{n}_h}=100$

$p=\frac{a}{n}=0.03,q=1-p=0.97$

$v\left(p\right)=\frac{1-f}{n-1}pq\approx \frac{pq}{n-1}\approx \text{0}\text{.0002939394},se\left(p\right)=\sqrt{v\left(p\right)}\approx \text{0}\text{.01714466}$

$L_{P_{srs}}=p-t\cdot se\left(p\right)\approx -0.003602918\text{,}{U}_{P_{srs}}=p+t\cdot se\left(p\right)\approx \text{0}\text{.06360292}$

2) 对于事后分层抽样，由理论公式，可以计算：

$p_{pst}={\displaystyle \underset{h}{\sum }{W}_h{p}_h\approx \text{0}\text{.02680539}}$

$v_{\text{1}}\left(p_{pst}\right)\approx \frac{1-f}{n}{\displaystyle \underset{h}{\sum }{W}_h{s}_h^2}+\frac{1}{n^2}{\displaystyle \underset{h}{\sum }\left(1-W_h\right)s_h^2}\approx \text{0}\text{.0002692841},s{e}_{\text{1}}\left(p_{pst}\right)=\sqrt{v_{\text{1}}\left(p_{pst}\right)}\approx \text{0}\text{.01640988}$

$L_{P_{pst,1}}=p_{pst}-t\cdot s{e}_{\text{1}}\left(p_{pst}\right)\approx -\text{0}\text{.005357385,}{U}_{P_{pst,1}}=p_{pst}+t\cdot s{e}_{\text{1}}\left(p_{pst}\right)\approx \text{0}\text{.05896816}$

$v_{\text{2}}\left(p_{pst}\right)={\displaystyle \underset{h}{\sum }{W}_h^2\left(1-f_h\right)\frac{p_h{q}_h}{n_h-1}}\approx \text{0}\text{.0003194205},s{e}_{\text{2}}\left(p_{pst}\right)=\sqrt{v_{\text{2}}\left(p_{pst}\right)}\approx \text{0}\text{.01787234}$

$L_{P_{pst,\text{2}}}=p_{pst}-t\cdot s{e}_{\text{2}}\left(p_{pst}\right)\approx -\text{0}\text{.008223753,}{U}_{P_{pst,\text{2}}}=p_{pst}+t\cdot s{e}_{\text{2}}\left(p_{pst}\right)\approx \text{0}\text{.06183452}$

代入数据，调用R函数Compute_P_srs_pst()进行计算，详细的R程序输入及输出结果请见附录A.6。

故用简单随机抽样估计的差错率为0.030，抽样标准误为0.0171，95%置信区间为(−0.004，0.064)。用事后分层估计的差错率为0.027，抽样标准误的第一种估计为0.0164 < 0.0171，95%置信区间的第一种估计为(−0.005，0.059)，抽样标准误的第二种估计为0.0179 > 0.0171，95%置信区间的第二种估计为(−0.008，0.062)。

R函数7：Compute_nh_given_n_Y_bar_st()

对于按比例分配和尼曼分配的分层随机抽样，给定样本单位的观察值组成的矩阵y_matrix和样本量n等信息，得到计算总体均值时所需的各层样本量的R函数(程序) Compute_nh_given_n_Y_bar_st()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中。

下面我们举一个例子来说明该R函数的使用方法。

例7 ( [16] 中例4.3)对某地区的2850户居民豆制品年消费支出进行调查，以居民户为抽样单位，根据收人水平将居民户划分为4层，每层按简单随机抽样抽取10户，调查获得以下数据，见表6所示。样本量为 $n=40$ ，按比例分配和尼曼分配时，各层的样本量分别应为多少？

Table 6. Annual consumption expenditure table of sample household soybean products

表6. 样本户豆制品年消费支出表

解：对于分层随机抽样，由理论公式，可以计算：

$L=4,N={\displaystyle \underset{h=1}{\overset{\text{4}}{\sum }}{N}_h}=\text{2850},{\bar{y}}_{st}={\displaystyle \underset{h=1}{\overset{\text{4}}{\sum }}{W}_h{\bar{y}}_h}\approx \text{78}\text{.77193},{\displaystyle \underset{h=1}{\overset{\text{4}}{\sum }}{W}_h{s}_h}\approx \text{40}\text{.61926}$

$W_h=\frac{N_h}{N},即W_1\approx 0.07018,W_2\approx 0.14035,W_3\approx 0.26316,W_4\approx 0.52632$

${\bar{y}}_h=\frac{1}{n_h}{\displaystyle \underset{i=1}{\overset{n_h}{\sum }}{y}_{hi}},即{\bar{y}}_1=43.5,{\bar{y}}_2=112.0,{\bar{y}}_3=180.0,{\bar{y}}_4=24.0$

$s_h^2=\frac{1}{n_h-1}{\displaystyle \underset{i=1}{\overset{n_h}{\sum }}{\left(y_{hi}-{\bar{y}}_h\right)}^2},即s_1^2\approx 1433.61,s_2^2\approx 1956.6\text{7},s_3^2\approx 8622.22,s_4^2\approx 193.33$

$W_1{s}_1\approx 2.6571,W_2{s}_2\approx 6.2083,W_3{s}_3\approx 24.4358,W_4{s}_4\approx 7.3181$

按比例分配的样本量为 $n_h=n{W}_h$ ，即

$n_{\text{1}}\approx 2.81,n_{\text{2}}\approx 5.61,n_{\text{3}}\approx 10.53,n_{\text{4}}\approx 21.05$

按尼曼分配的样本量为

$n_h=n\frac{W_h{s}_h}{{\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h{s}_h}}$

即

$n_{\text{1}}\approx 2.62,n_{\text{2}}\approx 6.11,n_{\text{3}}\approx 24.06,n_{\text{4}}\approx 7.21$

代入数据，调用R函数Compute_nh_given_n_Y_bar_st()进行计算，详细的R程序输入及输出结果请见附录A.7。

因此，按比例分配确定的各层样本量为3，6，10，21；按尼曼分配确定的各层样本量为3，6，24，7。

R函数8：Compute_n_nh_Y_bar_st()

对于按比例分配和尼曼分配的分层随机抽样，给定样本单位的观察值组成的矩阵y_matrix及相应的精度要求，得到计算总体均值时所需的总样本量及各层样本量的R函数(程序) Compute_n_nh_Y_bar_st()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中。

下面我们举一个例子来说明该R函数的使用方法。

例8 ( [16] 中例4.4和例4.3)在例7中，如果要求在95%置信度下，相对允许误差不超过10%，则按比例分配和尼曼分配时，总样本量及各层样本量分别为多少？

解：对于分层随机抽样，由理论公式，可以计算：

$t=Z_{\alpha /2}\approx \text{1}\text{.959964},L=4,N={\displaystyle \underset{h=1}{\overset{\text{4}}{\sum }}{N}_h}=\text{2850},{\bar{y}}_{st}={\displaystyle \underset{h=1}{\overset{\text{4}}{\sum }}{W}_h{\bar{y}}_h}\approx \text{78}\text{.77193}$

${\displaystyle \underset{h=1}{\overset{\text{4}}{\sum }}{W}_h{s}_h^2}\approx \text{2745}\text{.984,}{\displaystyle \underset{h=1}{\overset{\text{4}}{\sum }}{W}_h{s}_h}\approx \text{40}\text{.61926,}V={\left(\frac{\gamma {\bar{y}}_{st}}{t}\right)}^2\approx \text{16}\text{.15276}$

$W_h=\frac{N_h}{N},即W_1\approx 0.07018,W_2\approx 0.14035,W_3\approx 0.26316,W_4\approx 0.52632$

${\bar{y}}_h=\frac{1}{n_h}{\displaystyle \underset{i=1}{\overset{n_h}{\sum }}{y}_{hi}},即{\bar{y}}_1=43.5,{\bar{y}}_2=112.0,{\bar{y}}_3=180.0,{\bar{y}}_4=24.0$

$s_h^2=\frac{1}{n_h-1}{\displaystyle \underset{i=1}{\overset{n_h}{\sum }}{\left(y_{hi}-{\bar{y}}_h\right)}^2},即s_1^2\approx 1433.61,s_2^2\approx 1956.67,s_3^2\approx 8622.22,s_4^2\approx 193.33$

$W_1{s}_1^2\approx 100.604,W_2{s}_2^2\approx 274.620,W_3{s}_3^2\approx 2269.006,W_4{s}_4^2\approx 101.754$

$W_1{s}_1\approx 2.657,W_2{s}_2\approx 6.208,W_3{s}_3\approx 24.436,W_4{s}_4\approx 7.318$

按比例分配确定的

$n_0=\frac{{\displaystyle \sum W_h{s}_h^2}}{V}\approx \text{170}\text{.00}$

对 $n_0$ 进行修正，得到修正后的样本量为

$n=\frac{n_0}{1+\frac{n_0}{N}}\approx 160.43$

且各层应分配的样本量为 $n_h=n{W}_h$ ，即

$n_1\approx 11.30,n_2\approx 22.60,n_3\approx 42.37,n_4\approx 84.74$

按尼曼分配确定的总样本量为

$n=\frac{{\left({\displaystyle \sum W_h{s}_h}\right)}^2}{V+\frac{{\displaystyle \sum W_h{s}_h^2}}{N}}\approx 96.40$

且各层样本量为

$n_h=n\frac{W_h{s}_h}{{\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h{s}_h}}$

即

$n_1\approx 6.35,n_2\approx 14.83,n_3\approx 58.35,n_4\approx 17.48$

代入数据，调用R函数Compute_n_nh_Y_bar_st()进行计算，详细的R程序输入及输出结果请见附录A.8。

因此，按比例分配确定的总样本量为161，各层样本量为11，23，42，85；按尼曼分配确定的总样本量为97，各层样本量为6，15，58，18 (程序结果为17，但由于总样本量需达到97，故该层样本量需加1成18)。

R函数9：Compute_n_nh_P_st()

对于按比例分配和尼曼分配的分层随机抽样，给定各层的总体单位数N_h，各层的样本比例p_h及相应的精度要求，得到计算总体比例时所需的总样本量及各层样本量的R函数(程序) Compute_n_nh_P_st()。由于正文版面的限制，该R函数的内容及输入输出的解释放在了补充材料中。

下面我们举一个例子来说明该R函数的使用方法。

例9 ( [16] 中例4.5和例4.3)在例4中，如果要求在95%的置信度下，绝对误差不超过5%，则按比例分配和尼曼分配时，总样本量及各层样本量分别为多少？

解：对于分层随机抽样，由理论公式，可以计算：

$t=Z_{\alpha /2}\approx 1.959964,L=4,N={\displaystyle \underset{h=1}{\overset{4}{\sum }}{N}_h}=2850,p_{st}={\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h{p}_h}=0.2$

$V={\left(\frac{\Delta }{t}\right)}^2\approx 0.0006507944,{\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h{p}_h{q}_h}\approx 0.1442105,{\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h\sqrt{p_h{q}_h}}\approx 0.3710258$

$W_h=\frac{N_h}{N},即W_1\approx 0.07018,W_2\approx 0.14035,W_3\approx 0.26316,W_4\approx 0.52632$

$q_h=1-p_h,即q_1=0.8,q_2=0.8,q_3=0.6,q_4=0.9$

$W_1{p}_1\approx 0.014,W_2{p}_2\approx 0.028,W_3{p}_3\approx 0.105,W_4{p}_4\approx 0.053$

$W_1{p}_1{q}_1\approx 0.011,W_2{p}_2{q}_2\approx 0.022,W_3{p}_3{q}_3\approx 0.063,W_4{p}_4{q}_4\approx 0.047$

$W_1\sqrt{p_1{q}_1}\approx 0.028,W_2\sqrt{p_2{q}_2}\approx 0.056,W_3\sqrt{p_3{q}_3}\approx 0.129,W_4\sqrt{p_4{q}_4}\approx 0.158$

按比例分配确定的

$n_0=\frac{{\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h{p}_h{q}_h}}{V}\approx 221.59$

对 $n_0$ 进行修正，得到修正后的样本量为

$n=\frac{n_0}{1+\frac{n_0}{N}}\approx 205.61$

且各层样本量为 $n_h=n{W}_h$ ，即

$n_1\approx 14.46,n_2\approx 28.91,n_3\approx 54.21,n_4\approx 108.42$

按尼曼分配确定的总样本量为

$n=\frac{{\left({\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h\sqrt{p_h{q}_h}}\right)}^2}{V+\frac{{\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h{p}_h{q}_h}}{N}}\approx 196.27$

且各层样本量为

$n_h=n\frac{W_h\sqrt{p_h{q}_h}}{{\displaystyle \underset{h=1}{\overset{4}{\sum }}{W}_h\sqrt{p_h{q}_h}}}$

即

$n_1\approx 14.90,n_2\approx 29.81,n_3\approx 68.45,n_4\approx 83.84$

代入数据，调用R函数Compute_n_nh_P_st()进行计算，详细的R程序输入及输出结果请见附录A.9。

因此，按比例分配所需的总样本量为206，各层样本量为15 (程序结果为14，但由于总样本量需达到206，故该层样本量需加1成15)，29，54，108；按尼曼分配所需的总样本量为197，各层样本量为15，30，68，84。

3. 总结

本文就分层随机抽样的R软件实现方面自编了九个非常实用的R函数，分别是Compute_Y_bar_st() (用于分层随机抽样下总体均值的点估计和区间估计)、Compute_Y_bar_prop_from_y_bar_h_s_h_st() (用于给定各层的样本均值和各层的样本标准差等信息的按比例分配的分层随机抽样下总体均值的点估计和区间估计)、Compute_Y_bar_srs_pst() (用于事后分层抽样和简单随机抽样下总体均值的点估计和区间估计)、Compute_P_st() (用于分层随机抽样下总体比例的点估计和区间估计)、Compute_P_from_a_h_st() (用于给定各层样本中具有所考虑特征的单位数等信息的分层随机抽样下总体比例的点估计和区间估计)、Compute_P_srs_pst() (用于事后分层抽样和简单随机抽样下总体比例的点估计和区间估计)、Compute_ nh_given_n_Y_bar_st() (用于给定总样本量等信息的按比例分配和尼曼分配的分层随机抽样下计算总体均值时所需的各层样本量)、Compute_n_nh_Y_bar_st() (用于按比例分配和尼曼分配的分层随机抽样下计算总体均值时所需的总样本量及各层样本量)及Compute_n_nh_P_st() (用于按比例分配和尼曼分配的分层随机抽样下计算总体比例时所需的总样本量及各层样本量)。我们相信，这九个R函数一定可以给利用分层随机抽样以提高估计精度进行实际问题分析的使用者提供极大的方便。

基金项目

教育部人文社会科学研究西部和边疆地区项目：基于临床试验大数据的条件势的贝叶斯无效分析的基础研究(20XJC910001)，2020.1~2022.12。国家社科基金西部项目：基于贝叶斯的八种预测势在临床试验中用于节约新药研发成本的评价研究(21XTJ001)，2021.9~2024.12。国家自然科学基金面上项目：大数据驱动的中小微企业全息风险评估与介观调控机制研究(72071019)，2021.1~2024.12。

附录

本附录给出了九个例子的详细的R程序输入及输出的结果。

A.1. 例1的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

>n_h = c(10, 10)

>N_h = c(250, 500)

> alpha = 0.05

> y1 = c(3, 2, 3, 4, 3, 3, 4, 5, 2, 3); y1

> y2 = c(3, 4, 5, 5, 4, 3, 6, 2, 4, 4); y2

>y_matrix = rbind(y1, y2); y_matrix

y1 3 2 3 4 3 3 4 5 2 3

y2 3 4 5 5 4 3 6 2 4 4

>

> result = Compute_Y_bar_st(y_matrix, n_h, N_h, alpha); result

$df_0

tN y_bar_stv_y_bar_stse_y_bar_st Deltagamma

1 1.959964 750 3.733333 0.06708148 0.2590009 0.5076325 0.135973

L_Y_barU_Y_bar

1 3.225701 4.240966

$W_h

$f_h

$y_bar_h

$s_2_h

$W_h_y_bar_h

A.2. 例2的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

>W_h = c(0.082, 0.065, 0.137, 0.056, 0.118, 0.116, 0.17, 0.098, 0.088, 0.07)

>n_h = c(16, 13, 27, 11, 24, 23, 34, 20, 18, 14)

>y_bar_h = c(89, 56, 102, 76, 97, 79, 83, 52, 36, 52)

>s_h = c(105, 74, 186, 97, 106, 89, 112, 73, 44, 65)

> f = 0 # 各层抽样比可以忽略

> alpha = 0.05

>

> result = Compute_Y_bar_prop_from_y_bar_h_s_h_st(W_h, n_h, y_bar_h, s_h, f, alpha); result

t n y_bar_st v_y_bar_st se_y_bar_st L_Y_barU_Y_bar

1 1.959964 200 75.792 59.46035 7.711054 60.67861 90.90539

A.3. 例3的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

>W_h = c(0.7, 0.3)

>n_h = c(120, 80)

>y_bar_h = c(253.4, 329.4)

>s_h = c(231, 367)

> N = 8000

> alpha = 0.05

>

> result = Compute_Y_bar_srs_pst(W_h, n_h, y_bar_h, s_h, N, alpha); result

t n f s2 s y_bar_pst v_y_bar_pst

1 1.959964 200 0.025 86772.05 294.571 276.2 381.8343

se_y_bar_pst L_Y_bar_pst U_Y_bar_pst y_barv_y_bar

1 19.54058 237.9012 314.4988 283.8 423.0137

se_y_bar L_Y_bar_srs U_Y_bar_srs

120.5673 243.4888 324.1112

A.4. 例4的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

>n_h = c(10, 10, 10, 10)

>N_h = c(200, 400, 750, 1500)

> alpha = 0.05

> y1 = c(0, 0, 0, 1, 0, 0, 0, 1, 0, 0); y1

> y2 = c(0, 1, 0, 0, 0, 0, 0, 0, 1, 0); y2

> y3 = c(1, 1, 0, 0, 0, 0, 1, 0, 1, 0); y3

> y4 = c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0); y4

>y_matrix = rbind(y1, y2, y3, y4); y_matrix

y1 0 0 0 1 0 0 0 1 0 0

y2 0 1 0 0 0 0 0 0 1 0

y3 1 1 0 0 0 0 1 0 1 0

y4 1 0 0 0 0 0 0 0 0 0

>

>result_P = Compute_P_st(y_matrix, n_h, N_h, alpha); result_P

$df_0

t N p_st v_p_st se_p_st Delta

1 1.959964 28500.2 0.004998324 0.07069883 0.1385672

gamma L_P U_P

1 0.6928358 0.06143284 0.3385672

$a_h

y1 y2 y3 y4

2241

$W_h

$f_h

$p_h

y1y2y3y4

0.2 0.2 0.4 0.1

$q_h

y1y2y3y4

0.8 0.8 0.6 0.9

$v_p_h

y1 y2 y3 y4

0.016888889 0.017333333 0.026311111 0.009933333

此外，我们也可以利用总体均值的分层随机抽样的点估计和区间估计的R函数（程序）Compute_Y_bar_st()来进行计算。

>result_Y_bar = Compute_Y_bar_st(y_matrix, n_h, N_h, alpha); result_Y_bar

$df_0

t N y_bar_stv_y_bar_st se_y_bar_st Delta

1 1.959964 2850 0.2 0.0049983240.07069883 0.1385672

gamma L_Y_bar U_Y_bar

1 0.6928358 0.06143284 0.3385672

$W_h

$f_h

$y_bar_h

$s2_h

$W_h_y_bar_h

A.5. 例5的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

>W_h = c(0.18, 0.21, 0.14, 0.09, 0.16, 0.22)

>n_h = c(30, 30, 30, 30, 30, 30)

>a_h = c(27, 28, 27, 26, 28, 29)

> N = 1650000

>N_h = N * W_h; N_h

>f_h = n_h / N_h; f_h

>f_h = c(0, 0, 0, 0, 0, 0); f_h # 各层抽样比可以忽略

> alpha = 0.05

> result = Compute_P_from_a_h_st(W_h, n_h, a_h, f_h, alpha); result

$df_0

tp_st v_p_st se_p_st L_P U_P

1 1.959964 0.924 0.0003969808 0.01992438 0.8849489 0.9630511

$p_h

$q_h

$W_h_p_h

A.6. 例6的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

>W_h = c(0.7, 0.3)

>a_h = c(1, 2)

>n_h = c(43, 57)

>f_h = c(0, 0) # 各层抽样比可以忽略

> f = 0 # 抽样比可以忽略

> alpha = 0.05

>

> result = Compute_P_srs_pst(W_h, a_h, n_h, f_h, f, alpha); result

$df_srs

t a n p q v_p se_p L_P_srs

1 1.959964 3 100 0.03 0.97 0.0002939394 0.01714466 -0.003602918

U_P_srs

1 0.06360292

$df_pst

p_pst v_p_pst_1 se_p_pst_1 L_P_pst_1U_P_pst_1

1 0.02680539 0.0002692841 0.01640988 -0.005357385 0.05896816

v_p_pst_2 se_p_pst_2 L_P_pst_2U_P_pst_2

1 0.0003194205 0.01787234 -0.008223753 0.06183452

$p_h

$q_h

$s2_h

A.7. 例7的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

> n = 40

>N_h = c(200, 400, 750, 1500)

> y1 = c(10, 40, 40, 110, 15, 10, 40, 80, 90, 0); y1

> y2 = c(50, 130, 130, 80, 100, 55, 160, 85, 160, 170); y2

> y3 = c(180, 260, 260, 0, 140, 60, 200, 180, 300, 220); y3

> y4 = c(50, 35, 15, 0, 20, 30, 25, 10, 30, 25); y4

>y_matrix = rbind(y1, y2, y3, y4); y_matrix

y1 10 40 40110 15 10 40 80 90 0

y2 50130130 80100 55160 85160 170

y3180260260 0140 60200180300 220

y4 50 35 15 0 20 30 25 10 30 25

>

> result = Compute_nh_given_n_Y_bar_st(y_matrix, n, N_h); result

$df_0

L N y_bar_st sum_W_h_s_h

1 4 2850 78.77193 40.61926

$n_h_prop

n1 n2 n3 n4

2.8070185.614035 10.526316 21.052632

$n_h_Neyman

n1 n2 n3 n4

2.6165486.113663 24.0632337.206555

$W_h

$y_bar_h

y1 y2 y3 y4

43.5 112.0 180.024.0

$s2_h

y1 y2 y3 y4

1433.6111 1956.6667 8622.2222193.3333

$W_h_s_h

y1 y2 y3 y4

2.6570576.208312 24.4357697.318124

A.8. 例8的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

>N_h = c(200, 400, 750, 1500)

> alpha = 0.05

> gamma = 0.1

> y1 = c(10, 40, 40, 110, 15, 10, 40, 80, 90, 0); y1

> y2 = c(50, 130, 130, 80, 100, 55, 160, 85, 160, 170); y2

> y3 = c(180, 260, 260, 0, 140, 60, 200, 180, 300, 220); y3

> y4 = c(50, 35, 15, 0, 20, 30, 25, 10, 30, 25); y4

>y_matrix = rbind(y1, y2, y3, y4); y_matrix

y1 10 40 40110 15 10 40 80 90 0

y2 50130130 80100 55160 85160 170

y3180260260 0140 60200180300 220

y4 50 35 15 0 20 30 25 10 30 25

>

> result = Compute_n_nh_Y_bar_st(y_matrix, N_h, alpha, Given = gamma, input = gamma); result

$df_0

t L N y_bar_st sum_W_h_s2_h sum_W_h_s_h V

1 1.959964 4 2850 78.77193 2745.984 40.61926 16.15276

n0_prop n_prop n_Neyman

1 170.0009 160.4313 96.39512

$n_h_prop

n1 n2 n3 n4

11.29825 22.59649 42.36842 84.73684

$n_h_Neyman

n1 n2 n3 n4

6.34513 14.82563 58.35334 17.47590

$W_h

$y_bar_h

y1 y2 y3 y4

43.5 112.0 180.024.0

$s2_h

y1 y2 y3 y4

1433.6111 1956.6667 8622.2222193.3333

$W_h_s2_h

y1 y2 y3 y4

100.6043274.6199 2269.0058101.7544

$W_h_s_h

y1 y2 y3 y4

2.6570576.208312 24.4357697.318124

A.9. 例9的详细的R程序输入及输出的结果

>rm(list = ls(all = TRUE))

> source(subfunctions.R)

>n_h = c(10, 10, 10, 10)

>N_h = c(200, 400, 750, 1500)

>p_h = c(0.2, 0.2, 0.4, 0.1)

> alpha = 0.05

> Delta = 0.05

>

> result = Compute_n_nh_P_st(N_h, p_h, alpha, Given = Delta, input = Delta); result

$df_0

t L N p_st V sum_W_h_p_h_q_h

1 1.959964 4 28500.2 0.0006507944 0.1442105

sum_W_h_sqrt_p_h_q_h n0_prop n_prop n_Neyman

1 0.3710258 221.5915 205.6054 196.2663

$n_h_prop

n1 n2 n3 n4

14.4561428.9122854.21053 108.42105

$n_h_Neyman

n1 n2 n3 n4

14.90415 29.80830 68.45169 83.83586

$W_h

$q_h

$W_h_p_h

$W_h_p_h_q_h

$W_h_sqrt_p_h_q_h

NOTES

^*第一作者。

^#通讯作者。