31 R初等统计分析

这一部分讲授如何用R进行统计分析， 包括基本概括统计和探索性数据分析， 置信区间和假设检验， 回归分析与各种回归方法， 广义线性模型， 非线性回归与平滑， 判别树和回归树， 等等。

主要参考书：

  • (Venables and Ripley 2002)
  • (Kabacoff 2012)

31.1 概率分布

R中与xxx分布有关的函数包括：

  • dxxx(x)， 即xxx分布的分布密度函数(PDF)或概率函数(PMF)\(p(x)\)。
  • pxxx(q)， 即xxx分布的分布函数(CDF)\(F(q)=P(X \leq q)\)。
  • qxxx(p)， 即xxx分布的分位数函数\(q(p)\), \(p \in (0,1)\)， 对连续型分布，\(q(p) = F^{-1}(p)\)， 即\(F(x)=p\)的解\(x\)。
  • rxxx(n)， 即xxx的随机数函数，可以生成\(n\)个xxx的随机数。

dxxx(x)函数可以加选项log=TRUE， 用来计算\ln p(x)， 这在计算对数似然函数时有用， 比log(dxxx(x))更精确。

pxxx(q)可以加选项lower.tail=FALSE， 变成计算\(P(X>q) = 1 - F(q)\)。

qxxx(p)可以加选项lower.tail=TRUE， 表示求\(P(X>x)=p\)的解\(x\)； 可以加选项log.p=TRUE， 表示输入的\(p\)实际是\(\ln p\)。

这些函数都可以带有自己的分布参数， 有些分布参数有缺省值， 比如正态分布的缺省参数值为零均值单位标准差。

具体的分布类型可以在R命令行用?Distributions查看列表。 常用的分布密度有：

  • 离散分布有dbinom二项分布, dpois泊松分布， dgeom几何分布， dnbinom负二项分布， dmultinom多项分布， dhyper超几何分布。
  • 连续分布有 dunif均匀分布， dnorm正态分布， dchisq卡方分布， dt t分布(包括非中心t)， df F分布， dexp指数分布， dweibull 威布尔分布， dgamma 伽马分布， dbeta 贝塔分布， dlnorm 对数正态分布， dcauchy 柯西分布, dlogis 逻辑斯谛分布。

R的扩展包提供了更多的分布， 参见R网站的如下链接：

31.2 最大似然估计

R函数optim()、nlm()、optimize()可以用来求函数极值， 因此可以用来计算最大似然估计。 optimize()只能求一元函数极值。

31.2.1 一元正态分布参数最大似然估计

正态分布最大似然估计有解析表达式。 作为示例， 用R函数进行数值优化求解。

对数似然函数为：

\[ \ln L(\mu,\sigma^2) = -\frac{n}{2}\ln(2\pi) -\frac{n}{2}\ln\sigma^2 - \frac{1}{2\sigma^2} \sum(X_i - \mu)^2 \]

定义R的优化目标函数为上述对数似然函数去掉常数项以后乘以\(-2\)， 求其最小值点。 目标函数为：

objf.norm1 <- function(theta, x){
  mu <- theta[1]
  s2 <- exp(theta[2])
  n <- length(x)
  res <- n*log(s2) + 1/s2*sum((x - mu)^2)
  res
}

其中\(\theta_1\)为均值参数\(\mu\)， \(\theta_2\)为方差参数\(\sigma^2\)的对数值。 x是样本数值组成的R向量。

可以用optim函数来求极小值点。下面是一个模拟演示:

norm1d.mledemo1 <- function(n=30){
  mu0 <- 20
  sigma0 <- 2
  set.seed(1)
  x <- rnorm(n, mu0, sigma0)

  theta0 <- c(0,0)
  ores <- optim(theta0, objf.norm1, x=x)
  print(ores)
  theta <- ores$par
  mu <- theta[1]
  sigma <- exp(0.5*theta[2])
  cat('真实mu=', mu0, ' 公式估计mu=', mean(x), 
      ' 数值优化估计mu=', mu, '\n')
  cat('真实sigma=', sigma0, 
      '公式估计sigma=', sqrt(var(x)*(n-1)/n), 
      ' 数值优化估计sigma=', sigma, '\n')
}
norm1d.mledemo1()
## $par
## [1] 20.166892  1.193593
## 
## $value
## [1] 65.83709
## 
## $counts
## function gradient 
##      115       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## 真实mu= 20  公式估计mu= 20.16492  数值优化估计mu= 20.16689 
## 真实sigma= 2 公式估计sigma= 1.817177  数值优化估计sigma= 1.816291

也可以用nlm()函数求最小值点， 下面是正态分布最大似然估计的另一个演示：

norm1d.mledemo2 <- function(){
  set.seed(1)
  n <- 30
  mu <- 20
  sig <- 2
  z <- rnorm(n, mean=mu, sd=sig)
  neglogL <- function(parm) {
      -sum( dnorm(z, mean=parm[1],
                  sd=exp(parm[2]), log=TRUE) )
  }

  res <- nlm(neglogL, c(10, log(10)))
  print(res)
  sig2 <- exp(res$estimate[2])
  cat('真实mu=', mu, ' 公式估计mu=', mean(z), 
      ' 数值优化估计mu=', res$estimate[1], '\n')
  cat('真实sigma=', sig, 
      '公式估计sigma=', sqrt(var(z)*(n-1)/n), 
      ' 数值优化估计sigma=', sig2, '\n')
  invisible()
}
norm1d.mledemo2()
## $minimum
## [1] 60.48667
## 
## $estimate
## [1] 20.1649179  0.5972841
## 
## $gradient
## [1] 1.444576e-05 9.094805e-06
## 
## $code
## [1] 2
## 
## $iterations
## [1] 28
## 
## 真实mu= 20  公式估计mu= 20.16492  数值优化估计mu= 20.16492 
## 真实sigma= 2 公式估计sigma= 1.817177  数值优化估计sigma= 1.817177

函数optim()缺省使用Nelder-Mead单纯型搜索算法， 此算法不要求计算梯度和海色阵， 算法稳定性好， 但是收敛速度比较慢。 可以用选项method="BFGS"指定BFGS拟牛顿法。 这时可以用gr=选项输入梯度函数， 缺省使用数值微分计算梯度。 如：

norm1d.mledemo1b <- function(n=30){
  mu0 <- 20
  sigma0 <- 2
  set.seed(1)
  x <- rnorm(n, mu0, sigma0)

  theta0 <- c(1,1)
  ores <- optim(theta0, objf.norm1, x=x, method="BFGS")
  print(ores)
  theta <- ores$par
  mu <- theta[1]
  sigma <- exp(0.5*theta[2])
  cat('真实mu=', mu0, ' 公式估计mu=', mean(x), 
      ' 数值优化估计mu=', mu, '\n')
  cat('真实sigma=', sigma0, 
      '公式估计sigma=', sqrt(var(x)*(n-1)/n), 
      ' 数值优化估计sigma=', sigma, '\n')
}
norm1d.mledemo1b()
## $par
## [1] 20.164916  1.194568
## 
## $value
## [1] 65.83704
## 
## $counts
## function gradient 
##       41       17 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## 真实mu= 20  公式估计mu= 20.16492  数值优化估计mu= 20.16492 
## 真实sigma= 2 公式估计sigma= 1.817177  数值优化估计sigma= 1.817177

注意使用数值微分时， 参数的初值的数量级应该与最终结果相差不大， 否则可能不收敛， 比如上面程序取初值\((0,0)\)就会不收敛。

optim()函数支持的方法还有"CG"是一种共轭梯度法， 此方法不如BFGS稳健， 但是可以减少存储量使用， 对大规模问题可能更适用。 "L-BFGS-B"是一种区间约束的BFGS优化方法。 "SANN"是模拟退火算法， 有利于求得全局最小值点， 速度比其它方法慢得多。

