专栏首页优雅R数据科学18 | 统计推断-渐近性

数据科学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.6573

binom.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.099

poisson.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%。

编辑:李雪纯 冯文清

校审:张健 罗鹏

本文分享自微信公众号 - 优雅R(elegant-r)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-06-10

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 数据科学19 | 统计推断-t分布置信区间

    当样本量足够大,总体标准差已知时,根据中心极限定理可以用标准正态分布估计总体均值;t分布适用于小样本估计呈正态分布的总体均值。

    王诗翔呀
  • 「R」数据可视化13 : 相关性图

    在生物信息领域我们常常使用R语言对数据可视化。在对数据可视化的时候,我们需要明确想要展示的信息,从而选择最为合适的图突出该信息。本系列文章将介绍多种基于不同R包...

    王诗翔呀
  • 「R」Shiny:工作流(二)调试

    当你开始编写应用程序时,几乎可以确定会出错。导致大多数错误的原因是我们心里的 Shiny 设计模型与 Shiny 实际的运行情况的不匹配。当你阅读本文时,你的思...

    王诗翔呀
  • c++中两个类互相引用的问题

                    删除指向不完整“Q2DTorusNode”类型的指针;没有调用析构函数                 1> ...

  • c++中两个类互相引用的问题

  • c++中两个类互相引用的问题

                    删除指向不完整“Q2DTorusNode”类型的指针;没有调用析构函数                 1> ...

  • 一些关于python的小感想

    python是一门优秀的语言,但随之而来的是大量的知识,各种模块,相信一个人的大脑是很难记住如此多的内容。这时后的我们就应该想办法避免去记忆这么多的内容。 1。...

    用户1149564
  • Python中的collections模块

    接下来主要对collections模块中的常用集合类进行介绍,调用collections模块:

    oYabea
  • python笔记47-面试题:如何判断字典a在字典b

    已知一个dict 比如a = {“a”:1},另一个dict比如为b = {“a”:1,”b”:2},如何判断a是否在与b中。 一般在接口测试的时候,返回的参数...

    上海-悠悠
  • 升级为私有全栈云的腾讯云TStack究竟强在了哪儿?

    在刚刚结束的美国丹佛Open Infrastructure峰会上,腾讯云对全新升级为私有全栈云的TStack进行了详细的介绍,引起峰会现场众多OpenStack...

    腾讯技术工程官方号

扫码关注云+社区

领取腾讯云代金券