【抽样调查】实验

小白在进击

已于 2023-06-16 11:13:00 修改

阅读量896

点赞数 3

文章标签：矩阵线性代数

于 2023-06-15 21:04:50 首次发布

本文链接：https://blog.csdn.net/m0_69442905/article/details/131235713

版权

文章目录

7、往届试卷
实验答案

1、数组矩阵简单抽样

（1）构造数组

c(1, 2, 3, 4, 5, 6)
seq(from=1, to=10, by= 2)
seq(1,100,4)
1:10

（2）构造矩阵

matrix(1:12, 3, 4)
# 三行四列

（3）产生来自正态分布的随机数

rnorm(100, 0, 1) # 产生100个0 - 1之间服从正态分布的随机数

m = rnorm(100, 0, 1)
n1 = sample(m, 10)  # 从m中抽取10个数	 
n1 				   # 打印

（4）从正态总体中抽取若干个样本

m = rnorm(100, 0, 1)
n1 = sample(m, 10)
n2 = sampel(m, 10)
mean1 = mean(n1)     # 计算n1抽样的平均数
mean2 = mean(n2)	#  计算n2抽样的平均数
mean3 = mean(c(n1, n2))  # 计算n1和n2构成数组的元素平均数


# 对总体均值于样本均值的平均数均值进行比较
> mean3
[1] 0.52052
mean0 = mean(m) # 总体均值
mean0
[1] 0.05755668
# 总体均值 0.05755668 < 样本均值 0.52052

（5）对矩阵的行或列进行统计计算

2、R软件作图

require() # 要求 用于加载库
graphics  # 绘图 库名称
expression # 表示 显示列名称

（1）用curve函数画Gamma(1 / 2, 1)… Gamma(13 / 2, 1)

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-iqYJOXau-1686833845961)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230613162257583.png)]$

require(graphics) # 加载库
par(lwd=2)        # 打开绘图板
color=hsv(seq(0, 1, 1/ 7))   # 设置颜色
curve(dgamma(x, 1 / 2, 1), xlim = c(0, 10), ylim = c(0, 0.7),  # 坐标刻度设定
	ylab = expression(paste("Ga(",alpha,",",lambda,")")),    # 坐标名称设定
	col = color[1])
number<-7
for (i in 2 : 7){
	curve(dgamma(x, (i - 1) + 1 / 2, 1), col = color[i], add = T)
	legend("topright", c("Ga(1 / 2 , 1)", "Ga(3 / 2, 1)", "Ga(5 / 2, 1)", "Ga(7 / 2, 1)", "Ga(9 / 2, 1)", "Ga(11 / 2, 1)", "Ga(13 / 2, 1)"),
	col = color, bty = "n", cex = 1.2, pch = 15)
}

（2）用curve函数画Be(1 / 2, 1 / 2), Be(3 / 2.，3 / 2)

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-QInDbmMt-1686833746631)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230613163926996.png)]$

require(graphics) # 加载库
par(lwd = 2) 	# 打开绘图板
color = hsv(seq(0, 1, 1 / 7))  # 设置颜色

Be = function(x, a, b){
	dbeta(x, a, b)
}


curve(Be(x, 1 / 2, 1 / 2), xlim = c(0, 10), ylim = c(0, 2), # 画线+设置坐标刻度
	ylab = expression(paste("Ga(",alpha,",",lambda,")")),	# 设置y轴名称
	col = color[1],						# 设置颜色
	main = "Beta-Density", lty = 1)  			# 设置图名称

curve(Be(x, 3 / 2, 3 / 2), col = color[2], add = T)   # add 图坐标以创建，添加线
legend("topright", c("Be(1 / 2, 1 / 2)", "Be (3 / 2,3 / 2)"), 
	col = color, bty = "n", cex = 1.2, pch = 15) # 设置左上角图例

（3）用plot函数画N(0, 1), t(20), t(10), t(5), t(2)

3、简单抽样（一）

replicate  # 重复

（1）从R中装载数据集"cars"

data("cars")
Y = cars[, 2] # 将数据集cars中的第二列数据构成Y数组

（2）iteration = 100, N = 50, n = 10 时的误差