函数optimize()进行一元函数数值优化计算。 例如，在一元正态分布最大似然估计中， 在对数似然函数中代入\(\mu = \bar x\)，目标函数（负2倍对数似然函数）为 \[ h(\sigma^2) = \ln(\sigma^2) + \frac{1}{\sigma^2} \sum (X_i - \bar X)^2 \] 下面的例子模拟计算\(\sigma^2\)的最大似然估计：

norm1d.mledemo3 <- function(n=30){
  mu0 <- 20
  sigma0 <- 2
  set.seed(1)
  x <- rnorm(n, mu0, sigma0)

  mu <- mean(x)
  ss <- sum((x - mu)^2)/length(x)
  objf <- function(delta, ss) log(delta) + 1/delta*ss
  ores <- optimize(objf, lower=0.0001,
                   upper=1000, ss=ss)
  delta <- ores$minimum
  sigma <- sqrt(delta)

  print(ores)
  cat('真实sigma=', sigma0, 
      '公式估计sigma=', sqrt(var(x)*(n-1)/n), 
      ' 数值优化估计sigma=', sigma, '\n')
}
norm1d.mledemo3()
## $minimum
## [1] 3.30214
## 
## $objective
## [1] 2.194568
## 
## 真实sigma= 2 公式估计sigma= 1.817177  数值优化估计sigma= 1.817179

31.3 假设检验和置信区间

31.3.1 均值的假设检验和置信区间

31.3.1.1 单样本均值

31.3.1.1.1 大样本情形

设总体\(X\)期望为\(\mu\)，方差为\(\sigma^2\)。 \(X_1, X_2, \dots, X_n\)是一组样本。

在大样本时，令 \[ Z = \frac{\bar x - \mu}{S / \sqrt{n}} \stackrel{\bullet}{\sim} \text{N}(0,1) \] 其中\(\stackrel{\bullet}{\sim}\)表示近似服从某分布。 \(\bar x\)为样本均值，\(S\)为样本标准差。 \(\sigma\)未知且大样本时， \(\mu\)近似\(1-\alpha\)置信区间为 \[ \bar x \pm \text{qnorm}(1-\alpha/2) S / \sqrt{n} \]

用BSDA::z.test(x, sigma.x=sd(x), conf.level=1-alpha)可以计算如上的近似置信区间。

对于假设检验问题 \[\begin{aligned} & H_0: \mu = \mu_0 \longleftrightarrow H_a: \mu \neq \mu_0 \\ & H_0: \mu \leq \mu_0 \longleftrightarrow H_a: \mu > \mu_0 \\ & H_0: \mu \geq \mu_0 \longleftrightarrow H_a: \mu < \mu_0 \end{aligned}\] 定义检验统计量 \[ Z = \frac{\bar x - \mu_0}{S / \sqrt{n}} \] 当\(\mu=\mu_0\)时\(Z\)近似服从N(0,1)分布。 可以用BSDA::z.test(x, mu=mu0, sigma.x=sd(x))进行双侧Z检验。 加选项alternative="greater"作右侧检验， 加选项alternative="less"作左侧检验。

如果上述问题中标准差\(\sigma\)为已知， 应该在\(Z\)的公式用\(\sigma\)替换\(S\)，即 \[ Z = \frac{\bar x - \mu_0}{\sigma / \sqrt{n}} \] 计算\(\mu\)的置信区间程序为BSDA::z.test(x, sigma.x=sigma0)， 计算\(Z\)检验的程序为BSDA::z.test(x, mu=mu0, sigma.x=sigma0)， 仍可用选项alternative=指定双侧、右侧、左侧。

31.3.1.1.2 小样本正态情形

如果样本量较小（小于30）， 则需要假定总体服从正态分布N(\(\mu, \sigma^2\))。

当\(\sigma\)已知时， 有 \[ Z = \frac{\bar x - \mu_0}{\sigma / \sqrt{n}} \sim \text{N}(0,1) \] 计算\(\mu\)的置信区间程序为BSDA::z.test(x, sigma.x=sigma0)， 计算\(Z\)检验的程序为BSDA::z.test(x, mu=mu0, sigma.x=sigma0)。

可以自己写一个计算Z检验的较通用的R函数：

z.test.1s <- function(
  x, n=length(x), mu=0, sigma=sd(x), alternative="two.sided"){
  z <- (mean(x) - mu) / (sigma / sqrt(n))
  if(alternative=="two.sided"){   # 双侧检验
    pvalue <- 2*(1 - pnorm(abs(z)))
  } else if(alternative=="less"){ # 左侧检验
    pvalue <- pnorm(z)
  } else if(alternative=="greater"){ # 右侧检验
    pvalue <- 1 - pnorm(z)
  } else {
    stop("alternative unknown!")
  }
  
  c(stat=z, pvalue=pvalue)
}

这个函数输入样本x和\(\sigma=\sigma_0\)， 或者输入\(\bar x\)为x和\(\sigma=\sigma_0\)， 输出检验统计量值和p值。

如果\(\sigma\)未知，应使用t统计量检验 \[ T = \frac{\bar x - \mu}{S / \sqrt{n}} \] 当\(\mu=\mu_0\)时\(T \sim t(n-1)\)。 如果x保存了样本， 用t.test(x, conf.level=1-alpha)计算\(\mu\)的置信区间； 用t.test(x, mu=mu0, alternative="two.side")计算\(H_0: \mu=\mu_0\)的双侧检验， 称为单样本t检验。

如果已知的\(\bar x\)和\(S\)而非具体样本数据， 可以自己写一个R函数计算t检验：

t.test.1s <- function(
  x, n=length(x), sigma=sd(x), mu=0, 
  alternative="two.sided"){
  tstat <- (mean(x) - mu) / (sigma / sqrt(n))
  if(alternative=="two.sided"){   # 双侧检验
    pvalue <- 2*(1 - pt(abs(tstat), n-1))
  } else if(alternative=="less"){ # 左侧检验
    pvalue <- pt(tstat, n-1)
  } else if(alternative=="greater"){ # 右侧检验
    pvalue <- 1 - pt(tstat, n-1)
  } else {
    stop("alternative unknown!")
  }
  
  c(stat=tstat, pvalue=pvalue)
}

这个函数允许输入\(\bar x\)到x中，输入\(n\)，输入\(S\)到sigma中， 然后计算t检验。

31.3.1.1.3 统计假设检验显著性的解释

统计假设检验如果拒绝了\(H_0\)， 也称“结果显著”、“有显著差异”等。 这只代表样本与零假设的差距足够大， 可以以比较小的代价（第一类错误概率） 排除这样的差距是由于随机样本的随机性引起的， 但是并不表明这样的差距是有科学意义或者有现实意义的。

尤其是超大样本时， 完全没有实际意义的差距也会检验出来是统计显著的。

这时，可以设定一个“有实际意义的差距”\(\delta\)， 当差距超过\(\delta\)时才认为有实际差距，如： \[ H_0: |\mu - \mu_0| \leq \delta \longleftrightarrow |\mu - \mu_0| > \delta \]

统计假设检验如果不拒绝\(H_0\)， 也称为“结果不显著”、“没有显著差异”， 通常不认为\(H_0\)就是对的， 看法是没有充分理由推翻\(H_0\)。

进一步收集更多的数据， 有可能就会发现显著差异。

例如，\(H_0\)是嫌疑人无罪， 最后不拒绝\(H_0\)， 只能是没有充分证据证明有罪， 如果后来收集到了关键的罪证， 就还可以再判决嫌疑人有罪。

31.3.1.1.4 置信区间与双侧检验的关系

设总体为正态分布\(\text{N}(\mu, \sigma^2)\), \(\sigma\)已知。 随机抽取了\(n\)个样品，计算得\(\bar x\)。 \(\mu\)的\(1-\alpha\)置信区间为 \[ \bar x \pm z_{\alpha/2} \sigma/\sqrt{n} \]

