math - 在不存储所有数据点的情况下计算平均置信区间

Question

对于大n（请参阅下文了解如何确定足够大），通过中心极限定理将样本均值的分布视为正态（高斯）是安全的，但我想要一个为任何n. 这样做的方法是使用具有n-1自由度的学生 T 分布。

所以问题是，给定您一次收集或遇到一个数据点的流，您如何计算c（例如，c=.95）数据点均值的置信区间（不存储所有先前遇到的数据）？

另一种问这个问题的方法是：如何在不存储整个流的情况下跟踪数据流的第一时刻和第二时刻？

奖励问题：您可以在不存储整个流的情况下跟踪更高的时刻吗？

Answer 1

这是一篇关于如何在单次通过中计算均值和标准差的文章，而不存储任何数据。一旦你有了这两个统计数据，你就可以估计一个置信区间。[mean - 1.96*stdev, mean + 1.96*stdev]假设您的数据和大量数据点呈正态分布，则 95% 的置信区间为 。

对于较少数量的数据点，您的置信区间将取决于您的样本量和置信水平[mean - c(n)*stdev, mean + c(n)*stdev]。c(n)对于 95% 的置信水平，这里是您的值c(n)for n= 2, 3, 4, ..., 30

12.70620，4.302653，3.182446，2.776445，2.570582，2.446912，2.364624，2.306004，2.262157，2.228139，2.200985，2.178813，2.160369，2.144787，2.131450，2.119905，2.109816，2.100922，2.093024，2.085963，2.079614，2.073873，2.068658，2.063899，2.059539， 2.055529、2.051831、2.048407、2.045230

这些数字是具有g(0.025, n-1)ng个自由度的 t 分布的逆 CDF。如果您想要 99% 的置信区间，请将 0.025 替换为 0.005。通常，对于 的置信水平1-alpha，使用alpha/2。

这是生成上述常量的R命令。

n = seq(2, 30); qt(0.025, n-1)

这是一篇博客文章，解释了为什么上面的数字不像您预期​​的那样接近 1.96。

Answer 2

[非常感谢约翰·D·库克（John D Cook）在整理这个答案时学到了很多东西！]

首先，这是不使用平方和的原因：http: //www.johndcook.com/blog/2008/09/26/

你应该怎么做：

跟踪计数 (n)、平均值 (u) 和数量 (s)，从中可以确定样本方差和标准误差。（改编自http://www.johndcook.com/standard_deviation.html。）

初始化n = u = s = 0.

对于每个新数据点，x：

u0 = u;
n ++;
u += (x - u) / n;
s += (x - u0) * (x - u);

则样本方差为s/(n-1)，样本均值的方差为 ，样本均值s/(n-1)/n的标准误为SE = sqrt(s/(n-1)/n)。

仍然需要计算 Student-tc置信区间（c在 (0,1) 中）：

u [plus or minus] SE*g((1-c)/2, n-1)

其中g是均值为 0 方差为 1 的 Student-t 分布的逆 cdf（又名分位数），取概率和自由度（比数据点数少一）：

g(p,df) = sign(2*p-1)*sqrt(df)*sqrt(1/irib(1, -abs(2*p-1), df/2, 1/2) - 1)

其中irib是逆正则化不完全 beta 函数：

irib(s0,s1,a,b) = z such that rib(s0,z,a,b) = s1

哪里rib是正则化不完全 beta 函数：

rib(x0,x1,a,b) = B(x0,x1,a,b) / B(a,b)

其中B(a,b)是 Euler beta 函数，B(x0,x1,a,b)是不完全 beta 函数：

B(a,b) = Gamma(a)*Gamma(b)/Gamma(a+b) = integral_0^1 t^(a-1)*(1-t)^(b-1) dt
B(x0,x1,a,b) = integral_x0^x1 t^(a-1)*(1-t)^(b-1) dt

典型的数值/统计库将具有 beta 函数的实现（或直接的 Student-t 分布的逆 cdf）。对于 C，事实上的标准是 Gnu Scientific Library (GSL)。通常给出 beta 函数的 3 参数版本；对 4 个论点的概括如下：

B(x0,x1,a,b) = B(x1,a,b) - B(x0,a,b)
rib(x0,x1,a,b) = rib(x1,a,b) - rib(x0,a,b)

最后，这是 Mathematica 中的一个实现：

(* Take current {n,u,s} and new data point; return new {n,u,s}. *)
update[{n_,u_,s_}, x_] := {n+1, u+(x-u)/(n+1), s+(x-u)(x-(u+(x-u)/(n+1)))}

Needs["HypothesisTesting`"];
g[p_, df_] := InverseCDF[StudentTDistribution[df], p]

(* Mean CI given n,u,s and confidence level c. *)
mci[n_,u_,s_, c_:.95] := With[{d = Sqrt[s/(n-1)/n]*g[(1-c)/2, n-1]}, 
  {u+d, u-d}]

相比于

StudentTCI[u, SE, n-1, ConfidenceLevel->c]