N = length(Y)
n = 10
iteration=100
y = replicate(Iteration, sample(Y, n))
ybar = apply(y, 2, mean)   # 选择列来进行平局n
Eybar = mean(ybar)		   #计算期望
Dx = Eybar - mean(Y)
Dx

（3）iteration = 500, N = 50, n = 10 时的误差

data("cars")
Y = cars[, 2] 
N = length(Y)
n = 10
iteration = 500
y = replicate(iteration, sample(Y, n))
ybar = apply(y, 2, mean)
Eybar = mean(ybar)
Dx = Eybar - mean(Y)
Dx
[1] -0.1918
Eybar
[1] 42.7882
ybar

（4）iteration = 1000, N = 50, n = 10 时的误差

data("cars")
Y = cars[, 2]
N = length(Y)
n = 10
iteration = 1000
y = replicate(iteration, sample(Y, n))
ybar = apply(y, 2, mean)
Eybar = mean(ybar)
Dx = Eybar - mean(Y)
Dx
0.0226

（5）编写循环语句，计算iteration从100到2000，步长为100的误差

data("cars")
Y = cars[, 2]
N = length(Y)
n = 10

for (iteration in seq(100, 2000, 100)){
	y = replicate(iteration, sample(Y, n))
	ybar = apply(y, 2, mean)
	Eybar = mean(ybar)
	Dx = Eybar - mean(Y)
	print(Dx)
}

（6）编写程序，画出(5)中误差变化的情况

require(graphics)
data("cars")
Y = cars[, 2]
N = length(Y)
n = 10

for (iteration in seq(100, 2000, 100)){
	y = replicate(iteration, sample(Y, n))
	ybar = apply(y, 2, mean)
	Eybar = mean(ybar)
	Dx[iteration / 100] = Eybar - mean(Y)
	
}
plot(x = seq(100, 2000, 100), y = Dx, type = 'o', xlab = 'iteration', ylab = 'Dx', col = 'red')

（7）将（5）中的结果记录到"record"文件

write.table(Dx, file = 'record.txt', col.names=F)

4、简单抽样（二）

1、产生200个均值15，标准差1的正态随机数

mu = 15 	# 均值
sigma = 1   # 标准差
N = 200
r = rnorm(N, mu , sigma)
r

2、用简单随机抽样方法（无放回），抽取样本容量为20的样本

mu = 15 	# 均值
sigma = 1   # 标准差
N = 200
n = 20
Y = rnorm(N, mu, sigma)
r = sample(Y, n)
r

3、抽培养如2所示样本100个，分别用for循环、replicate（）实现

# for 循环
mu = 15
sigma = 1
N = 200
n = 20
Y = rnorm(N, mu, sigma)
number = 100
r = replicate(number, sample(Y, n))
r

# replicate()
mu = 15
sigma = 1
N = 200
n = 20
Y = rnorm(N, mu, sigma)
number = 100
for (i in 1 : number) r[, i] = sample(Y, n)
r

4、计算100个样本中每个样本的样本均值、样本标准差

mu = 15
sigma = 1
N = 200
n = 20
Y = rnorm(N, mu, sigma)
number = 100
y = replicate(number, sample(Y, n))
ybar = apply(y, 2, mean)       # 计算每一个样本均值
s = apply(y, 2, sd)			  # 计算每一个样本标准差
ybar

5、根据每个样本，计算总体均值的置信水平95%的置信区间

mu = 15
sigma = 1
N = 200
n = 20
Y = rnorm(N, mu, sigma)
number = 100
y = replicate(number, sample(Y, n))
ybar = apply(y, 2, mean)
s = apply(y, 2, sd)

f = n / N
ybar = apply(y, 2, mean)
s = apply(r, 2, sd)
yl = ybar - sqrt((1 - f) / n) * s
yu = ybar + sqrt((1 - f) / n) * s
yl		# 置信下区间
yu		# 置信上区间

6、在平面直角坐标系中画出100个置信区间

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CHjoIt1c-1686833915840)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230614215217276.png)]$

mu = 15
sigma = 1
N = 200
n = 20
Y = rnorm(N, mu, sigma)
number = 100
y = replicate(number, sample(Y, n))
s = apply(y, 2, sd)

