第5章：参数估计与假设检验

### 数据准备 ###
# 清空工作空间
rm(list = ls())

案例引入

本章主要包括三个简短的案例，以下对每个案例及数据进行简要介绍。

案例一：手机续航时间研究

手机续航能力是指手机在正常工作时的待机时间。它与手机本身的耗电量和所配电池容量的大小有关。在电池容量相同的情况下，耗电量越小，则续航能力越强；在耗电量相同的情况下，电池容量越大，待机时间越长。近年来，各个品牌厂家都致力于打造续航能力更强的新款手机，从而满足现代人们对手机设备的长时间使用需求，提高自身的市场竞争力。

battery <-c(20.2,19.4,18.6,18.2,19.8,21.4,18.2,19.8,17.4,19.4,21.4,19.0,19.0,17.0,23.0,22.6,21.4,20.2,22.6,20.2,21.4,20.2,18.2,20.2,19.0,19.8,19.8,20.6,20.6,18.6,19.0,21.4,20.6,21.8,19.4,19.8,21.8,19.4,17.4,20.2,19.8,17.8,17.4,19.8,20.2,19.0,20.2,18.2,21.0,19.8)
# 显示续航时间样本数据
print(battery)

##  [1] 20.2 19.4 18.6 18.2 19.8 21.4 18.2 19.8 17.4 19.4 21.4 19.0 19.0 17.0 23.0 22.6 21.4 20.2 22.6 20.2 21.4 20.2
## [23] 18.2 20.2 19.0 19.8 19.8 20.6 20.6 18.6 19.0 21.4 20.6 21.8 19.4 19.8 21.8 19.4 17.4 20.2 19.8 17.8 17.4 19.8
## [45] 20.2 19.0 20.2 18.2 21.0 19.8

本案例研究的问题是判断这一批手机的续航时间是否达到了声称的20小时。

案例二：减肥药疗效研究

随着生活水平的提高，油腻、高热量的食物成为了很多年轻人的便捷之选，肥胖患者日剧增多。然而由此引发的肥胖问题会极大危害到身体健康，仅在2015年，全球因肥胖造成死亡的人数就已超过400万，世界卫生组织已将超重、肥胖定义为一种慢性病。因此，减肥不仅是一个时尚与外貌相关的话题，更关乎身体的健康。为了应对肥胖带来的身体健康问题，制药公司纷纷推出减肥特效药。减肥药药效如何，需要大量试验验证，得出结论。

在减肥的案例中，将会使用以下减肥药药效的数据。一个合规的减肥药要在美国市场上合法上市，必须有美国食品药品管理局（FDA）的批准。为此，FDA对使用该药物的受试者做了一组试验，记录了其在使用减肥药前后的体重（单位：kg）。试验的数据如下所示：

受试者编号	1	2	3	4	5	6	7	8	9	10
服用药物前的体重x	43	55	49	62	59	49	57	54	55	48
服用药物后的体重y	50	59	55	60	58	54	56	53	61	51
差 d = x-y	-7	-4	-6	2	1	-5	1	1	-6	- 3

受试者在服用减肥药后是否体重发生了显著下降，这是本章所要探究的问题。

案例三：红楼梦作者争议

《红楼梦》是一部章回体长篇小说，被列为中国古代四大名著之首。《红楼梦》文学价值极高，围绕作品本身的争议也很多。其中一大争议便是作者归属问题。

本案例收集了《红楼梦》中有代表性的若干个文言虚词，并希望探究这些虚词在《红楼梦》前后的使用差异。在本章中，我们选择具有代表性的“之”和“亦”两个虚词的章回词频，通过对比词频在整本书前后出现的差异，可以了解作者用语习惯的改变。具体而言，统计每个虚词在每个章回的总词频，将全书1-40回、41-80回、81-120回切分为三个总体。全书1-40回、41-80回、81-120回中两个字的平均字频如下图所示：

众多大家各执一词，学术界仍无定论。本章从文本分析的角度，分析作者用语习惯的改变，探索《红楼梦》作者归属的问题。

5.2 参数估计

5.2.1 矩估计

电池续航时间