置信区间理论依据： \[ P\left(\frac{|\bar x - \mu|}{\sigma/\sqrt{n}} \leq z_{\alpha/2} \right) = 1 - \alpha \tag{*} \]

对双侧检验 \[ H_0: \mu = \mu_0 \longleftrightarrow H_a: \mu \neq \mu_0 \] 当 \[ \frac{|\bar x - \mu_0|}{\sigma/\sqrt{n}} \leq z_{\alpha/2}\] 时接受\(H_0\)。

接受\(H_0\)的条件，就是\(\mu_0\)满足(*)式的大括号中的条件， 即当且仅当\(\mu_0\)落入\(\mu\)的\(1-\alpha\)置信区间时接受\(H_0\)。

其它检验也与置信区间有类似关系， 但是单侧检验对应于单侧置信区间。

31.3.1.1.5 单样本均值比较例子

载入假设检验部分示例程序所用数据：

load("hyptest-data.RData")
31.3.1.1.5.1 咖啡标重的单侧检验

美国联邦贸易委员会定期检查商品标签是否与实际相符。 Hiltop Coffee大罐咖啡标签注明含量为3磅。 每罐具体多一些或少一些并不要紧， 只要总体平均值不少于3磅， 消费者利益就未受侵害。

只有在有充分证据时才应该做出分量不够的判决进而对厂家进行处罚。 所以零假设是\(H_0: \mu \geq 3\)， 对立假设是\(H_a: \mu < 3\)。

收集数据进行统计判决后， 如果不拒绝\(H_0\)， 就不需要采取行动； 如果拒绝了\(H_0\)， 就说明有足够证据认为厂商灌装分量不足， 应进行处罚。

随机抽取了36件样品， 计算平均净重\(\bar x\)作为总体均值\(\mu\)的估计。

如果\(\bar x < 3\)， 则\(H_0\)有理由被怀疑。

\(\bar x\)比\(\mu_0 = 3\)少多少， 才让我们愿意拒绝\(H_0\)? 做出拒绝\(H_0\)的决定可能会犯第一类错误， 如果差距很小， 这样的差距很可能是随机抽样的随机性引起的， 拒绝\(H_0\)有很大的机会犯错。

必须先确定一个能够容忍的第一类错误概率。 调查员认为，可以容忍1%的错误拒绝\(H_0\)的机会 （不该处罚时给出了处罚） 以换得能够在厂商分量不足时正确地做出处罚决定的机会。 即选检验水平\(\alpha=0.01\)。

当\(H_0\)实际成立时， 边界\(\mu=\mu_0=3\)处犯第一类错误的概率最大 （这里最不容易区分）， 所以计算第一类错误概率只要在边界处计算。

统计假设检验第一步是确定零假设与对立假设， 第二步是确定检验水平\(\alpha\)。 第三步是抽取随机样本，计算检验统计量。

以往经验证明净重的标准差是\(\sigma=0.18\)， 且净重的分布为正态分布。

则\(\bar x\)的抽样分布为\(\text{N}(\mu, \sigma^2/n)\)， \(\bar x\)的标准误差为\(\text{SE}(\bar x) = 0.18 / \sqrt{36} = 0.03\)。

令 \[ z = \frac{\bar x - \mu_0}{\text{SE}(\bar x)} \] 则\(H_0\)成立且\(\mu\)取边界处的\(\mu_0\)时\(z \sim \text{N}(0,1)\)。

当\(\bar x\)比\(\mu_0=3\)小很多时应拒绝\(H_0\)。 即检验统计量\(z\)很小时拒绝\(H_0\)。

如果\(z = -1\)， 标准正态分布取-1或比-1更小的值的概率是0.1586 (用R程序pnorm(-1)计算)， 这个概率不小， 不小于\(\alpha=0.01\)， 这时不应该拒绝\(H_0\)。

如果\(z = -2\)， 标准正态分布取-2或比-2更小的值的概率是0.0228， (用R程序pnorm(-2)计算)， 这个概率已经比较小了， 但还不小于预先确定的第一类错误概率\(\alpha=0.01\)， 这时不应该拒绝\(H_0\)。

如果\(z = -3\)， 标准正态分布取-3或比-3更小的值的概率是0.0013， (用R程序pnorm(-3)计算)， 这个概率已经很小了， 也小于预先确定的第一类错误概率\(\alpha=0.01\)， 这时应该拒绝\(H_0\)。

对左侧检验问题\(H_0: \mu \geq \mu_0 \longleftrightarrow H_a: \mu < \mu_0\), 若\(Z\)表示标准正态分布随机变量， \(z\)是检验统计量的值，在\(\mu=\mu_0\)时\(z\)的抽样分布是标准正态分布， 可以计算抽样分布中取到等于\(z\)以及比\(z\)更小的值的概率\(P(Z \leq z)\)， 称为检验的p值。

p值越小，说明\(\bar x\)比\(\mu_0\)小得越多， 拒绝\(H_0\)也越有充分证据。 当且仅当p值小于检验水平\(\alpha\)时拒绝\(H_0\)。

p值小于\(\alpha\)，说明如果\(H_0\)真的成立， 出现\(\bar x\)那么小的检验统计量值的概率也很小， 小于\(\alpha\)，所以第一类错误概率小于等于\(\alpha\)。

一般地，p值是从样本计算的一个在假定零假设成立情况下， 检验统计量取值更倾向于对立假设成立的概率。 p值越小，拒绝\(H_0\)越有依据。

设咖啡标称重量问题中抽取了36件样品，测得\(\bar x = 2.92\)磅， 计算得 \[ z = \frac{\bar x - \mu_0}{\sigma/\sqrt{n}} = \frac{2.92 - 3}{0.18/\sqrt{36}} = -2.67 \]

p值为\(P(Z \leq -2.67)\)(\(Z\)是标准正态分布随机变量)， 用R计算得pnorm(-2.67)\(=0.0038 \leq \alpha=0.01\)， 拒绝\(H_0\)， 应对该厂商进行处罚。

已经预先确定了检验水平\(\alpha=0.01\)就不能在看到检验结果后更改检验水平。 但是， 报告出p值后， 其他的人可以利用这样的信息， 不同的人或利益方可能有不同的检验水平的选择， 这通过同一个p值都可以给出需要的结果。

p值也称为“观测显著性水平”， 因为它是能够拒绝\(H_0\)的最小的可选检验水平， 选取比p值更小的检验水平就不能拒绝\(H_0\)了。

用自定义的z.test.1s()计算：

z.test.1s(x=2.92, mu=3, n=36, sigma=0.18, alternative="less")
##         stat       pvalue 
## -2.666666667  0.003830381

p值为0.004，在0.01水平下拒绝原假设， 认为咖啡平均重量显著地低于标称的3磅。

设Coffee数据框的Weight是这36个样本值， 也可以将程序写成：

z.test.1s(Coffee[["Weight"]], mu=3, sigma=0.18, alternative="less")
##         stat       pvalue 
## -2.666666667  0.003830381

或

BSDA::z.test(Coffee[["Weight"]], mu=3, sigma.x=0.18, alternative="less")
## 
##  One-sample z-Test
## 
## data:  Coffee[["Weight"]]
## z = -2.6667, p-value = 0.00383
## alternative hypothesis: true mean is less than 3
## 95 percent confidence interval:
##        NA 2.969346
## sample estimates:
## mean of x 
##      2.92
31.3.1.1.5.2 例：高尔夫球性能的双侧检验

美国高尔夫协会对高尔夫运动器械有一系列的规定。 MaxFlight生产高品质的高尔夫球， 平均的飞行距离是295码。

生产线有时会有波动， 当生产的球平均飞行距离不足295码时， 用户会感觉不好用； 当生产的球平均飞行距离超过295码时， 美国高尔夫协会禁止使用这样的球。