f = n / N
ymean = apply(y, 2, mean)
s = apply(r, 2, sd)
yl = ymean - sqrt((1 - f) / n) * s
yu = ymean + sqrt((1 - f) / n) * s

f = n / N
u = 1.96

# 制作图表
plot (1, xlim=c(0.5, number + 0.5), ylim = c(min(yl) - 0.5, max(yu) + 0.5),
	type = "n", xlab = "Trials", ylab = "Confidence Intervals",
	main = expression(paste("Confidence Interval of 1-", alpha, sep=" ")))

# 绘制线条
for (i in 1 : number){
	arrows(i, yl[i], i, yu[i], length = 0.1, angle = 90, code = 3,
	col = ifelse(ymean > yl[i] & ymean < yu[i], "blue", "red"))
	points(i, ymean[i])
}
#Sys.sleep(0.5)

7、计算100个置信区间的置信概率

mu = 15
sigma = 1
N = 200
n = 20
Y = rnorm(N, mu, sigma)
number = 100
y = replicate(number, sample(Y, n))
s = apply(y, 2, sd)

f = n / N
ymean = apply(y, 2, mean)
s = apply(r, 2, sd)
yl = ymean - sqrt((1 - f) / n) * s
yu = ymean + sqrt((1 - f) / n) * s

f = n / N
u = 1.96
plot (1, xlim=c(0.5, number + 0.5), ylim = c(min(yl) - 0.5, max(yu) + 0.5),
	type = "n", xlab = "Trials", main = expression(paste("Confidence Interval of 1-", alpha, sep=" ")),
	ylab = "Confidence Intervals")

cn = 0 # 1

for (i in 1 : number){
	arrows(i, yl[i], i, yu[i], length = 0.1, angle = 90, code = 3,
	col = ifelse(ymean > yl[i] & ymean < yu[i], "blue", "red"))
	points(i, ymean[i])
	cn = cn + as.numeric(ymean > yl[i] & ymean < yu[i]) # 1
}

abline(h = ymean, lty = 2,) # 1
cp = cn / number # 1 
cat ("Confidence Probability = ", cp) # 1

5、中心极限定理模拟

（1）产生来自常用分布的随机数

N = 1000
# 从 1-6 中随机抽取一个，重复 1000 次
Y1 = replicate(N, sample(1:6, 1)) # 次数, 样本
Y1

（2）画直方图，并设置图形参数

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-v2CRXamF-1686834005316)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230614223136724.png)]$

N = 1000
# 从 1-6 中随机抽取一个，重复 1000 次
Y1 = replicate(N, sample(1:6, 1)) # 次数, 样本

hist(Y1, probability = T, breaks = 0.5 + (0:6) * 1,
	xlab = expression(paste(Y[1], "=", X[1])),
	ylab = "直方图", ylim = c(0, 0.25),
	col = "light blue")

（3）画函数分布的核密度曲线，并设置图形参数

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-jb5Rn37z-1686834005319)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230614224003299.png)]$

### 核密度曲线
# 密度函数density
lines(density(Y1), lwd = 3, col = "blue")

（4）画正态分布的密度函数，并设置图形参数

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-QI92sZ9g-1686834005320)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230614224022561.png)]$

### 正态分布密度函数， 样式，粗细，颜色
# 正态分布密度函数curve
curve(dnorm(x, mean(Y1), sd(Y1)), lty = 3, lwd = 3, col = "red", add = T)

6、分层抽样

问题一

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-47q17ok8-1686834065957)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230614224750954.png)]$

# 1.1
# Correlation, Variance and Covariance (Matrices)
# 相关、方差和协方差（矩阵）
data = c(43, 84, 98, 0, 10, 44, 0, 124, 13, 0, 
	   50, 147, 62, 87, 84, 158, 170, 104, 56, 160,
	   228, 262, 110, 232, 139, 178, 334, 0, 63, 220, 
	   17, 34, 25, 34, 36, 0, 25, 7, 15, 31)
sam = matrix(data, nrow = 4, ncol = 10, byrow = TRUE) # byrow 指定数组横竖成矩阵

