数据科学18 | 统计推断-渐近性
渐近性（asymptopia）是样本量接近于无穷大时统计行为的一个术语。渐近统计即大样本统计主要研究当样本量n→∞时统计方法的有关渐进性质。渐近性有助于简单的统计推断和估计，也是频率解释概率的基础。
1. 大数定律
大数定律(Law of Large Numbers)：随着样本量的增加，样本均值收敛于总体均值。
随机变量服从正态分布
n <- 10000
means <- cumsum(rnorm(n))/(1:n)
#生成10000个标准正态分布随机数，求累积平均值
#即第1个观测值的平均值、前2个观测值的平均值、前3个观测值的平均值，以此类推
library(ggplot2)
g <- ggplot(data.frame(x = 1:n, y = means), aes(x = x, y = y))
g <- g + geom_hline(yintercept = 0) + geom_line(size = 2)
g <- g + labs(x = "Number of obs", y = "Cumulative mean")
g
样本量增大，样本均值越接近总体均值0。
随机变量服从伯努利分布
means <- cumsum(sample(0:1, n, replace = TRUE))/(1:n)
#生成10000个服从伯努利分布的随机样本，求累积平均值
g <- ggplot(data.frame(x = 1:n, y = means), aes(x = x, y = y))
g <- g + geom_hline(yintercept = 0.5) + geom_line(size = 2)
g <- g + labs(x = "Number of obs", y = "Cumulative mean")
g
样本量增大，样本均值越接近总体均值0.5。
- 当一个估计量收敛于想要估计的总体参数时，这个估计量满足一致性（相合性）。
- 大数定律表明，IID样本的样本均值与总体均值是一致的，样本方差和样本标准差也满足一致性。
2. 中心极限定理
中心极限定理(Central Limit Theorem)：随着样本量的增加，IID样本的样本均值的分布收敛于正态分布。
根据中心极限定理，当n→∞时，随机变量X的样本均值 的分布近似正态分布 ，随机变量 ，近似标准正态分布 。
例：假设 为第 次抛不规则硬币的结果，取值为0或1，取值为1的概率为 。
样本均值为 ；总体均值 ；总体方差 ；均值的标准误为 ；则n→∞时，变量 近似标准正态分布。
假设硬币是规则的，p=0.5，Y的分布：
随着n增大，样本均值的分布近似于正态分布，Y的分布近似于标准正态分布。
假设硬币不规则，p=0.9，Y的分布：
抛30次硬币时，Y的分布不是很好的近似正态分布，但是当n增大时，分布将近似正态分布。
CLT应用：估计量的置信区间
置信区间估计用一个区间来估计参数值。如果多次抽取样本量为n的样本集，每次计算1个估计量的置信区间，其中95%的置信区间包含总体参数，则对于一个样本集中计算的95%置信区间，有95%的信心认为该区间包含总体参数。
根据中心极限定理，样本均值 近似正态分布，均值为?，标准差为 。
样本均值 在区间(?-2?/√?, ?+2?/√?)内的概率约为95%， ±2?/√?为均值?的95%的置信区间，标准正态分布第97.5百分位数约为1.96，接近于2。标准正态分布的第95百分位数约为1.645， ±1.6452?/√?为均值?的90%的置信区间。
例：UsingR包的father.son数据集
library(UsingR)
data(father.son)
x <- father.son$sheight
(mean(x) + c(-1, 1) * qnorm(0.975) * sd(x)/sqrt(length(x)))/12
[1] 5.710 5.738儿子身高均值的95%置信区间为5.710到5.738。
二项分布的参数置信区间
若 为第 次抛不规则硬币的结果，取值为0或1，取值为1的概率为 ， ，样本均值为 。
p的置信区间为 ，这个置信区间称为Wald置信区间。
p的95%的置信区间可以用 ，快速计算。
例：假设竞选中，随机抽样的100名选民有56人打算投你一票，能否保证获得超过50%的选票赢得竞选？即 ，计算赢得竞选概率p的置信区间。
快速计算： ，p的95%的置信区间为0.56±0.1，即(0.46,0.66)，不能排除低于50%的可能性。
一般来说，二项分布试验中，小数点后1位的变化需要样本量为100，2位需要10 000，3位需要1000 000。
round(1/sqrt(10^(1:6)), 3)
[1] 0.316 0.100 0.032 0.010 0.003 0.001计算Wald置信区间：
0.56 + c(-1, 1) * qnorm(0.975) * sqrt(0.56 * 0.44/100)
[1] 0.4627 0.6573binom.test( )计算置信区间：
binom.test(56, 100)$conf.int
[1] 0.4572 0.6592
attr(,"conf.level")#默认为95%
[1] 0.95若 为第 次抛不规则硬币的结果，取值为0或1，取值为1的概率为 ， ，样本均值为 。
n<-20
pvals <- seq(0.1, 0.9, by = 0.05)#设置参数p以0.05的变化从0.1增大到0.9
nosim <- 1000
coverage <- sapply(pvals, function(p) {
phats <- rbinom(nosim, prob = p, size = n)/n #计算1000次模拟的样本均值
ll <- phats - qnorm(0.975) * sqrt(phats * (1 - phats)/n) #计算置信区间的下限
ul <- phats + qnorm(0.975) * sqrt(phats * (1 - phats)/n) #置信区间的上限
mean(ll<p&ul>p) #计算置信区间覆盖真实p值的比例
})对于每一个p值，进行1000次模拟，每次模拟抛20次硬币，计算每次模拟得到的样本均值 以及相应的95%的置信区间，再求出1000次模拟中置信区间覆盖真实p值的次数占的比例。
#画出估计的p值的95%置信区间覆盖真实p值的比例
g <- ggplot(data.frame(x = pvals, y = coverage), aes(x = x, y = y))
g <- g + scale_y_continuous(limits = c(0.75, 1.00))
g <- g + geom_hline(yintercept = 0.95) + geom_line(size = 2)
g <- g + labs(x = "pvals", y = "coverage")
g
p=0.5时， 得到的置信区间覆盖p值的比例比95%要高；但是大部分情况下，没有得到接近95%的覆盖率。由于n不够大，根据中心极限定理计算置信区间的公式不适用。
一种快速解决的方法： ，取值为1的次数X加上2，取值为0的次数也加上2，得到的置信区间称为Agresti-Coull置信区间。
#Wald置信区间覆盖p值的比例
n<-100
pvals <- seq(0.1, 0.9, by = 0.05)
nosim <- 1000
coverage <- sapply(pvals, function(p) {
phats <- rbinom(nosim, prob = p, size = n)/n
ll <- phats - qnorm(0.975) * sqrt(phats * (1 - phats)/n)
ul <- phats + qnorm(0.975) * sqrt(phats * (1 - phats)/n)
mean(ll<p&ul>p)
})
n=100时，得到的Wald置信区间覆盖p值的比例接近95%。
#Agresti-Coull置信区间覆盖p值的比例
n<-20
pvals <- seq(0.1, 0.9, by = 0.05)
nosim <- 1000
coverage <- sapply(pvals, function(p) {
phats <- (rbinom(nosim, prob = p, size = n) + 2)/(n + 4)
ll <- phats - qnorm(0.975) * sqrt(phats * (1 - phats)/n)
ul <- phats + qnorm(0.975) * sqrt(phats * (1 - phats)/n)
mean(ll<p&ul>p)
})
Agresti-Coull置信区间覆盖真实p值的比例往往会高于95%，但是覆盖率过高有时可能由于区间过宽，过于保守。尽管如此，考虑本例建议使用Agresti-Coull置信区间代替Wald置信区间。
泊松分布的参数置信区间
例：一台核泵在94.32天内发生5次故障，给出每天故障率的95%置信区间。
假设发生故障的次数服从X ～ Poisson(?t)，?是故障率，t为天数。对?的估计是 ， ， 为方差估计。
x<-5
t <- 94.32
lambda <- x/t
round(lambda + c(-1, 1) * qnorm(0.975) * sqrt(lambda/t), 3)
[1] 0.007 0.099poisson.test( )可计算95%的置信区间，但区间范围较大，较保守。
poisson.test(x, T = 94.32)$conf
[1] 0.01721 0.12371
attr(,"conf.level")
[1] 0.95重复上文所做的模拟：
lambdavals <- seq(0.005, 0.1, by = 0.01)
nosim <- 1000
t<-100
coverage <- sapply(lambdavals, function(lambda) {
lhats <- rpois(nosim, lambda = lambda * t)/t
ll <- lhats - qnorm(0.975) * sqrt(lhats/t)
ul <- lhats + qnorm(0.975) * sqrt(lhats/t)
mean(ll < lambda & ul > lambda)
})
#画出95%置信区间覆盖真实?值的比例
g <- ggplot(data.frame(x = lambdavals, y = coverage), aes(x = x, y=y))
g <- g + scale_y_continuous(limits = c(0.00, 1.00))
g <- g + geom_hline(yintercept = 0.95) + geom_line(size = 2)
g <- g + labs(x = "lambdavals", y = "coverage")
g 
随着?值的增大，覆盖率接近95%，而对于较小的?值，不应该使用这个区间。
lambdavals <- seq(0.005, 0.1, by = 0.01)
nosim <- 1000
t<-1000
coverage <- sapply(lambdavals, function(lambda) {
lhats <- rpois(nosim, lambda = lambda * t)/t
ll <- lhats - qnorm(0.975) * sqrt(lhats/t)
ul <- lhats + qnorm(0.975) * sqrt(lhats/t)
mean(ll < lambda & ul > lambda)
})
随着监测时间t的延长，覆盖率将收敛于95%。
编辑：李雪纯 冯文清
校审：张健 罗鹏