该厂商定期抽检生产线上下来的产品， 每次抽取50只。

当平均飞行距离与295码没有显著差距时， 不需要采取措施； 如果有显著差距，就需要调整生产线。

所以取\(H_0: \mu = 295\), \(H_a: \mu \neq 295\)。 选检验水平\(\alpha=0.05\)。

随机抽取\(n\)个样品， 在可控条件下测试出每个的飞行距离， 计算出平均值\(\bar x\)。

当\(\bar x\)比\(\mu_0\)(这里是295)小很多或者 \(\bar x\)比\(\mu_0\)大很多时拒绝\(H_0\)， 否则不拒绝\(H_0\)。

设已知总体标准差\(\sigma=12\), 且飞行距离服从正态分布。

取统计量 \[z = \frac{\bar x - \mu_0}{\sigma/\sqrt{n}}\] 在\(H_0\)成立时， \(z\)的抽样分布为标准正态分布。

  • 设测得\(\bar x = 297.6\)，超过了\(\mu_0=295\)。

\[ z = = \frac{\bar x - \mu_0}{\sigma/\sqrt{n}} = \frac{297.6 - 295}{12/\sqrt{50}} = 1.53 \]

\(H_0\)成立时\(z\)抽样分布为标准正态分布， 设\(Z\)是标准正态分布随机变量， \(Z\)取到1.53或者比1.53更倾向于反对\(H_0\)的值概率？ 这包括\(P(Z \geq 1.53)\)， 还应包括\(P(Z \leq -1.53)\)。

所以p值是\(P(|Z| \geq |z|)\)， R中计算为2*(1 - pnorm(abs(z)))。 p值等于0.1260， 超过检验水平， 不拒绝\(H_0\)， 不需要调整生产线。

用自定义的z.test.1s()检验：

z.test.1s(GolfTest[['Yards']], mu=295, sigma=12, alternative="two.sided")
##      stat    pvalue 
## 1.5320647 0.1255065

或用BSDA::z.test():

BSDA::z.test(GolfTest[['Yards']], mu=295, sigma.x=12, alternative="two.sided")
## 
##  One-sample z-Test
## 
## data:  GolfTest[["Yards"]]
## z = 1.5321, p-value = 0.1255
## alternative hypothesis: true mean is not equal to 295
## 95 percent confidence interval:
##  294.2738 300.9262
## sample estimates:
## mean of x 
##     297.6

p值为0.13，在0.05水平下不拒绝\(H_0\)， 生产的球的平均飞行距离与标准的295码之间没有显著差异。

31.3.1.1.5.3 Heathrow机场打分的检验

某旅行杂志希望根据乘机进行商务旅行的意见给跨大西洋门户机场评分。 乘客评分取0到10分，超过7分算作是服务一流的机场。 在每个机场随机选取60位商务旅行乘客打分评价。 英国Heathrow机场的打分样本\(\bar x = 7.25\), \(S = 1.052\)。 英国Heathrow机场可以算是服务一流吗？

由于随机样本的随机性， 要判断的是\(\mu > 7\)是否成立， 这里\(\bar x = 7.25\)并不能很确定第说\(\mu > 7\)成立。

因为将机场评为一流需要有充分证据， 所以将对立假设\(H_a\)设为\(\mu > 7\)， 然后取\(H_0: \mu \leq 7\)。

取检验水平\(\alpha=0.05\)。

计算检验统计量 \[ t = \frac{\bar x - \mu_0}{S / \sqrt{n}} = \frac{7.25 - 7}{1.052/\sqrt{60}} = 1.84 \] 设\(T\)是服从\(t(n-1)\)分布的随机变量， p值为\(P(T \geq 1.84)\)。 用R软件计算为1 - pt(1.84, 60-1)\(=0.0354\)。 p值小于\(\alpha\)，拒绝\(H_0\)， 认为平均评分显著高于7， 可以认为Heathrow机场是服务一流。

从原始样本数据用t.test()检验：

t.test(AirRating[["Rating"]], mu=7, alternative="greater")
## 
##  One Sample t-test
## 
## data:  AirRating[["Rating"]]
## t = 1.8414, df = 59, p-value = 0.03529
## alternative hypothesis: true mean is greater than 7
## 95 percent confidence interval:
##  7.023124      Inf
## sample estimates:
## mean of x 
##      7.25

用自定义的t.test.1s()和样本统计量计算：

t.test.1s(x=7.25, sigma=1.052, n=60, mu=7, alternative="greater")
##       stat     pvalue 
## 1.84077155 0.03534255
31.3.1.1.5.4 玩具厂商订货量的假设检验

玩具厂商Holiday Toys通过超过1000家玩具商店售货。 为了准备即将到来的冬季销售季节， 销售主管需要确定今年某新款玩具的生产量。 凭经验判断平均每家商店需要40件。

为了确定，随机抽取了25家商店， 提供了新款玩具的特点、进货价和建议销售价， 让他们给出订货量估计。 目的是判断按每家40件来生产是否合适。

因为多生产或少生产都不合适， 所以取零假设为\(H_0: \mu=40\)。 对立假设为\(H_a: \mu \neq 40\)。 只有拒绝\(H_0\)时才需要更改生产计划。

取检验水平\(\alpha=0.05\)。 样本计算得\(n=25\)，\(\bar x = 37.4\), \(S = 11.79\)。

检验统计量 \[ t = \frac{\bar x - \mu_0}{S / \sqrt{n}} = \frac{37.4 - 40}{11.79 / \sqrt{25}} = -1.10 \]

设\(T\)为服从t\((n-1)\)分布的随机变量， p值为\(2 P(|T| > |-1.10|)\)， 用R计算为2*(1 - pt(abs(-1.10), 25-1))\(=0.2822\)。

不拒绝\(H_0\)，所以不需要更改生产计划。

从原始样本数据用t.test()检验：

t.test(Orders[["Units"]], mu=40, alternative="two.sided")
## 
##  One Sample t-test
## 
## data:  Orders[["Units"]]
## t = -1.1026, df = 24, p-value = 0.2811
## alternative hypothesis: true mean is not equal to 40
## 95 percent confidence interval:
##  32.5334 42.2666
## sample estimates:
## mean of x 
##      37.4

用自定义的t.test.1s()和样本统计量计算：

t.test.1s(x=37.4, sigma=11.79, n=25, mu=40, alternative="two.sided")
##       stat     pvalue 
## -1.1026293  0.2811226

31.3.1.2 均值比较

31.3.1.2.1 独立两样本Z检验

假设有两个独立的总体\(X, Y\)，例如，男性与女性的寿命。 要比较两个总体的期望\(\mu_1=EX\)和\(\mu_2=EY\)。 设两个总体的方差分别为\(\sigma_1^2\)和\(\sigma_2^2\)。 可以对均值作双侧、左侧、右侧检验。 比如，双侧检验为 \[ H_0: \mu_1 = \mu_2 \longleftrightarrow H_a: \mu_1 \neq \mu_2 \]

如果两个总体都服从正态分布， 而且分别的方差\(\sigma_1^2\)和\(\sigma_2^2\)已知， 设\(X\)的样本量为\(n_1\)的样本得到样本均值\(\bar x\)， 设\(Y\)的样本量为\(n_2\)的样本得到样本均值\(\bar y\)， 取统计量 \[ Z = \frac{\bar x - \bar y} {\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}} \] 其中分母是分子的标准误差。 当\(\mu_1 = \mu_2\)时\(Z \sim \text{N}(0,1)\)。

如果分布不一定是正态分布但是样本量足够大， 可以令 \[ Z = \frac{\bar x - \bar y} {\sqrt{\frac{S_x^2}{n_1} + \frac{S_y^2}{n_2}}} \] 当\(\mu_1 = \mu_2\)时\(Z\)近似服从\(\text{N}(0,1)\)。