2030年，苹果公司发布新一代手机iPhone SSR。该手机声称续航能力比上一版本（20小时）有较大提升。为调查市面上售卖的苹果手机是否满足这样的续航能力，我们收集并测试了50部苹果手机的续航时间，得到如下数据。假设苹果手机的续航时间服从指数分布，那么使用矩估计推断该苹果手机的平均续航时间是多久？

# 例1
battery <- c(20.2,19.4,18.6,18.2,19.8,21.4,18.2,19.8,17.4,19.4,21.4,19,19,17,
             23,22.6,21.4,20.2,22.6,20.2,21.4,20.2,18.2,20.2,19,19.8,19.8,20.6,
             20.6,18.6,19,21.4,20.6,21.8,19.4,19.8,21.8,19.4,17.4,20.2,19.8,17.8,
             17.4,19.8,20.2,19,20.2,18.2,21,19.8)
mean(battery)  # 样本均值

## [1] 19.824

答：可以计算得到以上50部手机的续航时间均值为19.824。

正态总体参数估计

使用模拟的方法，估计来自于正态总体的均值\(\mu\)和方差\(\sigma^2\)。

# 例2
set.seed(123)
data1 <- rnorm(n = 1000, mean = 5, sd = 2)  # 产生服从正态分布的随机数
mean(data1)  # 样本均值

## [1] 5.032256

var(data1)  # 样本方差

## [1] 3.933836

产生1000个服从\(N(5,4)\)的随机数，分别求这个样本的一阶原点矩（均值）和二阶中心矩（方差），作为总体均值和方差的估计。结果显示，对总体均值的估计为5.03，总体方差的估计为3.93。

5.2.2 极大似然估计

苹果手机续航时间

假设苹果公司发布的新一代手机iPhone SSR续航时间是服从正态分布\(N(\mu, \sigma^2)\)，其中，\(\mu\) 是总体均值，\(\sigma^2\) 是总体方差。请根据案例数据中的50个续航时间样本使用最大似然估计的方法，估计总体的均值和方差。

cat('苹果手机电池续航数据为: ', battery, '\n')

## 苹果手机电池续航数据为:  20.2 19.4 18.6 18.2 19.8 21.4 18.2 19.8 17.4 19.4 21.4 19 19 17 23 22.6 21.4 20.2 22.6 20.2 21.4 20.2 18.2 20.2 19 19.8 19.8 20.6 20.6 18.6 19 21.4 20.6 21.8 19.4 19.8 21.8 19.4 17.4 20.2 19.8 17.8 17.4 19.8 20.2 19 20.2 18.2 21 19.8

# 续航时间均值的极大似然估计是
mu_mle <- mean(battery)
# 续航时间方差的极大似然估计是
sigma_mle <- mean((battery - mean(battery)) ** 2)
cat('均值的最大似然估计值是: ', mu_mle, 
    '方差的最大似然估计值是: ', sigma_mle, '\n')

## 均值的最大似然估计值是:  19.824 方差的最大似然估计值是:  1.948224

根据极大似然估计的方法，可以得到手机电池续航时间均值的最大似然估计值是19.824 方差的最大似然估计值是1.948224。

均匀分布总体参数估计

设\(x_1, ..., x_n\)是从均匀分布总体 \(U(0,\theta)\) 中抽出的一个样本，观测值如下：

6.0627 9.3764 2.6435 3.8009 8.0748 9.7808 9.5793 7.6273 5.0965 0.6448

6.4357 9.1591 0.9523 2.9537 7.6993 2.5589 5.179 6.7785 1.4723 7.0053

试求解\(\theta\)的最大似然估计。

set.seed(6)
theta = 10
# 生成的样本是
samples = runif(n=20, 0, theta)
samples = round(samples, 4)
cat('样本是: ', samples, '\n')

## 样本是:  6.0627 9.3764 2.6435 3.8009 8.0748 9.7808 9.5793 7.6273 5.0965 0.6448 6.4357 9.1591 0.9523 2.9537 7.6993 2.5589 5.179 6.7785 1.4723 7.0053

# theta 的最大似然估计值是
theta_mle = max(samples)
cat('theta的最大似然估计值是: ', theta_mle, '\n')

## theta的最大似然估计值是:  9.7808

在此案例数据中，均匀分布总体的参数\(\theta\)极大似然估计值是9.7808。

5.2.3 区间估计

