Calcolo di un intervallo di confidenza medio senza memorizzazione di tutti i punti dati

Question

Per n di grandi dimensioni (vedi sotto per sapere cosa è abbastanza grande), è sicuro trattare, secondo il teorema del limite centrale, la distribuzione del campione medio come normale (gaussiano) ma io ' d piace una procedura che dia un intervallo di confidenza per qualsiasi n . Il modo per farlo è usare una distribuzione Student T con gradi di libertà n-1 .

Quindi la domanda è, dato un flusso di punti dati che raccogli o incontri uno alla volta, come calcolare un c (ad esempio, c = .95 ) intervallo di confidenza sulla media dei punti dati (senza memorizzare tutti i dati rilevati in precedenza)?

Un altro modo per chiederlo è: come tenere traccia del primo e del secondo momento di un flusso di dati senza archiviare l'intero flusso?

DOMANDA BONUS: puoi tenere traccia dei momenti più elevati senza memorizzare l'intero flusso?

Answer 1

[Grazie enormi a John D Cook per gran parte di ciò che ho imparato nel mettere insieme questa risposta!]

Innanzitutto, ecco il motivo per non utilizzare la somma dei quadrati: http: //www.johndcook.com/blog/2008/09/26/

Cosa invece dovresti fare:

Tieni traccia del conteggio (n), della media (u) e delle quantità da cui è possibile determinare la varianza del campione e l'errore standard. (Adattato da http://www.johndcook.com/standard_deviation.html .)

Inizializza n = u = s = 0 .

Per ogni nuovo punto dati, x :

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

La varianza del campione è quindi s / (n-1) , la varianza della media del campione è s / (n-1) / n e lo standard l'errore della media del campione è SE = sqrt (s / (n-1) / n) .

Resta da calcolare l'intervallo di confidenza Student-t c ( c in (0,1)):

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

dove g è il cdf inverso (aka quantile) della distribuzione Student-t con media 0 e varianza 1, prendendo una probabilità e i gradi di libertà (uno in meno del numero di punti dati ):

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

dove irib è la funzione beta incompleta regolarizzata inversa:

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

dove rib è la funzione beta incompleta regolarizzata:

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

dove B (a, b) è la funzione beta di Eulero e B (x0, x1, a, b) è la funzione beta incompleta:

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

Tipiche librerie numeriche / statistiche avranno implementazioni della funzione beta (o direttamente il cdf inverso della distribuzione Student-t). Per C, lo standard di fatto è la Gnu Scientific Library (GSL). Spesso viene fornita una versione a 3 argomenti della funzione beta; la generalizzazione in 4 argomenti è la seguente:

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)

Infine, ecco un'implementazione in 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}]

Confronta con

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

oppure, quando è disponibile l'intero elenco di punti dati,

MeanCI[list, ConfidenceLevel->c]

Infine, se non vuoi caricare librerie matematiche per cose come la funzione beta, puoi hardcodificare una tabella di ricerca per -g ((1-c) / 2, n-1) . Qui è per c = .95 e 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.98

Answer 2

Ecco un articolo su come calcolare la media e la deviazione standard in un unico passaggio , non memorizzare alcun dato. Dopo aver ottenuto queste due statistiche, è possibile stimare un intervallo di confidenza. Un intervallo di confidenza del 95% sarebbe [media - 1,96 * stdev, media + 1,96 * stdev], ipotizzando una distribuzione normale per i tuoi dati e un gran numero di punti dati.

Per un numero inferiore di punti dati, l'intervallo di confidenza sarebbe [mean - c (n) * stdev, mean + c (n) * stdev] dove c (n) dipende dalle dimensioni del campione e dal livello di confidenza. Per un livello di confidenza del 95%, ecco i tuoi valori di c (n) per 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.0 2.0924 2.085 2.0 , 2.059539, 2.055529, 2.051831, 2.048407, 2.045230

Questi numeri sono g (0,025, n-1) dove g è il CDF inverso della distribuzione t con n gradi di libertà. Se si desidera un intervallo di confidenza del 99%, sostituire 0,025 con 0,005. In generale, per un livello di confidenza di 1-alpha, usa alpha / 2.

Ecco il R che ha generato le costanti sopra.

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

Ecco un post di blog che spiega perché i numeri sopra non sono vicini a 1,96 come ci si potrebbe aspettare.

Answer 3

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

Dove t (x, n-1) è la distribuzione t con n-1 gradi di libertà. se stai utilizzando gsl

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

Non è necessario memorizzare alcun dato oltre la somma dei quadrati. Ora, potresti avere un problema numerico perché la somma dei quadrati può essere piuttosto grande. Puoi usare la definizione alternativa di s

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

e fai due passaggi. Vorrei incoraggiarvi a prendere in considerazione l'uso di gnu scientific library o un pacchetto come R per aiutarti ad evitare problemi numerici. Inoltre, fai attenzione al tuo uso del teorema del limite centrale. L'abuso di questo è in parte responsabile di tutta la crisi finanziaria in atto proprio ora.

Answer 4

Non vuoi accumulare la somma dei quadrati. Le statistiche risultanti sono numericamente inaccurate - finirai per sottrarre due grandi numeri simili. Vuoi mantenere la varianza, o (n-1) * varianza, o qualcosa del genere.

Il modo più semplice è quello di accumulare i punti dati in modo incrementale. La formula non è complicata o difficile da ottenere (consultare il link di John D. Cook).

Un modo ancora più accurato per farlo è combinare i punti dati in modo ricorsivo a coppie. Puoi farlo con la memoria logaritmica in n: register k contiene statistiche per 2 ^ k punti dati più vecchi, che sono combinati con statistiche per 2 ^ k punti più recenti per ottenere statistiche per 2 ^ (k + 1) punti ...

Answer 5

Penso che non devi preoccuparti così tanto della dimensione di n perché presto supererà il numero di 30, dove la distribuzione può essere considerata normale. Usare la ricorsione bayesiana per fare inferenza posteriore sulla media della popolazione e i parametri di varianza, assumendo un modello normale, è il modo migliore, se non si desidera memorizzare punti di dati di campioni precedenti. Puoi dare un'occhiata a questo documento per un'inferenza congiunta per la media e la varianza, e in particolare le equazioni 38a, 38b e 38c.

Answer 6

Penso che tu possa. Dovrei Google / Wikipidia per questo, quindi lo lascerò come esercizio per il lettore.