total = c(1411, 4700, 2558, 14997)
Wh = c()
Sh = c()
ybar = c()
nh = c(10, 10, 10, 10)
N = sum(total)

for (i in (1 : length(total))){
	Wh[i] = total[i] / N      # 每一层占总体的比例
	Sh[i] = var(sam[i, ])     # 计算每一层的方差
	ybar[i] = mean(sam[i, ])  # sam[i,] 输出第i行的数据,每一层的平均值
}

Y = sum(total * ybar) 		  # 牛总数
Vss = sum(Wh ^ 2 * Sh / nh) - sum(Wh * Sh / N)
sd = sqrt(Vss)
sd 
Mean = sum(Wh * ybar)
Mean
cv = sd / Mean # 变异系数
cv

# 1.2
s2 = 1 / (sum(nh) - 1) * sum ((sam - mean(sam)) ^ 2)
s2 
Vsrs = (1 - sum(nh) / N) * (s2 / sum(nh))
Vsrs
deff = Vss / Vsrs
deff

# 1.3
Vyst = 1 / sum(nh) * (sum(Wh * sqrt(Sh))) ^ 2 - 1 / N * (sum (Wh * Sh))
change = Vyst / Vss
change

问题二

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GTAWA4M4-1686834065957)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230614224759788.png)]$

# 1.1
pre = c(550, 720, 1500, 1020, 620, 980, 928, 1200, 1350, 1750, 670, 729, 1530)
now = c(610, 780, 1600, 1030, 600, 1050, 977, 1440, 1570, 2210, 980, 865, 1710)
y_tilde = sum(now) / sum(pre) * 128200
y_tilde
t.test(now / pre * 128200)

data:  now/pre * 128200
t = 32.17, df = 12, p-value = 5.142e-13
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:  # 95%置信区间
 135639 155347
sample estimates:			 # 样本估计
mean of x 					
   145493

# 1.2
pre = c(550, 720, 1500, 1020, 620, 980, 928, 1200, 1350, 1750, 670, 729, 1530)
now = c(610, 780, 1600, 1030, 600, 1050, 977, 1440, 1570, 2210, 980, 865, 1710)
y_tilde = sum(now) / sum(pre) * 128200
y_tilde
t.test(now / pre * 128200)

xbar = mean(pre)
ybar = mean(now)
N = 123
n = 13
X = 128200
sxy = cov(pre, now)
varx = var(pre)
vary = var(now)
R = ybar / xbar
Yrbar = R * X
a = sxy / varx
Xbar = X / N
ylr = N ^ 2 * (1 - n / N) 
se = sqrt((n - 1) * (vary - a * sxy) / (n - 2))
vYlr = N ^ 2 * ( 1 - n / N) * se ** 2 / n
sd = sqrt(vYlr)
alpha = 0.05
(C = Q(Yrbar, sd, alpha))

问题三

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AfpLaI78-1686834065958)(C:\Users\Lenovo\AppData\Roaming\Typora\typora-user-images\image-20230614224812351.png)]$

N1=120
 N2=180
 n1=6
 n2=9
 y1=c(210,160,75,280,300,190)
 x1=c(204,143,82,256,275,198)
 y2=c(150,200,125,60,110,100,180,75,90)
 x2=c(137,189,119,63,103,107,159,63,87)
 X1=24500
 X2=21200
 f1=n1/N1
 f2=n2/N2
 y1bar<-mean(y1) 
 x1bar<-mean(x1) 
 y2bar<-mean(y2) 
 x2bar<-mean(x2)
 y1var=var(y1)
 x1var=var(x1)
 y2var=var(y2)
 x2var=var(x2)
 sxy1=cov(x1,y1)
 sxy2=cov(x2,y2)
 # (1)分别比估计
 R1=y1bar/x1bar
 R2=y2bar/x2bar
 yrs=X1*R1+X2*R2
 yrs#总体估计量
[1] 48206.45
 yrsvar=(N1^2*(1-f1)/n1*(y1var+R1^2*x1var-2*R1*sxy1)
 +N2^2*(1-f2)/n2*(y2var+R2^2*x2var-2*R2*sxy2))
 yrssd=sqrt(yrsvar)
 yrssd#标准误