手机续航时间的区间估计

已知上述案例中苹果手机的续航时间服从正态分布，其标准差为2小时，试求该型号手机续航时间的0.95置信区间。可以直接使用z.test()函数，也可以根据公式进行计算。

# 例5：续航时间区间估计
# 方法一：按照公式计算
mu <- mean(battery)
sigma <- 2
n <- 50
spread <- qnorm(1-0.05/2, 0, 1) * sigma / sqrt(n)
cat(sprintf('续航时间的置信区间是: [%s, %s]', 
            round(mu - spread, 4),
            round(mu + spread, 4), '\n'))

## 续航时间的置信区间是: [19.2696, 20.3784]

# 方法二：使用z.test()计算
library(BSDA)
interval_1 <- z.test(battery, sigma.x = sigma, conf.level = 0.95)$conf.int
cat(sprintf('续航时间的置信区间是: [%s, %s]', 
            round(interval_1[1], 4),
            round(interval_1[2], 4), '\n'))

## 续航时间的置信区间是: [19.2696, 20.3784]

无论按照公式计算，还是直接使用R中自带的函数，续航时间的置信区间估计都是 [19.2696, 20.3784]。

红楼梦词频区间估计

为了了解红楼梦在前、中、后三个部分对虚词使用的习惯，拟对不同章节的“之”字词频进行分析。假设小说对文言虚词的使用频率服从正态分布，试求第41-80回和第81-120回虚词“之”词频差的0.95置信区间。（注意：请区分两总体方差是否相等两种情况讨论）

# 例6
wenyan <- read.csv("./data/红楼梦虚词词频统计.csv", stringsAsFactors = F, fileEncoding = "utf-8")
freqx <- wenyan$`之`[41:80]
freqy <- wenyan$`之`[81:120]
m = length(freqx); n = length(freqy)
u_x = mean(freqx); u_y = mean(freqy)
s_x = sd(freqx); s_y = sd(freqy)

如果两部分数据的方差相等:

## 方法一：使用公式计算
s_w = sqrt(((m-1)*s_x^2 + (n-1)*s_y^2) / (m+n-2))
spread = qt(1 - 0.05/2, m+n-2) * s_w * sqrt(1/m + 1/n)
cat(sprintf('如果两部分数据方差相等，第41-80回和第81-120回“之”词频差的置信区间是: [%s, %s]', 
            round(u_x - u_y - spread, 4),
            round(u_x - u_y + spread, 4), '\n'))

## 如果两部分数据方差相等，第41-80回和第81-120回“之”词频差的置信区间是: [8.9824, 20.3176]

## 方法二：使用t.test()计算
interval_3 <- t.test(freqx, freqy, var.equal = T, conf.level = 0.95)$conf.int
cat(sprintf('如果两部分数据方差相等，第41-80回和第81-120回“之”词频差的置信区间是: [%s, %s]', 
            round(interval_3[1], 4),
            round(interval_3[2], 4), '\n'))

## 如果两部分数据方差相等，第41-80回和第81-120回“之”词频差的置信区间是: [8.9824, 20.3176]

如果两部分数据方差不等：

## 方法一：使用公式计算
s_o = sqrt(s_x^2 / m + s_y^2 / n)
freedom = s_o^4 / (s_x^4/(m^2 * (m-1)) + s_y^4/(n^2 * (n-1)))
freedom = round(freedom, 0)
spread = qt(1 - 0.05/2, freedom) * s_o
cat(sprintf('如果两部分数据方差不等，第41-80回和第81-120回“之”词频差的置信区间是: [%s, %s]', 
            round(u_x - u_y - spread, 4),
            round(u_x - u_y + spread, 4), '\n'))

## 如果两部分数据方差不等，第41-80回和第81-120回“之”词频差的置信区间是: [8.9291, 20.3709]

## 方法二：使用t.test()计算
interval_4 <- t.test(freqx, freqy, var.equal = F, conf.level = 0.95)$conf.int
cat(sprintf('如果两部分数据方差不等，第41-80回和第81-120回“之”词频差的置信区间是: [%s, %s]', 
            round(interval_4[1], 4),
            round(interval_4[2], 4), '\n'))

## 如果两部分数据方差不等，第41-80回和第81-120回“之”词频差的置信区间是: [8.928, 20.372]

