Вычисления среднего доверительного интервала, не храня все точки данных
Это потребует немного борьбы. Взгляните сюда:
http://jnicookbook.owsiak.org/recipe-No-020/
также, посмотрите здесь образцы, связанные с прохождением Карта в нативный код
https://github.com/mkowsiak/jnicookbook/tree/master/recipeNo037 https://github.com/mkowsiak/jnicookbook/tree/ ведущий / recipeNo038
6 ответов
[Огромный благодаря John D Повар для большой части того, что я изучил в соединении этого ответа!]
Во-первых, вот причина не использовать сумму квадратов: 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 в (0,1)):
u [plus or minus] SE*g((1-c)/2, n-1)
где g инверсия cdf (иначе квантиль) Студенческого-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) Euler бета-функция и 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 Студенческого-t распределения непосредственно). Для C фактический стандарт является Гну научной библиотекой (GSL). Часто версия с 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]
Наконец, если Вы не хотите загружать математические библиотеки для вещей как бета-функция, Вы можете hardcode справочная таблица для -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 точек данных.
Как показывает опыт, необходимо использовать Студента-t вместо нормального приближения для n по крайней мере до 300, не 30 на расхожее мнение. Cf. http://www.johndcook.com/blog/2008/11/12/.
См. также "Улучшение Сжатый подсчет" Ping Li Cornell.
Вот статья о том, как вычислить среднее и стандартное отклонение в единственной передаче, не храня данных. После того как у Вас есть эти два статистических данных, можно оценить доверительный интервал. 95%-й доверительный интервал был бы [средний - 1.96*stdev, средний + 1.96*stdev], приняв нормальное распределение для Ваших данных и большого количества точек данных.
Для меньшего числа точек данных Ваш доверительный интервал был бы [средний - c (n) *stdev, средний + c (n) *stdev], где c (n) зависит от Вашего объема выборки и Вашего доверительного уровня. Для 95%-го доверительного уровня вот Ваши значения c (n) для 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), где g является обратный CDF t распределения с n степенями свободы. Если Вы хотели 99%-й доверительный интервал, замените 0.025 0,005. В целом, для доверительного уровня 1 альфы, используйте альфу/2.
Вот команда R, которая генерировала константы выше.
n = seq(2, 30); qt(0.025, n-1)
Вот сообщение в блоге, объясняющее, почему числа выше не как близко к 1,96, как Вы могли бы ожидать.
sigma = sqrt( (q - (s*s/n)) / (n-1) )
delta = t(1-c/2,n-1) * sigma / sqrt(n)
Где t (x, n-1) является t-распределением с n-1 степенями свободы. если Вы используете gsl
t = gsl_cdf_tdist_Qinv (c/2.0, n-1)
Нет никакой потребности хранить любые данные вне суммы квадратов. Теперь, у Вас могла бы быть числовая проблема, потому что сумма квадратов может быть довольно большой. Вы могли использовать альтернативное определение s
sigma = sqrt(sum(( x_i - s/n )^2 / (n-1)))
и сделайте две передачи. Я поощрил бы Вас рассматривать использование гну научная библиотека или пакет как R, чтобы помочь Вам избежать числовых проблем. Кроме того, быть осторожным относительно Вашего использования центральной предельной теоремы. Злоупотребление им частично виновато в целом финансовом кризисе, продолжающемся прямо сейчас.
Вы не хотите накапливать сумму квадратов. Получающиеся статистические данные численно неточны - Вы закончите тем, что вычли два больших, подобных числа. Вы хотите поддержать различие или (n-1) *variance, или что-то как этот.
Простой путь состоит в том, чтобы накопить точки данных инкрементно. Формула не является сложной или трудно произойти (см. ссылку John D. Cook).
Еще более точный способ сделать это состоит в том, чтобы объединить точки данных попарно рекурсивно. Можно сделать это с памятью, логарифмической в n: зарегистрируйтесь k содержит статистику для 2^k более старые точки данных, которые объединены со статистикой для 2^k более новые точки для получения статистики для 2^ (k+1) точки...
Я думаю, что Вы можете. Я имел бы к Google/Wikipidia для него так, я оставлю это как осуществление для читателя.
Я думаю, что не стоит так беспокоиться о размере n, потому что он скоро превысит число 30, где распределение можно считать нормальным. Использование Байесовской рекурсии, чтобы сделать апостериорный вывод о среднем заселении и параметрах дисперсии, предполагая нормальную модель, я думаю, лучший способ, если вы не хотите хранить никаких точек данных из предыдущих выборок.
Вы можете взглянуть на этот документ для совместного вывода о среднем значении и дисперсии, а именно на уравнения 38a, 38b и 38c.