6

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

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

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

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

4

6 回答 6

5

这是一篇关于如何在单次通过中计算均值和标准差的文章,而不存储任何数据。一旦你有了这两个统计数据,你就可以估计一个置信区间。[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。

于 2008-11-12T00:37:26.870 回答
4

[非常感谢约翰·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=.95n=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/

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

于 2008-11-12T06:25:51.610 回答
2
   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之类的包来帮助您避免数值问题。此外,请注意使用中心极限定理。滥用它部分归咎于目前正在发生的整个金融危机。

于 2008-11-12T00:40:02.120 回答
1

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

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

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

于 2008-11-12T00:46:04.090 回答
1

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

于 2010-02-05T18:00:18.993 回答
-7

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

于 2008-11-12T00:37:22.877 回答