如果两部分数据方差相等，第41-80回和第81-120回“之”词频差的置信区间是: [8.9824, 20.3176]；如果两部分数据方差不等，第41-80回和第81-120回“之”词频差的置信区间是: [8.9291, 20.3709]。在两总体方差不等的情况中，对t统计量的自由度做了近似处理，因此两种方法计算得到的结果稍有不同。

5.3 假设检验

苹果手机续航时间

假设苹果手机SSR的续航时间服从正态分布\(N(\theta, 0.8^2)\)，其中\(\theta\)的设计值不低于20小时。根据案例引入中给出的50次抽样检测数据，判断该苹果手机的续航时间是否满足出厂设计的要求？除了使用公式计算p值的方法。对于这个问题，你还可以尝试使用BSDA包中的z.test()和t.test()函数（假设总体方差未知的情况）进行检验。

battery <- c(20.2,19.4,18.6,18.2,19.8,21.4,18.2,19.8,17.4,19.4,21.4,19,19,17,
             23,22.6,21.4,20.2,22.6,20.2,21.4,20.2,18.2,20.2,19,19.8,19.8,20.6,
             20.6,18.6,19,21.4,20.6,21.8,19.4,19.8,21.8,19.4,17.4,20.2,19.8,17.8,
             17.4,19.8,20.2,19,20.2,18.2,21,19.8)
a <- 20
s <- 0.8
n <- 50
xbar <- mean(battery)
p_value <- pnorm(xbar, mean = a, sd = s/sqrt(n))
print(paste0("检验的p值为", round(p_value, 4)))

## [1] "检验的p值为0.0599"

检验的p值为0.0599。 当\(\alpha<0.0599\)时，\(z_\alpha<-1.56\)，由于拒绝域为\(W=\{z≤z_\alpha\}\)，于是观测值\(z=-1.56\)不在拒绝域里，应接受原假设； 当\(\alpha≥0.0599\)时，\(z_\alpha≥-1.56\)，由于拒绝域为\(W=\{z≤z_\alpha\}\)，于是观测值\(z=-1.56\)落在拒绝域里，应拒绝原假设； 由此可以看出，0.0599是能用观测值\(z=-1.56\)做出“拒绝\(H_0\)”的最小的显著性水平。

# 使用检验函数`z.test()`和`t.test()`
library(BSDA)
# 单样本均值的单边检验
z.test(battery, mu = 20, sigma.x = 0.8, alternative = "less", conf.level = 0.95)

## 
##  One-sample z-Test
## 
## data:  battery
## z = -1.5556, p-value = 0.0599
## alternative hypothesis: true mean is less than 20
## 95 percent confidence interval:
##        NA 20.01009
## sample estimates:
## mean of x 
##    19.824