[1] 751.4072
 
 # （2）联合比估计
 N=N1+N2
 W1=N1/N
 W2=N2/N
 xstbar=W1*x1bar+W2*x2bar
 ystbar=W1*y1bar+W2*y2bar
 RC=ystbar/xstbar
 yrc=X*RC
 yrc#总体估计量
[1] 135240.7
 yrcvar=(N1^2*(1-f1)/n1*(y1var+RC^2*x1var-2*RC*sxy1)
 +N2^2*(1-f2)/n2*(y2var+RC^2*x2var-2*RC*sxy2))
 yrcsd=sqrt(yrcvar)
 yrcsd#标准误
[1] 746.388

7、往届试卷

1、R软件基本操作

读取数据

getwd();
data = read.csv("file_name")

（2）用三种方法构造向量（1， 1，3，3， 3， 7， 7， 7）

# 方法1
a = c(1, 1, 3, 3, 3, 7, 7, 7)
# 方法2
b = c(1, 1)
c = c(3, 3, 3)
d = c(7, 7, 7)
e = c(b, c, d)
# 方法3
f = c(1, 1)
g = c(3, 7)
h = c(f ,rep(g, each=3))

将两个向量合并成一个矩阵

a = c(4, 1, 7)
b = c(5, 3, 2)
d = c(a, b)
my_matrix = matrix(d, nrow = 2, ncol = 3, byrow = TRUE)
my_matrix
     [,1] [,2] [,3]
[1,]    4    1    7
[2,]    5    3    2

用1-12生成4*3的矩阵

my_matrix = matrix(1:12, nrow = 4, ncol = 3, byrow = TRUE)

产生来自正态分布N(5, 2)的100个随机数

sigma = sqrt(2)
N = 100
r = rnorm(N, 5, sigma)
r

2、简单抽样

# (1)
Y = c(2, 3, 4, 4, 5, 5, 5, 5, 7, 7, 7, 8, 9, 11, 12)
N = length(Y)
s = var(Y)
s
[1] 8.066667


# (2)
n = 3 
y = sample(Y, n)
f = n / N
vy = ((1 - f) / n) * s
vy

3、分层抽样

N1=120
 N2=180
 n1=6
 n2=9
 y1=c(210,160,75,280,300,190)
 x1=c(204,143,82,256,275,198)
 y2=c(150,200,125,60,110,100,180,75,90)
 x2=c(137,189,119,63,103,107,159,63,87)
 X1=24500
 X2=21200
 f1=n1/N1
 f2=n2/N2
 y1bar<-mean(y1) 
 x1bar<-mean(x1) 
 y2bar<-mean(y2) 
 x2bar<-mean(x2)
 y1var=var(y1)
 x1var=var(x1)
 y2var=var(y2)
 x2var=var(x2)
 sxy1=cov(x1,y1)
 sxy2=cov(x2,y2)
 # (1)分别比估计
 R1=y1bar/x1bar
 R2=y2bar/x2bar
 yrs=X1*R1+X2*R2
 yrs#总体估计量
[1] 48206.45
 yrsvar=(N1^2*(1-f1)/n1*(y1var+R1^2*x1var-2*R1*sxy1)
 +N2^2*(1-f2)/n2*(y2var+R2^2*x2var-2*R2*sxy2))
 yrssd=sqrt(yrsvar)
 yrssd#标准误
[1] 751.4072
 
 # （2）联合比估计
 N=N1+N2
 W1=N1/N
 W2=N2/N
 xstbar=W1*x1bar+W2*x2bar
 ystbar=W1*y1bar+W2*y2bar
 RC=ystbar/xstbar
 yrc=X*RC
 yrc#总体估计量
[1] 135240.7
 yrcvar=(N1^2*(1-f1)/n1*(y1var+RC^2*x1var-2*RC*sxy1)
 +N2^2*(1-f2)/n2*(y2var+RC^2*x2var-2*RC*sxy2))
 yrcsd=sqrt(yrcvar)
 yrcsd#标准误
[1] 746.388

4、不等概率抽样