检验问题还可以变成\(\mu_1 - \mu_2\)与某个\(\delta\)的比较，如 \[ H_0: \mu_1 \leq \mu_2 - \delta \longleftrightarrow H_a: \mu_1 > \mu_2 - \delta \] 这个对立假设是“不次于”的结论。

可以写一个大样本方差未知或已知， 小样本独立两正态总体方差已知情况做Z检验的R函数：

z.test.2s <- function(
  x, y, n1=length(x), n2=length(y), delta=0,
  sigma1=sd(x), sigma2=sd(y), alternative="two.sided"){
  z <- (mean(x) - mean(y) - delta) / sqrt(sigma1^2 / n1 + sigma2^2 / n2)
  if(alternative=="two.sided"){   # 双侧检验
    pvalue <- 2*(1 - pnorm(abs(z)))
  } else if(alternative=="less"){ # 左侧检验
    pvalue <- pnorm(z)
  } else if(alternative=="greater"){ # 右侧检验
    pvalue <- 1 - pnorm(z)
  } else {
    stop("alternative unknown!")
  }
  
  c(stat=z, pvalue=pvalue)
}

函数允许输入两组样本到x和y中， 输入\(\sigma_1\)和\(\sigma_2\)到sigma1和sigma2中， 这时不输入sigma1和sigma2则从样本中计算样本统计量代替。 也允许输入\(\bar x\)和\(\bar y\)到x和y中， 输入\(\sigma_1\)和\(\sigma_2\)到sigma1和sigma2中， 或者输入S_x$和$S_y$到sigma1和sigma2中， 输入$n_1$到n1中， 输入$n_2$到n2`中， 计算独立两样本均值比较的\(Z\)检验。

31.3.1.2.2 独立两样本t检验

如果样本是小样本，但两个总体分别服从正态分布， 两个总体的方差未知， 但已知\(\sigma_1^2 = \sigma_2^2\)。 这时\(\sigma_1^2 = \sigma_2^2\)可以统一估计为 \[ S_p^2 = \frac{1}{n_1 + n_2 - 2}\left( (n_1 - 1)S_x^2 + (n_2 - 1)S_y^2 \right) \] 取检验统计量 \[ T = \frac{\bar x - \bar y}{S_p \sqrt{\frac{1}{n_1} + \frac{1}{n_2}}} \] 当\(\mu_1 = \mu_2\)时\(T \sim \text{t}(n_1 + n_2 - 2)\)。

当x和y存放了两组独立样本时， 可以用t.test(x, y, var.equal=TRUE)计算两样本t检验， 用alternative=选项指定双侧、右侧、左侧检验。

可以自己写一个这样的R函数作两样本t检验， 允许只输入统计量而非具体样本值：

t.test.2s <- function(
  x, y, n1=length(x), n2=length(y),
  sigma1=sd(x), sigma2=sd(y), delta=0,
  alternative="two.sided"){
  sp <- sqrt(1/(n1+n2-2) * ((n1-1)*sigma1^2 + (n2-1)*sigma2^2))
  t <- (mean(x) - mean(y) - delta) / (sp * sqrt(1 / n1 + 1 / n2))
  if(alternative=="two.sided"){   # 双侧检验
    pvalue <- 2*(1 - pt(abs(t), n1+n2-2))
  } else if(alternative=="less"){ # 左侧检验
    pvalue <- pt(t, n1+n2-2)
  } else if(alternative=="greater"){ # 右侧检验
    pvalue <- 1 - pt(t, n1+n2-2)
  } else {
    stop("alternative unknown!")
  }
  
  c(stat=t, pvalue=pvalue)
}

如果两个独立正态总体的方差不相等，样本量不够大， 可以用如下的在\(H_0\)下近似服从t分布的检验统计量： \[ T = \frac{\bar x - \bar y}{\sqrt{ \frac{S_x^2}{n_1} + \frac{S_y^2}{n_2} }} \] 自由度为 \[ m = \frac{\left( S_x^2 / n_1 + S_y^2 / n_2 \right)^2} { \frac{\left( S_x^2 / n_1\right)^2}{n_1 - 1} + \frac{\left(S_y^2 / n_2 \right)^2}{n_2 - 1} } \] 这个检验称为Welch两样本t检验，或Satterthwaite 两样本t检验 当x和y存放了两组独立样本时， 可以用t.test(x, y, var.equal=FALSE)计算Welch两样本t检验。

如果方差是否相等不易判断， 可以直接选用Welch检验。

因为大样本时t分布与标准正态分布基本相同， 所以非正态总体大样本Z检验也可以用t.test()计算。

31.3.1.2.3 成对均值比较问题

设总体包含两个\(X, Y\)分量， 比如， 一组病人的服药前血压与服药后血压。 这两个变量是相关的， 不能用独立两总体的均值比较方法比较\(\mu_1=EX\)与\(\mu_2=EY\)。

若\(\xi = X - Y\)服从正态分布， 则可以对\(\xi\)进行均值为零的单样本t检验。 这样的检验称为成对均值的t检验。

若x和y分别为\(X\)和\(Y\)的样本值， 可以用t.test(x, y, paired=TRUE)计算计算成对均值t检验。

31.3.1.2.4 均值比较的例子
31.3.1.2.4.1 顾客平均年龄差别比较

HomeStyle是一个连锁家具店， 两家分店一个在城区，一个在郊区。 经理发现同一种家具有可能在一个分店卖的很好， 另一个分店则不好。 怀疑是顾客人群的人口学差异（年龄、性别、种族等）。

独立地抽取两个分店的顾客样本，比较平均年龄。 设城内的分店顾客年龄\(\sigma_1=9\)已知， 郊区的分店顾客年龄\(\sigma_2=10\)已知。

城内调查了\(n_1=36\)位顾客，年龄的样本平均值\(\bar x=40\); 郊区调查了\(n_2=49\)位顾客，年龄的样本平均值\(\bar y=35\)。 \(\bar x - \bar y = 5\)。

作双侧检验，水平\(\alpha=0.05\)。

计算\(Z\)统计量： \[ Z = \frac{40 - 35}{\sqrt{ \frac{9^2}{36} + \frac{10^2}{49} }} = 2.4138 \]

p值为\(P(|V| \geq |2.4138|)\)，其中\(V\)为标准正态分布随机变量。 用R计算为2*(1 - pnorm(abs(2.4138)))=0.016， 两个分店的顾客年龄有显著差异， 城里的顾客平均年龄显著地高。

基于样本统计量用自定义的z.test.2s()计算：

z.test.2s(n1=36, x=40, sigma1=9, 
          n2=49, y=35, sigma2=10,
          alternative="two.sided")
##       stat     pvalue 
## 2.41379310 0.01578742
31.3.1.2.4.2 两个银行营业所顾客平均存款比较

某银行要比较两个营业所的支票账户的平均存款额。 作双侧检验，水平0.05。没有方差信息时可以直接选用Welch检验。

设第一个营业所随机抽查了\(n_1=28\)个支票账户， 均值为\(\bar x=1025\)(美元)，样本标准差为\(S_1=150\)， 第二个营业所随机抽查了\(n_2=22\)个支票账户， 均值为\(\bar y=910\), 样本标准差为\(S_2=125\)。

从原始数据出发计算检验：

t.test(CheckAcct[[1]], CheckAcct[[2]], var.equal = FALSE, alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  CheckAcct[[1]] and CheckAcct[[2]]
## t = 2.956, df = 47.805, p-value = 0.004828
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   36.77503 193.25224
## sample estimates:
## mean of x mean of y 
## 1025.0000  909.9864

p值为0.005，在0.05水平下显著， 两个营业所的账户平均余额有显著差异。

31.3.1.2.4.3 两种工具软件的比较

考虑开发信息系统的两种工具软件的比较。 新的软件声称比旧的软件能加快进度。 设使用旧软件时平均项目完成时间为\(\mu_1\)， 使用新软件时为\(\mu_2\)。

检验： \[ H_0: \mu_1 \leq \mu_2 \longleftrightarrow H_a: \mu_1 > \mu_2 \] (对立假设是新软件的工期更短) 取0.05水平。

t.test(SoftwareTest[['Current']], SoftwareTest[['New']],
       var.equal = FALSE, alternative = 'greater')
## 
##  Welch Two Sample t-test
## 
## data:  SoftwareTest[["Current"]] and SoftwareTest[["New"]]
## t = 2.2721, df = 21.803, p-value = 0.01665
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  9.514309      Inf
## sample estimates:
## mean of x mean of y 
##       325       286

p值为0.017, 使用新的工具软件的平均项目完成时间显著地少于使用旧工具软件。

31.3.1.2.4.4 两种工艺所需时间的比较

工厂希望比较同一生产问题两种工艺各自需要的时间， 将选用时间较短的工艺。 可以一组工人用工艺1，一组工人用工艺2， 随机分组。 两组之间工人的个体差异会使得差异的比较误差较大。

所以，让同一个工人分别采用两种工艺， 次序先后随机。 这样不会引入个体差异的影响， 精度较高。 使用成对检验。 用双侧检验，水平0.05。

t.test(Matched[["Method 1"]], Matched[["Method 2"]], 
       paired=TRUE, alternative = "two.sided")
## 
##  Paired t-test
## 
## data:  Matched[["Method 1"]] and Matched[["Method 2"]]
## t = 2.1958, df = 5, p-value = 0.07952
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.05120834  0.65120834
## sample estimates:
## mean of the differences 
##                     0.3

31.3.1.3 多元均值置信域

待补充。

31.3.2 比例的假设检验和置信区间

31.3.2.1 单个比例的问题

设要比较的是总体的某个比例， 如选民中对候选人甲的支持率。 这时总体是两点分布总体， 支持为1， 不支持为0，参数为\(p=P(X=1)\)。 其它的比例问题类似， \(p\)是“成功概率”。

关于单总体的比例，也有双侧、右侧、左侧检验问题： \[\begin{aligned} & H_0: p=p_0 \longleftrightarrow H_a: p \neq p_0 \\ & H_0: p \leq p_0 \longleftrightarrow H_a: p > p_0 \\ & H_0: p \geq p_0 \longleftrightarrow H_a: p < p_0 \end{aligned}\]

在R中， 大样本情况下可以用prop.test(x, n, p=p0)作单个比例的检验， n是样本量，x是样本中符合条件（如支持）的个数， 检验使用近似卡方统计量， 用alternative选项选择双侧、右侧、左侧检验。 这个函数也给出比例\(p\)的近似\(1-\alpha\)置信区间， 用conf.level选项指定置信度。

大样本时还有 \[ \frac{\hat p - p}{\sqrt{p(1-p)/n}} \] 近似服从标准正态分布， 可以据此计算Z检验。 自定义的R函数如下：

prop.test.1s <- function(x, n, p=0.5, alternative="two.sided"){
  phat <- x/n
  zstat <- (phat - p)/sqrt(p*(1-p)/n)
  if(alternative=="two.sided"){   # 双侧检验
    pvalue <- 2*(1 - pnorm(abs(zstat)))
  } else if(alternative=="less"){ # 左侧检验
    pvalue <- pnorm(zstat)
  } else if(alternative=="greater"){ # 右侧检验
    pvalue <- 1 - pnorm(zstat)
  } else {
    stop("alternative unknown!")
  }
  
  c(stat=zstat, pvalue=pvalue)
}

因为比例问题本质上是一个二项分布推断问题， 所以R还提供了binom.test(x, n, p=p0)函数进行精确检验并计算精确置信区间， 但其结果略保守。 用法与prop.test()类似。

31.3.2.2 两个比例的比较

设有\(X\)和\(Y\)两个独立的总体， 如男性和女性， 比较\(p_1 = P(X=1)\)和\(p_2 = P(Y=1)\)。

检验问题包括： \[\begin{aligned} H_0: p_1 = p_2 & \longleftrightarrow H_\text{a}: p_1 \neq p_2 \\ H_0: p_1 \leq p_2 & \longleftrightarrow H_\text{a}: p_1 > p_2 \\ H_0: p_1 \geq p_2 & \longleftrightarrow H_\text{a}: p_1 < p_2 \end{aligned}\]

在大样本情况下， 可以用R函数 prop.test(c(nsucc1, nsucc2), c(n1,n2), alternative=...) 进行独立两样本比例的比较检验， 其中n1和n2是两个样本量， nsucc1和nsucc2是两个样本中符合特征的个数（个数除以样本量即样本中的比例）， alternative的选择也是"two.sided"、"less"、"greater"。

还可以比较\(p_1 - p_2\)与某个\(\delta\)。 设两个总体中分别抽取了\(n_1\)和\(n_2\)个样本点， 样本比例分别为\(\hat p_1\)和\(\hat p_2\)。 估计共同的样本比例 \[ \hat p = \frac{n_1 \hat p_1 + n_2 \hat p_2}{n_1 + n_2} \] 当\(\delta=0\)时取统计量 \[ Z = \frac{\hat p_1 - \hat p_2} {\sqrt{\hat p (1 - \hat p) \left( \frac{1}{n_1} + \frac{1}{n_2} \right)}} \] 当\(p_1 = p_2\)时\(Z\)在大样本情况下近似服从标准正态分布。 当\(\delta\neq 0\)时取统计量 \[ Z = \frac{\hat p_1 - \hat p_2 - \delta} {\sqrt{\frac{\hat p_1(1 - \hat p_1)}{n_1} + \frac{\hat p_2(1 - \hat p_2)}{n_2}}} \] 当\(p_1 - p_2 = \delta\)时\(Z\)在大样本情况下近似服从N(0,1)分布。

可以定义Z检验R函数为：

prop.test.2s <- function(x, n, delta=0.0, alternative="two.sided"){
  phat <- sum(x)/sum(n)
  p <- x / n
  if(delta==0.0){
    zstat <- (p[1] - p[2])/sqrt(phat*(1-phat)*(1/n[1] + 1/n[2]))
  } else {
    zstat <- (p[1] - p[2] - delta)/sqrt(p[1]*(1-p[1])/n[1] + p[2]*(1-p[2])/n[2])
  }
  if(alternative=="two.sided"){   # 双侧检验
    pvalue <- 2*(1 - pnorm(abs(zstat)))
  } else if(alternative=="less"){ # 左侧检验
    pvalue <- pnorm(zstat)
  } else if(alternative=="greater"){ # 右侧检验
    pvalue <- 1 - pnorm(zstat)
  } else {
    stop("alternative unknown!")
  }
  
  c(stat=zstat, pvalue=pvalue)
}

独立两总体比例比较的小样本情形意义不大， 可考虑使用Fisher精确检验， R函数为fisher.test()。

31.3.2.3 比例检验的例子

31.3.2.3.1 高尔夫培训女生比例检验例子

橡树溪高尔夫培训机构的女性学员比例较少，只有20%。 为了增加女性学员，设计了促销措施。一个月后调查， 希望证明女性学员比例增加了。

设女性学员比例为\(p\)， \[ H_0: p \leq 0.20 \longleftrightarrow H_a: p > 0.20 \]

取检验水平0.05。 抽查了400名学员，其中100名是女性学员，\(\hat p = 0.25\)。

用prop.test()检验：

prop.test(100, 400, p=0.20, alternative = "greater")
## 
##  1-sample proportions test with continuity correction
## 
## data:  100 out of 400, null probability 0.2
## X-squared = 5.9414, df = 1, p-value = 0.007395
## alternative hypothesis: true p is greater than 0.2
## 95 percent confidence interval:
##  0.2149649 1.0000000
## sample estimates:
##    p 
## 0.25

检验p值为0.0074，小于检验水平0.05， 所以女性学员比例已经显著地高于过去的20%了。

用自定义的大样本Z检验函数prop.test.1s():

prop.test.1s(100, 400, p=0.20, alternative = "greater")
##        stat      pvalue 
## 2.500000000 0.006209665

检验p值为0.0062，与基于卡方检验的prop.test()结果略有差别。

用基于二项分布的binom.test()检验:

binom.test(100, 400, p=0.20, alternative = "greater")
## 
##  Exact binomial test
## 
## data:  100 and 400
## number of successes = 100, number of trials = 400, p-value = 0.008595
## alternative hypothesis: true probability of success is greater than 0.2
## 95 percent confidence interval:
##  0.2146019 1.0000000
## sample estimates:
## probability of success 
##                   0.25

p值为0.0086。

31.3.2.3.2 报税代理分理处的错误率比较

某个代理报税的机构希望比较下属的两个分理处申报纳税返还的出错率。 各自随机选取取了\(n_1=250\)和\(n_2=300\)件申报， 第一个分理处错了35件，错误率\(\hat p_1 = 0.14\); 第二个分理处错了27件，错误率\(\hat p_2 = 0.09\)。

作双侧检验，水平0.10。

用prop.test():

prop.test(c(35,27), c(250,300), alternative = "two.sided")
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(35, 27) out of c(250, 300)
## X-squared = 2.9268, df = 1, p-value = 0.08712
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.007506845  0.107506845
## sample estimates:
## prop 1 prop 2 
##   0.14   0.09

p值0.087，有显著差异，第一处的错误率高。 因为比例需要较大样本量， 所以样本量不太大时， 有人用0.10这样的检验水平。

用自定义的prop.test.2s():

prop.test.2s(c(35,27), c(250,300), alternative = "two.sided")
##       stat     pvalue 
## 1.84618928 0.06486473

31.3.3 方差的假设检验和置信区间

31.3.3.1 单总体方差的假设检验和置信区间

方差和标准差是总体的重要数字特征， 在金融应用中标准差代表波动率， 在工业生产中标准差也是重要的质量指标， 如机床加工精度可以用标准差表示。

对总体方差\(\sigma^2\)可以做双侧、左侧、右侧检验， 与某个\(\sigma_0^2\)比较。

设总体为正态总体，检验问题： \[\begin{aligned} H_0:& \sigma^2 = \sigma_0^2 \longleftrightarrow H_\text{a}: \sigma^2 \neq \sigma_0^2 \\ H_0:& \sigma^2 \geq \sigma_0^2 \longleftrightarrow H_\text{a}: \sigma^2 < \sigma_0^2 \\ H_0:& \sigma^2 \leq \sigma_0^2 \longleftrightarrow H_\text{a}: \sigma^2 > \sigma_0^2 \end{aligned}\]

可以用如下的自定义函数计算单正态总体方差的检验：

var.test.1s <- function(x, n=length(x), var0, alternative="two.sided"){
  if(length(x)==1){ # 输入的是方差
    varx <- x
  } else {
    varx <- var(x)
  }
  xi <- (n-1)*varx/var0
  
  if(alternative=="less"){
    pvalue <- pchisq(xi, n-1)
  } else if (alternative=="right"){
    pvalue <- pchisq(xi, n-1, lower.tail=FALSE)
  } else if(alternative=="two.sided"){
    pvalue <- 2*min(pchisq(xi, n-1), pchisq(xi, n-1, lower.tail=FALSE))
  }
  
  c(statistic=xi, pvalue=pvalue)
}

可以输入样本数据到x中， 输入\(sigma_0^2\)到var0中， 进行检验； 也可以输入\(S^2\)到x中，输入样本量\(n\)到n中， 输入\(sigma_0^2\)到var0中， 进行检验。

31.3.3.2 独立两总体方差的比较

设两个独立的正态总体的方差分别为\(\sigma_1^2\), \(\sigma_2^2\)， 有双侧、左侧、右侧检验： \[\begin{aligned} H_0:& \sigma_1^2 = \sigma_2^2 \longleftrightarrow H_\text{a}: \sigma_1^2 \neq \sigma_2^2 \\ H_0:& \sigma_1^2 \geq \sigma_2^2 \longleftrightarrow H_\text{a}: \sigma_1^2 < \sigma_2^2 \\ H_0:& \sigma_1^2 \leq \sigma_2^2 \longleftrightarrow H_\text{a}: \sigma_1^2 > \sigma_2^2 \end{aligned}\]

设两个总体分别抽取\(n_1\)和\(n_2\)个样本点， 样本方差分别为\(S_x^2\)和\(S_y^2\)。

检验统计量 \[ F = \frac{S_x^2}{S_y^2} \] 在\(\sigma_1^2 = \sigma_2^2\)时\(F\)服从\(F(n_1 - 1, n_2 - 1)\)分布。

如果输入x, y为两个样本的具体值， 可以用var.test(x, y, alternative=...)检验， 其中alternative的选择是"two.sided"、"less"、"greater"。

var.test()假定总体都服从正态分布。 R中mood.test()是一种不要求正态分布假定的检验方法。 bartlett.test()检验多个独立正态总体方差是否全相等。

31.3.3.3 方差检验例子

(待补充)

31.3.4 拟合优度检验

31.3.4.1 各类比例相等的检验

设类别变量\(X\)有\(m\)个不同类别， 抽取\(n\)个样本点后， 各类的频数分别为\(f_1, f_2, \dots, f_m\)。

零假设为所有类的比例相同； 对立假设为各类比例不完全相同。

在零假设下，每个类的期望频数为\(n/m\)。

检验统计量为 \[ \chi^2 = \sum_{i=1}^m \frac{(f_i - \frac{n}{m})^2}{\frac{n}{m}} \] 当\(\chi^2\)超过某一临界值时拒绝零假设。 在\(H_0\)成立且大样本情况下\(\chi^2\)统计量近似服从\(\chi^2(m-1)\)分布 （自由度为组数减1）。

设x中包含各类的频数， R函数chisq.test(x)作各类的总体比例相等的拟合优度卡方检验。

31.3.4.2 各类比例为指定值的检验

设分类变量\(X\)的各类每类有一个总体比例的假设， 希望作为零假设检验。

设共有\(m\)个类，在零假设中， 各类的比例分别为\(p_1, p_2, \dots, p_m\)。

取了\(n\)个样本点组成的样本， 各类的频数分别为\(f_1, f_2, \dots, f_m\)， 期望频数分别为\(np_1, np_2, \dots, np_m\)。

检验统计量为 \[ \chi^2 = \sum_{i=1}^m \frac{(f_i - n p_i)^2}{n p_i} \] 在\(H_0\)下成立且大样本情况下\(\chi^2\)近似服从\(\chi^2(m-1)\)， 自由度为组数减1。 因为是大样本检验法， 每个组的期望频数不能低于5， 否则可以合并较小的类。

设x中包含各类的频数， p中包含\(H_0\)假定的各类的概率， R函数chisq.test(x, p)作各类的总体比例为指定概率的拟合优度卡方检验。

例: 市场占有率的检验例子

某类产品中公司A占30%, 公司B占50%，公司C占20%。 近期C公司推出了新产品代替本公司的老产品。 Scott Marketing Research是一个调查公司， 公司C委托Scott Marketing Research调查分析新产品是否导致了市场占有率的变动。

设\(p_1, p_2, p_3\)分别为三个公司现在的市场占有率，零假设为: \[ H_0: p_1 = 0.30, p_2 = 0.50, p_3 = 0.20 \] 取检验水平0.05。 抽查了200位消费者， 购买三个公司的产品的人数分配为： \[(48, 98, 54)\]

卡方统计量

\[ \chi^2 = \frac{(48 - 200\times 0.3)^2}{200\times 0.3} + \frac{(98 - 200\times 0.5)^2}{200\times 0.5} + \frac{(54 - 200\times 0.2)^2}{200\times 0.2} = 7.34 \]

p值\(=P(\chi^2(3-1) > 7.34) = 0.02548\)， 拒绝零假设， 认为市场占有率有变化。

用chisq.test()计算：

chisq.test(c(48, 98, 54), p=c(0.3, 0.5, 0.2))
## 
##  Chi-squared test for given probabilities
## 
## data:  c(48, 98, 54)
## X-squared = 7.34, df = 2, p-value = 0.02548

31.3.4.3 带有未知参数的单分类变量的拟合优度假设检验

如果\(H_0\)下的概率有未知参数， 先用最大似然估计法估计未知参数， 然后将未知参数代入计算各类的概率， 从而计算期望频数， 计算卡方统计量。 求\(p\)值时自由度要改为“组数减1减去未知参数个数”。

例： 设某次选举有5位候选人， 其中呼声最高的是甲和乙。 问：甲、乙的支持率相等吗？

这不能用独立的两个比例的比较来解决， 因为这两个比例是互斥的， 不属于独立情况。

将类别简化为：甲，乙，其他。 零假设为“甲、乙的支持率都等于\(p\)， 其他三位候选人的支持率为\(1-2p\)，\(p\)未知”； 对立假设是甲、乙的支持率不相等。

假设调查了1000位选民， 其中300名支持甲，200名支持乙，500名支持其他三位候选人。

在假定零假设成立的情况下先用最大似然估计法估计未知参数\(p\)： \(\hat p = (300 + 200)/1000/2 = 0.25\)。 这样，三个类的概率估计为\((0.25, 0.25, 0.50)\)。

然后，计算检验统计量：

chisq.test(c(300, 200, 500), c(0.25, 0.25, 0.50))
## Warning in chisq.test(c(300, 200, 500), c(0.25, 0.25, 0.5)): Chi-squared
## approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  c(300, 200, 500) and c(0.25, 0.25, 0.5)
## X-squared = 3, df = 2, p-value = 0.2231

其中的\(p\)值所用的自由度错误，我们需要自己计算p值：

res <- chisq.test(c(300, 200, 500), c(0.25, 0.25, 0.50))
## Warning in chisq.test(c(300, 200, 500), c(0.25, 0.25, 0.5)): Chi-squared
## approximation may be incorrect
c(statistic=res$statistic, pvalue=pchisq(res$statistic, res$parameter - 1, lower.tail=FALSE))
## statistic.X-squared    pvalue.X-squared 
##          3.00000000          0.08326452

检验\(p\)值为0.08而非原来的0.22。

31.3.5 检验分布类型

vcd包提供了一个goodfit函数， 可以用来拟合指定的某种理论分布(包括泊松、二项、负二项分布）， 并检验服从该理论分布的零假设。

例如，我们生成一组速率参数为2的泊松随机数， 检验其分布是否泊松分布：

library(vcd)
## Warning: 程辑包'vcd'是用R版本3.6.3 来建造的
## 载入需要的程辑包：grid
set.seed(101)
datax <- rpois(100, 2)
summary(goodfit(datax, "poisson"))
## 
##   Goodness-of-fit test for poisson distribution
## 
##                       X^2 df  P(> X^2)
## Likelihood Ratio 4.289456  5 0.5085374

检验其是否服从二项分布，取二项分布试验数为10：

summary(goodfit(datax, "binomial", par = list(size = 10)))
## 
##   Goodness-of-fit test for binomial distribution
## 
##                       X^2 df  P(> X^2)
## Likelihood Ratio 5.052451  5 0.4095126

与二项分布也相容。 在检验分布类型时， 这种两种分布都不能拒绝的情况是常见的。

31.3.6 列联表独立性卡方检验

设分类变量\(X\)有\(m\)个类， 分类变量\(Y\)有\(k\)个类， 抽取了样本量为\(n\)的样本， 设\(X\)的第\(i\)类与\(Y\)的第\(j\)类交叉的频数为\(f_{ij}\)。

零假设为\(X\)与\(Y\)相互独立， 即行变量与列变量相互独立。

交叉频数列成列联表形式，第\(i\)行第\(j\)列为\(f_{ij}\)。

第\(i\)行的行和为\(f_{i\cdot}\), \(X\)的第\(i\)类的百分比为\(r_{i} = f_{i\cdot}/n\);

第\(j\)列的列和为\(f_{\cdot j}\), \(Y\)的第\(j\)类的百分比为\(c_{j} = f_{\cdot j}/n\)。

当\(X\), \(Y\)独立时， \((i,j)\)格子的期望频数为 \(E_{ij} = n \times r_i \times c_j\)。

列联表独立性检验统计量

\[ \chi^2 = \sum_{i=1}^m \sum_{j=1}^k \frac{(f_{ij} - E_{ij})^2}{E_{ij}} \]

其中\(E_{ij} = n r_i c_j\)， \(r_i\)为\(X\)第\(i\)类的百分比， \(c_j\)为\(Y\)第\(j\)类的百分比。

在\(H_0\)下\(\chi^2\)近似服从\((m-1)(k-1)\)自由度的卡方分布。 设统计量值为\(c\)， p值为\(P(\chi^2((m-1)(k-1)) > c)\)。

如果x, y分别是两个变量的原始观测值， chisq.test(x, y)可以做列联表独立性检验。 如果x保存了矩阵格式的列联表， 矩阵行名是\(X\)各个类的名称， 矩阵列名是\(Y\)各个类的名称， 则chisq.test(x)可以做列联表独立性检验。

列联表卡方检验法的检验统计量在零假设下的卡方分布是大样本情况的近似分布。 如果每个变量仅有两个类，每个类的期望频数不能少于5; 如果有多个单元格，期望频数少于5的单元格的个数不能超过20%， 否则应该合并较小的类。

31.3.6.1 列联表独立性卡方检验例子

31.3.6.1.1 性别与啤酒种类的独立性检验

Alber’s Brewery of Tucson, Arizona是啤酒制造与销售商。 有三类啤酒产品：淡啤酒，普通啤酒，黑啤酒。

了解不同顾客喜好有利于制定更精准的销售策略。 希望了解男女顾客对不同类型的偏好有没有显著差异， 实际就是检验性别与啤酒类型偏好的独立性。

抽查了150位顾客，得到如下的列联表：

ctab.beer <- rbind(c(
  20, 40, 20),
  c(30,30,10))
colnames(ctab.beer) <- c("Light", "Regular", "Dark")
rownames(ctab.beer) <- c("Male", "Female")
addmargins(ctab.beer)
##        Light Regular Dark Sum
## Male      20      40   20  80
## Female    30      30   10  70
## Sum       50      70   30 150

列联表独立性检验：

chisq.test(ctab.beer)
## 
##  Pearson's Chi-squared test
## 
## data:  ctab.beer
## X-squared = 6.1224, df = 2, p-value = 0.04683

在0.05水平下认为啤酒类型偏好与性别有关。 男性组的偏好分布、女性组的偏好分布、所有人的偏好分布：

tab2 <- round(prop.table(addmargins(ctab.beer, 1), 1), 3)
rownames(tab2)[3] <- "All"
tab2
##        Light Regular  Dark
## Male   0.250   0.500 0.250
## Female 0.429   0.429 0.143
## All    0.333   0.467 0.200

可以看出男性中对淡啤酒偏好偏少， 女性中对淡啤酒偏好偏多。

31.3.7 非参数检验

待完成。

References

Kabacoff, Robert I. 2012. R语言实战. 人民邮电出版社.

Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with S. 4th Ed. Springer.