# 单样本均值的单边检验（方差未知）
t.test(battery, mu = 20, alternative = "less", conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  battery
## t = -0.88266, df = 49, p-value = 0.1909
## alternative hypothesis: true mean is less than 20
## 95 percent confidence interval:
##     -Inf 20.1583
## sample estimates:
## mean of x 
##    19.824

使用R中自带的检验函数，在总体方差已知的情况下，检验p值为0.0599；而在总体方差未知的情况下，检验的p值为0.1909。

红楼梦虚词频率

假设小说对文言虚词的使用频率服从正态分布，且词频方差在不同阶段保持不变，在显著性水平\(\alpha=0.05\)下，尝试判断41-80回中“亦”使用的频数均值是否比81-120回中高？

# 红楼梦
wenyan <- read.csv("./data/红楼梦虚词词频统计.csv", stringsAsFactors = F, fileEncoding = "UTF-8")
freqx <- wenyan$`亦`[41:80]
freqy <- wenyan$`亦`[81:120]
t.test(freqx, freqy, alternative = "greater", var.equal = T, conf.level = 0.95)

## 
##  Two Sample t-test
## 
## data:  freqx and freqy
## t = 4.4123, df = 78, p-value = 1.621e-05
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  1.961594      Inf
## sample estimates:
## mean of x mean of y 
##      3.90      0.75

由于检验的P值小于显著性水平0.05，因此拒绝原假设，认为41-80回“亦”使用的平均数量显著高于81-120回，后期语言向白话文靠拢。

狗熊制药成对样本检验

假设狗熊会开了一个狗熊制药公司，专业研究减肥药。如果这个药品要在美国市场上合法上市，必须有美国食品药品管理局（FDA）的批准。为此，FDA对使用该药物的受试者做了一组试验，记录了其在使用减肥药前后的体重（单位：kg）。试验的数据分别为

受试者编号	1	2	3	4	5	6	7	8	9	10
服用药物前的体重x	43	55	49	62	59	49	57	54	55	48
服用药物后的体重y	50	59	55	60	58	54	56	53	61	51
差 d = x-y	-7	-4	-6	2	1	-5	1	1	-6	- 3

的狗熊制药案例数据中，假定受试者的体重数据服从正态分布，试问：试验前后的受试者体重在显著性水平\(\alpha=0.05\)上有无差异？

# 成对数据检验
x <- c(43,55,49,62,59,49,57,54,55,48)
y <- c(50,59,55,60,58,54,56,53,61,51)
# 成对数据t检验
t.test(x, y, alternative = "two.sided", paired = T, var.equal = F, conf.level = 0.95)

## 
##  Paired t-test
## 
## data:  x and y
## t = -2.3475, df = 9, p-value = 0.04348
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.10545182 -0.09454818
## sample estimates:
## mean of the differences 
##                    -2.6

检验的p值为0.0435，故应拒绝原假设\(H_0:\mu=0\)，即可认为试验前后的受试者体重有显著差异。进一步，平均含量差值的估计量为\(\hat \mu = \bar x - \bar y = -2.6\)，可见服用药物后的体重要高于服用药物前的体重。

习题答案

题目 5.2

假设某产品在制造时分为合格品与不合格品两类，我们用一个随机变量X来表示某个产品经检查后是否合格，X=0表示合格品，X=1表示不合格品，则X服从两点分布\(b(1,p)\)，其中p是未知的不合格品率。现抽取20个产品，看其是否合格，得到样本如下所示：

1 0 0 0 1 0 1 1 1 0 1 0 0 1 0 0 1 1 1 1

求p的最大似然估计，并在R中进行实现。

set.seed(4)
p = 0.5
# 生成的样本是
samples = rbinom(n=20, size=1, prob = p)
print(samples)

##  [1] 1 0 0 0 1 0 1 1 1 0 1 0 0 1 0 0 1 1 1 1

# p 的最大似然估计值是
p_mle = mean(samples)
cat('均值的最大似然估计值是: ', p_mle, '\n')

## 均值的最大似然估计值是:  0.55

均值的最大似然估计值是:0.55。

题目 5.3

假设灯泡的寿命服从正态分布。为了估计某种灯泡的平均寿命，现随机地抽取15只灯泡测试，得到它们的寿命（单位：小时）如下：

998 1021 988 986 1018 997 996

973 925 985 981 1007 1011 968 1002

试求灯泡平均寿命的0.95置信区间，并在R中进行实现。

set.seed(8)
n = 15
u = 1000
samples = rnorm(n, u, 25)
samples = round(samples, 0)
cat('寿命样本: ', samples, '\n')

## 寿命样本:  998 1021 988 986 1018 997 996 973 925 985 981 1007 1011 968 1002

# 区间估计
# 方法一：按照公式计算
weight_u = mean(samples)
weight_var = sum((samples - weight_u) ^ 2) / (n-1)
spread = qt(1-0.05/2, n-1) * sqrt(weight_var / n)

cat(sprintf('灯泡寿命的置信区间是: [%s, %s]', 
            round(weight_u - spread, 4),
            round(weight_u + spread, 4), '\n'))

## 灯泡寿命的置信区间是: [977.2537, 1003.5463]

# 方法二：使用t.test()计算
interval_2 <- t.test(samples, conf.level = 0.95)$conf.int
cat(sprintf('灯泡寿命的置信区间是: [%s, %s]', 
            round(interval_2[1], 4),
            round(interval_2[2], 4), '\n'))

## 灯泡寿命的置信区间是: [977.2537, 1003.5463]

无论使用公式还是使用R中的自带函数计算，灯泡寿命的置信区间都是: [977.2537, 1003.5463]。