m = c(12, 11, 20, 9, 17, 22, 10, 10, 8, 32)
N = 10
size = c()
while (length(size) < 4){
	i = sample(1 : N, 1)
	Mi = sample(1, m[i], 1)
	mi = sample(1 : 32, 1)
	if (Mi >= mi) {
		size = append(size, i)
	}
}
size

实验答案


getwd() # 将.csv文件复制到获取到的路径下
data = read.csv("car.test.csv")
price = data[,2]
weight = data[,7]
price                # 观察变量price
weight		   # 观察变量weight

mean_price = mean(price)          # 计算price 平均值
mean_wiegith = mean(weight) 	   # 计算weight 平均值
mean_price 
mean_weight
	
sd_price = sd(price)        # 计算price 标准差
sd_weight = sd(weight)      # 计算weight 标准差

sd_price
sd_weight

N = length(price)
n = 10
y = sample(price, n)	    # 抽取10个简单随机样本
sub_mean_price = mean(y)    # 样本价格的简单估计
sub_mean_price 

mean_price - sub_mean_price # 总体均值 - 样本均值 = 误差


# 方法1
a = c(1, 1, 2, 2, 2, 7, 7)
a
# 方法2
b = c(1, 1)
c = c(2, 2, 2)
d = c(7, 7)
e = c(b, c, d)
e
# 方法3
f = c(1, 1, 2)
g = c(2, 7)
h = c(f ,rep(g, each=2))
h


a = c(2, 1, 5)
b = c(1, 3, 4)
d = c(a, b)
my_matrix = matrix(d, nrow = 2, ncol = 3, byrow = TRUE)
my_matrix


my_matrix = matrix(1:12, nrow = 4, ncol = 3, byrow = TRUE)
my_matrix

sigma = 1
N = 100
r = rnorm(N, 5, sigma)
r

y=rnorm(100,5,1)	#rnorm(100,mean=5,sd=1)
# 计算其概率密度在 x = 3 处的函数值；
dnorm(3,y)

# 作图
x<-ecdf(y)
plot(x,col="dark blue",lwd=3)
lines(y,pnorm(y,5,1),col="red")

N1=120
N2=180
n1=6
n2=9
y1=c(210,160,75,280,300,190)
x1=c(204,143,82,256,275,198)
y2=c(150,200,125,60,110,100,180,75,90)
x2=c(137,189,119,63,103,107,159,63,87)
X1=24500
X2=21200
f1=n1/N1
f2=n2/N2
y1bar<-mean(y1) 
x1bar<-mean(x1) 
y2bar<-mean(y2) 
x2bar<-mean(x2)
y1var=var(y1)
x1var=var(x1)
y2var=var(y2)
x2var=var(x2)
sxy1=cov(x1,y1)
sxy2=cov(x2,y2)
# (1)分别比估计
R1=y1bar/x1bar
R2=y2bar/x2bar
yrs=X1*R1+X2*R2
yrs#总体估计量
yrsvar=(N1^2*(1-f1)/n1*(y1var+R1^2*x1var-2*R1*sxy1)
+N2^2*(1-f2)/n2*(y2var+R2^2*x2var-2*R2*sxy2))
yrssd=sqrt(yrsvar)
yrssd#标准误

# （2）联合比估计
N=N1+N2
X=X1+X2
W1=N1/N
W2=N2/N
xstbar=W1*x1bar+W2*x2bar
ystbar=W1*y1bar+W2*y2bar
RC=ystbar/xstbar
yrc=X*RC
yrc#总体估计量
yrcvar=(N1^2*(1-f1)/n1*(y1var+RC^2*x1var-2*RC*sxy1)
+N2^2*(1-f2)/n2*(y2var+RC^2*x2var-2*RC*sxy2))
yrcsd=sqrt(yrcvar)
yrcsd#标准误

m = c(13, 10, 21, 9, 15, 23, 9, 9, 7, 33)
N = 10
n = 4
size = c()
while (length(size) < 4){
	i = sample(1 : N, 1)
	Mi = sample(1 : m[i], 1)
	mi = sample(1 : 32, 1)
	if (Mi >= mi) {
		size = append(size, i)
	}
}
size