或者，当整个数据点列表可用时，

MeanCI[list, ConfidenceLevel->c]

最后，如果您不想为 beta 函数之类的东西加载数学库，您可以为-g((1-c)/2, n-1). 这里是c=.95和n=2..100：

12.706204736174698，4.302652729749464，3.182446305283708，2.7764451051977934，2.570581835636314，2.4469118511449666，2.3646242515927853，2.306004135204168，2.262157162798205，2.2281388519862735，2.2009851600916384，2.178812829667226，2.1603686564627917，2.1447866879178012，2.131449545559774，2.1199052992212533，2.1098155778333156，2.100922040241039，2.093024054408307，2.0859634472658626，2.0796138447276835，2.073873067904019，2.0686576104190477，2.0638985616280254，2.0595385527532963， 2.05552943864287，2.051830516480281，2.048407141795243，2.0452296421327034，2.042272456301236，2.039513446396408，2.0369333434600976，2.0345152974493392，2.032244509317719，2.030107928250338，2.0280940009804462，2.0261924630291066，2.024394163911966，2.022690920036762，2.0210753903062715，2.0195409704413745，2.018081702818439，2.016692199227822，2.0153675744437627，2.0141033888808457，2.0128955989194246，2.011740513729764，2.0106347576242314，2.0095752371292335，2.0085591121007527，2.007583770315835，2.0066468050616857，2.005745995317864，2.0048792881880577，2.004044783289136，2.0032407188478696，2.002465459291016，2.001717484145232，2.000995378088259，2.0002978220142578，1.9996235849949402，1.998971517033376，1.9983405425207483，1.997729654317692，1.9971379083920013，1.9965644189523084，1.996008354025304，1.9954689314298386，1.994945415107228， 1.9944371117711894，1.9939433678456229，1.993463566661884，1.9929971258898527，1.9925434951809258，1.992102154002232，1.9916726096446793，1.9912543953883763，1.9908470688116922，1.9904502102301198，1.990063421254452，1.989686323456895，1.9893185571365664，1.9889597801751728，1.9886096669757192，1.9882679074772156，1.9879342062390228，1.9876082815890748，1。9872898648311672，1.9869786995062702，1.986674540703777，1.986377154418625，1.9860863169510985，1.9858018143458114，1.9855234418666061，1.9852510035054973，1.9849843115224508，1.9847231860139618，1.98446745450849，1.9842169515863888

渐近地逼近 的正态 (0,1) 分布的逆 CDF c=.95，即：

-sqrt(2)*InverseErf(-c) = 1.959963984540054235524594430520551527955550...

有关反函数，请参见http://mathworld.wolfram.com/InverseErf.html 。erf()请注意，g((1-.95)/2,n-1)在至少有 474 个数据点之前，它不会四舍五入到 1.96。当有 29 个数据点时，它四舍五入为 2.0。

根据经验，您应该使用 Student-t 而不是正常的近似值，n最多至少 300，而不是按照传统观点的 30。参照。http://www.johndcook.com/blog/2008/11/12/。

另见康奈尔大学李平的“改进压缩计数”。

Answer 3

   sigma = sqrt( (q - (s*s/n)) / (n-1) )
   delta = t(1-c/2,n-1) * sigma / sqrt(n)

其中 t(x, n-1) 是具有 n-1 个自由度的 t 分布。如果您使用的是gsl

t = gsl_cdf_tdist_Qinv (c/2.0, n-1)

除了平方和之外，无需存储任何数据。现在，您可能会遇到数值问题，因为平方和可能非常大。您可以使用 s 的替代定义

sigma = sqrt(sum(( x_i - s/n )^2 / (n-1)))

并通过两次。我鼓励您考虑使用gnu 科学库或R之类的包来帮助您避免数值问题。此外，请注意使用中心极限定理。滥用它部分归咎于目前正在发生的整个金融危机。

Answer 4

您不想累积平方和。得到的统计数据在数字上是不准确的——你最终会减去两个相似的大数字。你想保持方差，或 (n-1)*variance，或类似的东西。

直接的方法是增量地累积数据点。该公式并不复杂或难以推导（参见 John D. Cook 的链接）。

一种更准确的方法是成对递归地组合数据点。您可以使用 n 中的内存对数来执行此操作：寄存器 k 保存 2^k 个旧数据点的统计信息，这些数据与 2^k 个较新点的统计信息相结合以获得 2^(k+1) 个点的统计信息...

Answer 5

我认为你不必太担心 的大小，n因为它很快就会超过 30 的数量，这里的分布可以认为是正常的。如果您不想存储来自先前样本的任何数据点，我认为使用贝叶斯递归对总体均值和方差参数进行后验推断（假设是正常模型）是最好的方法。您可以查看此文档以了解均值和方差的联合推断，特别是方程 38a、38b 和 38c。

Answer 6

我想你可以。我必须在谷歌/维基百科上找到它，所以我会把它作为练习留给读者。