[Огромное спасибо John D Cook за многое из того, что я узнал, составив этот ответ!]
Во-первых, вот причина не использовать сумму квадратов: http://www.johndcook.com/blog/2008/09/26/
Что вы должны сделать вместо этого:
Следите за количеством (n), средним (u) и количеством (величинами), по которым можно определить дисперсию выборки и стандартную ошибку. (Адаптировано с 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).
Осталось вычислить доверительный интервал Стьюдента-t c (c in (0,1)):
u [plus or minus] SE*g((1-c)/2, n-1)
, где g - обратный cdf (он же квантиль) распределения Student-t со средним 0 и дисперсией 1, принимая вероятность и степени свободы (на единицу меньше, чем число точек данных):
g(p,df) = sign(2*p-1)*sqrt(df)*sqrt(1/irib(1, -abs(2*p-1), df/2, 1/2) - 1)
где irib - обратная регуляризованная неполная бета-функция:
irib(s0,s1,a,b) = z such that rib(s0,z,a,b) = s1
, где rib - регуляризованная неполная бета-функция:
rib(x0,x1,a,b) = B(x0,x1,a,b) / B(a,b)
, где B(a,b) - бета-функция Эйлера, а B(x0,x1,a,b) - неполная бета-функция:
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
Типичные числовые / статистические библиотеки будут иметь реализации бета-функции (или обратного cdf дистрибутива Student-t напрямую). Для C стандартом де-факто является Научная библиотека Gnu (GSL). Часто дается бета-функция с тремя аргументами; обобщение на 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]
Наконец, если вы не хотите загружать математические библиотеки для таких вещей, как бета-функция, вы можете жестко закодировать таблицу поиска для -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
, который асимптотически приближается к обратному CDF нормального (0,1) распределения для c=.95, то есть:
-sqrt(2)*InverseErf(-c) = 1.959963984540054235524594430520551527955550...
См. http://mathworld.wolfram.com/InverseErf.html для обратной функции erf().
Обратите внимание, что g((1-.95)/2,n-1) не округляется до 1,96, пока не будет хотя бы 474 точек данных. Он округляется до 2,0, когда имеется 29 точек данных.
Как правило, вы должны использовать Student-t вместо нормального приближения для n до не менее 300, а не 30 для обычной мудрости. Ср http://www.johndcook.com/blog/2008/11/12/.
См. Также «Улучшение счета сжатых» Пинг Ли из Корнелла.