Test sequenziali con applicazioni al moto browniano

UNIVERSITÀ DEGLI STUDI DI TORINO DIPARTIMENTO DI MATEMATICA GIUSEPPE PEANO SCUOLA DI SCIENZE DELLA NATURA Corso di Laurea in Matematica per la Finanza e l’Assicurazione Tesi di Laurea Triennale Test sequenziali con applicazioni al moto browniano Relatore: Cristina Zucca Candidato: Fabrizio Lasaponara ANNO ACCADEMICO 2014/2015 Indice 1 Teoria classica sui test di ipotesi: concetti di base 3 2 Test di ipotesi sequenziale 2.1 Nozione iniziale di test sequenziale . . . . . . . . . . . . . . . . . 2.2 Conseguenze della scelta di ogni particolare test . . . . . . . . . . 2.2.1 La funzione Operating Characteristic (OC) . . . . . . . . 2.2.2 La funzione Average Sample Number (ASN) . . . . . . . 2.3 Principi di selezione di un sequential test . . . . . . . . . . . . . . 2.3.1 Livello di preferenza per accettazione e rifiuto dell’ipotesi nulla . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Condizioni imposte dall’OC function . . . . . . . . . . . 2.3.3 ASN come base per la scelta di un test sequenziale . . . . 2.4 Ipotesi semplici per un test sequenziale . . . . . . . . . . . . . . . 2.4.1 Efficienza di un test sequenziale . . . . . . . . . . . . . . 11 13 14 14 14 3 SPRT con ipotesi semplici 3.1 Definizione di sequential probability ratio test . . . 3.2 Disuguaglianze utili per la scelta A e B in un SPRT 3.3 Scelta di A e B in casi pratici . . . . . . . . . . . . 3.4 Funzione OC per un sequential probability ratio test 3.5 ASN di un SPRT . . . . . . . . . . . . . . . . . . . . . . . 16 16 17 18 19 21 4 Test sequenziali di ipotesi unilaterali 4.1 Test di ipotesi semplici . . . . . . . . . . . . . . . . . . . . . . . 4.2 Test di ipotesi semplici contro alternative in un’unica direzione . . 23 23 24 5 SPRT caso variabile di Bernoulli 5.1 Presentazione del problema . . . . . . . . . . . . . . . . . . . . . 5.2 Implementazione del SPRT caso Bernoulli . . . . . . . . . . . . . 5.3 Funzione OC caso Bernoulli . . . . . . . . . . . . . . . . . . . . 25 25 26 29 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 7 9 9 10 11 6 Caso variabile Normale 6.1 Presentazione del problema . . . . . . . . . . . . . . . . . . . . . 6.2 Implementazione del SPRT caso normale . . . . . . . . . . . . . 6.3 Considerazioni sul numero di iterate necessarie al termine del test 6.4 Confronto SPRT con test di Neyman Pearson caso Normale . . . . 6.5 Applicazioni del SPRT al moto browniano . . . . . . . . . . . . . 6.6 Alcuni richiami sul moto browniano . . . . . . . . . . . . . . . . 6.7 Presentazione del problema . . . . . . . . . . . . . . . . . . . . . 6.8 Individuazione del cambiamento nel drift . . . . . . . . . . . . . 6.9 Tempi di individuazione . . . . . . . . . . . . . . . . . . . . . . 6.10 Analisi del metodo . . . . . . . . . . . . . . . . . . . . . . . . . 6.10.1 Sottostima . . . . . . . . . . . . . . . . . . . . . . . . . 6.10.2 Dipendenza dai parametri . . . . . . . . . . . . . . . . . 30 30 31 31 33 35 35 35 36 38 39 39 42 7 Conclusioni 46 8 Appendice 8.1 Codice matlab: Sequential Probability Ratio Test caso Bernoulli 8.2 Codice matlab: sequential probability ratio test caso Normale . . 8.3 Neyman test . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4 Test sul Moto Browniano . . . . . . . . . . . . . . . . . . . . . 47 47 48 49 49 3 . . . . Introduzione Questa tesi nasce motivata dal problema dell’individuazione del momento in cui un moto browniano cambia il suo drift. Ne deriva la necessità di studiare in modo più approfondito la distribuzione dei tempi necessari all’identificazione del cambiamento. La questione viene affrontata attraverso gli strumenti della teoria rigurdante la "Sequential Analysis". La teoria classica del disegno di esperimenti è principalmente basata su metodi a numerosità campionaria fissata. In questi test i dati vengono analizzati una sola volta e la dimensione del campione è prederminata. In alternativa a tali metodi si possono utilizzare test sequenziali. Quando si adotta questo tipo di analisi i dati sono, invece, monitorati man mano che vengono raccolti. Lo scopo principale è di trarre informazioni utili nel minor tempo possibile al fine di anticipare i tempi di scelta o di intervento nella modifica di design del test. Questa teoria permette, appunto, di effettuare test che abbiano il valore atteso del numero di campioni esaminati minore rispetto a test di ipotesi classici con i medesimi errori di prima e seconda specie. L’albore di questo tipo di approccio prende piede negli anni ’20 motivato dalla necessità del controllo qualità nelle linee di produzione dalle idee di Walter Shewhart, che successivamente diventerà il padre della più moderna teoria statistica sul controllo qualità. Il metodo sequenziale diventa di centrale importanza durante la seconda guerra mondiale poiché la produzione di massa di equipaggiamenti militari e navali richiedeva procedure di controllo qualità sempre più rapide e che garantissero comunque errori molto contenuti. In quegli anni Abraham Wald, con i membri del gruppo di ricerca statistica della Columbia University, stava affrontando questi problemi. Egli per primo introdusse la nozione iniziale di Sequential probability ratio test (SPRT) ispirato al lemma di Neyman Pearson. Ci sono numerose ragioni per cui si predilige un’analisi intervallare dei dati; esse si possono riassumere in economiche, etiche e amministrative. I sequential test nascono proprio per ragioni economiche nell’ambito del controllo qualità. L’esito positivo o negativo di un test si riflette, economicamente, nell’immissione immediata di un prodotto nel mercato o nell’arresto anticipato di esso dalla produzione. I metodi sequenziali conducono, generalmente, ad un 1 risparmio nel numero di campioni esaminati, quindi tempo e costi (anche costo opportunità). Le ragioni amministrative riguardano il bisogno di assicurare che l’esperimento sia condotto come pianificato, che i soggetti soddisfino i criteri di selezione e che i trattamenti siano applicati come indicato. A volte un controllo iniziale dei risultati potrebbe rivelare la presenza di problemi, dando la possibilità di rimediare prima di sprecare troppe risorse o denaro. Un altro aspetto è la necessità di controllare la validità delle assunzioni fatte nel pianificare l’esperimento; eventuali assunzioni inappropriate possono quindi essere corrette in tempo. Negli esperimenti su persone, l’etica spinge a tenere sotto controllo costante i risulati, per evitare che i pazienti siano sottoposti a trattamenti poco sicuri, inferiori o non efficaci. Ragioni etiche spingono inoltre a concludere l’esperimento il prima possibile, in modo da poter utilizzare la terapia migliore, se viene rilevata la superiorità, o per poter concentrare le energie su altri trattamenti promettenti in attesa di essere valutati, nel caso in cui non vengano rilevate differenze. Il contributo che questa tesi desidera apportare è quello di presentare i modo semplice i pilastri della moderna teoria "Sequential Analysis" e alcune applicazioni dei Sequential Probability Ratio Test. Nel primo capitolo daremo alcune definizioni basilari per lo sviluppo degli argomenti successivi. Nel secondo capitolo illustreremo la teoria sequenziale partendo dalla nozione di sequential test. Nel terzo ci soffermeremo su un particolare test sequenziale chiamato "Sequentisl Probability Ratio Test". Il quarto è un breve capitolo in cui si effettuano considerazioni sui test di ipotesi semplici contro alternative unilaterali. Nei capitoli finali vengono presentate alcune applicazioni del Sequential Probability Ratio Test. Tutti gli algoritmi sviluppati sono stati implementati facendo uso del software Matlab e in appendice si possono trovare i codici più importanti. 2 Capitolo 1 Teoria classica sui test di ipotesi: concetti di base Prima di introdurre la teoria riguardante l’analisi sequenziale è necessario fornire alcune definizioni ed importanti risultati statistici propedeutici alla comprensione dell’argomento. Per eventuali approfondimenti si rimanda al testo G. Casella, R.L. Berger "Statistical Inference", Duxbury Press, 2001. Definizione 1. Nel calcolo delle probabilità lo spazio campionario o insieme universo (generalmente indicato dalle lettera Ω) è l’insieme dei possibili risultati di un esperimento casuale. Definizione 2. Sia X una variabile aleatoria discreta, la sua densità discreta pX pxq “ PpX “ xq @x, ovvero una funzione di variabile reale che assegna ad ogni valore possibile di X la probabilità dell’evento elementare pX “ xq. Nel caso in cui la variabile casuale X sia continua, allora tale probabilità è sempre nulla, quindi questa definizione è inutile. Si utilizza la funzione di densità di probabilità. pX pxq descrive nel caso continuo la "densità" di probabilità in ogni punto nello spazio campionario ed è tale che ż PpX P Aq “ pX pxq d x A Definizione 3. Date due variabili aleatorie X e Y, definite sullo stesso spazio di probabilità, si definisce la loro distribuzione congiunta Fpx, yq “ PpX ď x, Y ď yq Nel caso di due sole variabili, si parla di distribuzione bivariata, mentre nel caso di più variabili si parla di distribuzione multivariata. La funzione di ripartizione di una distribuzione congiunta è definita come più generalmente Fpx1 , ..., xn q “ PpX1 ď x1 , ..., Xn ď xn q. 3 Nel caso di variabili aleatorie discrete, la densità discreta congiunta è data da PpX “ x, Y “ yq “ PpY “ y | X “ xq ¨ PpX “ xq “ PpX “ x | Y “ yq ¨ PpY “ yq Siccome la densità congiunta è anch’essa una densità, è soddisfatto ÿÿ PpX “ x , Y “ yq “ 1. x y Nel caso di variabili aleatorie continue, la densità congiunta è data da fX,Y px, yq “ fY|X py|xq fX pxq “ fX|Y px|yq fY pyq Anche in questo caso, è soddisfatto ż ż fX,Y px, yq dy dx “ 1. x y Definizione 4. Due variabili aleatorie si dicono indipendenti se PpX “ x, Y “ yq “ PpX “ xqPpY “ yq. Definizione 5. Sia n P R, un campione casuale di taglia n corrisponde a n realizzazioni px1 , x2 , ..., xn q delle variabili aleatorie Xi „ f p¨, θq indipendenti e identicamente distribuite. In questo caso chiameremo X „ f p¨, θq popolazione o variabile aleatoria generatrice. Definizione 6. Sia px1 , x2 , ..., xn q un campione casuale di taglia n estratto da una popolazione X „ f p¨, θq, la funzione di verosimiglianza Lpθq è definita come la densità congiunta delle variabili aleatorie X1 , X2 , .., Xn calcolata in px ś1 , x2 , ..., xn q con θ parametro libero. Dunque Lpθq “ fX1 ,X2 ,...,Xn px1 , x2 , ..., xn ; θq “ ni“1 fp xi ; θq In statistica quando si vuole vericare la bontà di un’affermazione che ha come oggetto accadimenti nel mondo reale(ipotesi), che si presta ad essere confermata o smentita dai dati osservati sperimentalmente, si valuta l’attendibilità di essa con il metodo del test di ipotesi statistico. Definizione 7. Un test di ipotesi statistico è una procedura statistica volta a decidere se accettare o meno un’ipotesi. Nella pratica si ha: • un’ipotesi nulla H0 che rappresenta lo status quo, cioè che viene rifiutata solo se c’è una evidenza in favore di H1 • un’ipotesi alternativa H1 che è quella che voglio effettivamente testare 4 Un test di ipotesi semplici è del tipo # H0 : θ “ θ0 H1 : θ “ θ1 un test di ipotesi composte è del tipo # H0 : θ “ θ0 H1 : θ ď θ1 un test di ipotesi bilaterali è del tipo # H0 : θ “ θ0 H1 : θ , θ1 La procedura del test consiste nell’associare ciascun campione appartenente allo spazio campionario totale alla decisione di scegliere se accettare o rifiutare l’ipotesi nulla. In sostanza viene scelta una regola che divide l’intero spazio campionario in due regioni: la regione di rifiuto e la regione di accettazione. • @x “ px1 , x2 , .., xn q P Cpregione criticaq Ă χn rifiuto l’ipotesi nulla, • @x “ px1 , x2 , .., xn q P C c pregione d1 accettazioneq Ă χn la accetto. Il test condurrà ad una decisione che non può essere in assoluto quella corretta. Vi sono due tipi di errori che possiamo commettere: di prima e seconda specie. Definizione 8. L’errore di prima specie rappresenta l’errore che si commette nel rifiutare l’ipotesi nulla quando questa è vera, mentre quello di seconda specie consiste nell’accettare H0 quando questa è falsa. Denotiamo con αpθq la Pperrore di prima specieq “ PpX P C; θ P Θ0 q e con βpθq la Pperrore di seconda specieq “ PpX P C c ; θ P Θ1 q Definizione 9. Si dice potenza del test la funzione Πpθq “ PpX P C; θq che rappresenta la probabilità di rifiuto al variare di θ, dunque è: # αpθq θ P Θ0 Πpθq “ 1 ´ βpθq θ P Θ1 L’estremo superiore di Πpθq per θ P Θ0 si dice ampiezza del test e si denota con α. Un buon test deve avere un’ampiezza molto piccola, generalmente compresa tra 0, 1 e 0, 5, per poter rifiutare lo status quo. 5 Definizione 10. Un test y di ampiezza α per H0 : θ “ Θ0 contro H1 : θ “ Θ1 si dice uniformemente più potente se @ test y1 di ampiezza α1 ě α si ha β1 pθq ě βpθq @ θ cioè la potenza Πpθq ě Π1 pθq @ θ Lemma 11. Lemma di Neyman Pearson [2] sia px1 , x2 , .., xn q un campione casuale estratto da una popolazione X „ f p¨, θqq il test uniformemente più potente di ampiezza α per il sistema di ipotesi semplici H0 : θ “ Θ0 contro H1 : θ “ Θ1 è quello con l seguente regione di rifiuto: C “ tpx1 , x2 , .., xn q | Lpθ0 ; px1 , x2 , .., xn qq ď ku Lpθ1 ; px1 , x2 , .., xn qq dove k è una costante determinata in modo che Ppx P Cq “ α 6 Capitolo 2 Test di ipotesi sequenziale 2.1 Nozione iniziale di test sequenziale Un test è condotto sequenzialmente quando i dati vengono sottoposti ad analisi intermedie, man mano che sono disponibili. Si contrappone a questa analisi quella a dimensione campionaria fissata, la quale esamina i dati una sola volta quando l’esperimento è concluso. Il disegno sequenziale è particolarmente indicato quando la raccolta dati viene effettuata gradualmente in un intervallo di tempo. Ad ogni analisi si decide se continuare o terminare l’esperimento in base ad una regola di arresto definita in fase di pianificazione del test. Quando si ottiene evidenza in favore di una ipotesi tale da rendere superflua la continuazione dell’esperimento, è possibile terminarlo anticipatamente. Il numero di iterate e dunque la dimensione del campione del test sequenziale non è predeterminato, ma una variabile aleatoria dipendente dalla regola d’arresto scelta e dall’outcome delle osservazioni. Più nel dettaglio: una regola viene data e ad ogni step dell’esperimento viene effettuata una di tre decisioni: 1. accettare l’ipotesi nulla 2. rifiutare l’ipotesi nulla 3. continuare l’esperimento facendo una nuova osservazione Il processo continua fino a quando non si effettua una scelta tra 1 o 2. Più formalmente: @ m P N D Mm ” totalitá di tutti i possibli campioni px1 , x2 , ..., xm q Pensiamo px1 , x2 , ..., xm q come ad uno spazio m-dimensionale e @ m partizioniamo questo spazio in tre regioni: R0m , R1m , Rm . Dopo la prima osservazione decidiamo se: 7 1. accettare l’ipotesi nulla se x1 giace in R0m 2. rifiutare l’ipotesi nulla se x1 giace in R1m 3. effettuare una seconda osservazione se x1 giace in Rm e decidiamo se: (a) accettare l’ipotesi nulla se px1 , x2 q giace in R0m (b) rifiutare l’ipotesi nulla se px1 , x2 q giace in R1m (c) effettuare una terza osservazione se px1 , x2 q giace in Rm Il processo, come si può intuire, termina se ci troviamo nelle regioni di spazio campionario R0m e R1m , mentre continua ad iterarsi se ci troviamo nella regione Rm . Dunque il test sequenziale è costituito dalla definizione di queste tre regioni. Se R0m , R1m , Rm sono in somma diretta è sufficiente definire due delle tre regioni. Definizione 12. px1 , x2 , ..., xm q è detto campione inefficace seŤ contiene un segmen1 R0m . Un campione to iniziale px1 , x2 , ..., xn q dove n ă m | px1 , x2 , ..., xn q P Rm non inefficace è detto campione efficace. Proposizione 13. Per effettuare un test sequenziale bisogna avere un campione efficace per ogni step dell’esperimento. Nella definizione delle tre regioni vanno dunque scartati tutti i campioni inefficaci. Esempio 14. Consideriamo un lotto di n prodotti. Ogni prodotto può essere difettoso o non difettoso e p (valore sconosciuto) è la proporzione di prodotti difettosi all’interno del lotto. Sia p1 un valore fissato, possiamo accettare il lotto se p ď p1 e rifiutarlo se p ą p1 . Sia n0 P N fissato: • se le prime n0 unita sono non difettose allora fermiamo l’ispezione e accettiamo il lotto, • se per qualche m ď n0 l’m-esima unità e difettosa allora fermiamo l’ispezione e rifiutiamo il lotto. Assegniamo il valore 1 se il prodotto è difettoso e 0 se non difettoso. Si continua l’ispezione se e solo se m ď n0 e inoltre, x1 “ x2 “ .... “ xm´1 “ 0. Dalle definizioni precedenti possiamo dire che: • R0m non contiene campioni efficaci per m ă n0 ; l’unico campione efficace si ha con m “ n0 , dunque R0n0 “ p0, 0, ..., 0q dove il valore 1 si trova nel posto n0 -esimo. 8 • R1m contiene n0 campioni di taglia m, ovvero { p0, 0, ..., 1q | il valore 1 si trova all’m-esimo posto con m “ 1, 2, ...n0 } • Rm contiene n0 campioni di taglia m, ovvero {p0, .., 0q | con m “ 1, 2, .., n0 } Possiamo scegliere i set delle tre regioni in molteplici modi. Il fondamentale problema, difatti, dei test di ipotesi sequenziali è effettuare una corretta scelta di questi set. 2.2 Conseguenze della scelta di ogni particolare test 2.2.1 La funzione Operating Characteristic (OC) Sia X la variabile aleatoria da cui estraggo il campione, X „ f px, θ1 , θ2 , ....θk q la cui forma funzionale è conosciuta, mentre un numero finito di parametri è sconosciuto. Effettuiamo una particolare scelta dei set di R0m , R1m e Rm Definizione 15. La funzione operating characteristic (OC) equivale alla probabilità che il processo termini con l’accettazione dell’ipotesi nulla H0 . Questa funzione dipende solo dalla distribuzione della variabile aleatoria X, dunque solamente dai parametri pθ1 , θ2 , ....θk q. Per fare riferimento a questa funzione useremo la notazione Cpθq “funzione OC e se θ è unico sarà possibile rappresentarla sugli assi cartesiani. Si può rilevare che la funzione operating charateristic è il coniugato della funzione di potenza del test richiamata nella definizione (9): Cpθq “ 1 ´ Πpθq Difatti, considerando solo i test che terminano con P “ 1 abbiamo che la probabilità di rifiuto è esattamente 1 ´ Cpθq. Esempio 16. Riferiamoci ora all’esempio 14. Abbiamo θ “ p (con p “ proporzione prodotti difettosi) unico parametro. Inoltre, immaginiamo che il lotto sia infinitamente più grande della taglia n0 del campione; in questo modo Cpθq “ p1 ´ pqn0 . Supponiamo di far variare θ lungo l’asse orizzontale e Cpθq lungo quello verticale e otteniamo la rappresentazione in Figura 1 9 1 Cpθq 0.8 0.6 0.4 0.2 0 θ 0.2 0.4 0.6 0.8 1 Figura 2.1: Questa funzione dà alcune informazioni importanti: la probabilità di effettuare una corretta decisione per ogni parametro θ. Se θ è consistente con H0 varrà Cpθq, mentre per ogni parametro non consistente 1 ´ Cpθq. Dunque per θ consistenti con l’ipotesi nulla vogliamo Cpθq vicino a 1, per quelli non consistenti, invece, Cpθq vicino a 0. 2.2.2 La funzione Average Sample Number (ASN) Precedentemente abbiamo rilevato che il numero di osservazioni da effettuare per portare a termine un test sequenziale è una variabile aleatoria. Ora vogliamo esprimere attraverso una funzione la costosità dell’esperimento in termini di numero di osservazioni e dunque di tempo. Naturalmente è preferibile costruire un test che esca in un tempo finito non eccessivamente grande. Definizione 17. Sia θ il vettore di parametri della variabile aleatoria X che si vuole testare, la funzione average sample number(ASN)“ Eθ pnq equivale al valore atteso del numero di osservazioni; essa dipenderà, difatti, da un numero finito di parametri della distribuzione che prendiamo in considerazione e ci servirà proprio per capire, al variare dei parametri, quante osservazioni saranno effettuate. Esempio 18. Facciamo riferimento sempre all’Esempio precedentemente discusso. Il numero di osservazioni può variare da 1 a n0 . Chiamiamo n il numero di osservazioni e Eθ pnq il valore atteso di n. Sappiamo che: • se durante l’ispezione non trovo neanche un campione difettoso allora terminerò il test dopo n0 osservazioni, dunque la probabilità di terminare all’n0 esimo step vale p1 ´ pqn0 ´1 • se, invece, le prime m ´ 1 unità sono non difettose e l’m-esima è difettosa con m ď n0 allora m è il numero di osservazioni. La probabilità di terminare all’m-esimo step @ m vale p1 ´ pqm´1 . 10 10 9 8 7 6 5 4 3 2 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 2.2: rappresenta la funzione ASN. Mentre la funzione OC determinava quanto bene il test riusciva a fare una correta decisione, quest’ultima, invece, ci dice quanto il test è costoso in termini di numero di osservazioni al variare del parametro θ. Eθ pnq “ nÿ o ´1 m“1 mpp1 ´ pqm´1 ` n0 p1 ´ pqn0 ´1 pFig.2q 2.3 Principi di selezione di un sequential test 2.3.1 Livello di preferenza per accettazione e rifiuto dell’ipotesi nulla Sia θ il parametro sconosciuto della variabile aleatoria su cui stiamo effettuando il test, la preferenza per il rifiuto o l’accettazione è una funzione che dipende da θ. Chiamiamo ω la regione di tutti i parametri consistenti con l’ipotesi nulla H0 . In altri termini H0 è l’affermazione: il parametro reale θ P ω. Se consideriamo il caso in cui il parametro sconosciuto sia solo uno e l’ipotesi nulla è H0 : θ ⩽ θ0 , dunque ω “ tθ ⩽ θ0 u ñ è preferibile accettare H0 se θ P ω, mentre è preferibile rifiutare se θ P ω̄. Inizieremo, ora, ad introdurre livelli di preferenza di accettazione o rifiuto. Se il valore reale giace in ω, ma è vicino la frontiera, la preferenza per l’accettazione di H0 è piccola; se giace in ω̄, ma vicino la frontiera, la preferenza per il rifiuto di H0 è piccola anche in questo caso. Possiamo notare che non si commetterebbe un grosso errore nel rifiutare H0 quando θ P ω o nell’accettare H0 quando θ P ω̄. 11 Definizione 19. La frontiera è l’insieme di tutti i punti appartenenti ad un intorno di θ e che quindi giaciono sia in ω che in ω̄. Possiamo, in generale, dividere lo spazio campionario anche qui in tre regioni: • forte preferenza per l’accettazione dell’ipotesi nulla ( regione Ă ω ), • forte preferenza per il rifiuto dell’ipotesi nulla ( regione Ă ω̄ ), • indifferenza per rifiuto o accettazione ( regione vicino ai punti di frontiera ). La separazione in tre regioni non si tratta come un problema statistico, ma si effettua sulla base di considreazioni riguardanti le conseguenze delle scelte sbagliate. Più precisamente: # ω0 pθq “ importanza dell’errore causato dall’accettazione di H0 quando θ é vero. ω1 pθq “ importanza dell’errore causato dal rifiuto di H0 quando θ é vero. Si ottiene che: ω0 pθq “ 0 @θ Pω e ω1 pθq “ 0 @ θ P ω̄ ω0 pθq ě 0 @ θ P ω̄ e ω1 pθq ě 0 @θ Pω Inoltre crescono i valori ω0 e ω1 all’allontanarsi dalla frontiera, mentre @ θ P alla regione di indifferenza ω0 “ ω1 “ 0. Facciamo ora qualche esempio per comprendere meglio il significato di preferenza. Esempio 20. Consideriamo un’ispezione su un lotto di prodotti, classificabili in difettosi e non difettosi, la cui proporzione di prodotti difettosi sia definita p. Generalmente è possibile trovare p0 e p1 tali che se p ă p0 il livello di preferenza per l’accettazione del lotto è alto e se p ą p1 il livello di preferenza per il rifiuto è elevato. Dunque la regione compresa tra p0 e p1 sarà proprio la regione di indifferenza. 12 Esempio 21. Esaminiamo un prodotto che deve essere immesso sul mercato; la durezza di questo prodotto si distribuisce come una variabile aleatoria normale X „ Normpµ, σq con σ conosciuta. µ0 è il nostro valore desiderato per poterlo immettere. Per prendere una corretta decisione dobbiamo verificare che la media campionaria θ non discosti molto dal valore desiderato µ0 : |θ ´ θ0 | ă c, con c costante. E’ possibile anche in questo caso trovare un ∆| accetto con un livello di preferenza significativo se |θ ´ θ0 | ă c ´ ∆ e rifiuto se |θ ´ θ0 | ą c ` ∆. Se mi trovo nella zona compresa tra c ´ ∆ e c ` ∆ sono proprio nella regione di indifferenza. Esempio 22. 2.3.2 Condizioni imposte dall’OC function Dobbiamo affrontare alcune considerazioni teoriche che riguardano l’OC function; naturalmente quello che vogliamo che accada è trovare le seguenti condizioni: # H0 pθq : θ P ω Cpθq “ 1 @ θ H1 pθq : θ P ω̄ Cpθq “ 0 @ θ Questo è solamente un caso ideale che non può avere riscontri nella realtà, ma più l’OC function è simile alla funzione rappresentata in Figura 2.3 migliore sarà il test. Oltre a voler una OC function che sia simile a quella ideale, vogliamo che il valore atteso del numero di campioni da estrarre sia il minimo possibile. Molto spesso queste due richieste sono in rapporto di trade-off. Possiamo agire in questo modo: formuliamo prima i requisiti di vicinanza dall’OC ideale e consideriamo solo quei test che soddisfano tali requisiti, poi selezioniamo quello che ha valore atteso di campionamento minore. Affrontiamo il problema di formulare i requisiti di vicinanza. Dopo aver suddiviso lo spazio campionario in tre regioni, imponiamo condizioni sul comportamento di Cpθq nelle regioni di accettazione e di rifiuto, mentre non c’è necessità di apporre alcuna condizione nella regione di indifferenza. Le condizioni da imporre sono: # 1 ´ Cpθq ă α @θ P zona di accettazione (2.1) Cpθq ě β @θ P zona di rifiuto La scelta delle costanti arbitrarie α e β è da trattare non come un problema statistico, ma di design del test. Un test sequenziale si dice ammissibile se soddisfa le condizioni precedenti. 13 1.5 Cpθq 1 0.5 0 θ 0.2 0.4 0.6 0.8 1 Figura 2.3: funzione OC ideale . 2.3.3 ASN come base per la scelta di un test sequenziale Dopo aver diviso lo spazio dei parametri in tre zone, scelto α e β, selezioniamo esclusivamente i test che soddisfano le condizioni (2.1). Tra questi vogliamo utilizzare il test con Eθ pnq il piú piccolo possibile. Denotiamo con Eθ pn|S q ( definito in precedenza [pag.10] ) il valore atteso del numero di campioni da esaminare per terminare il test sequenziale S. Sia minS Eθ pn|S q il minimo valore di Eθ pn|S q rispetto a S fissato un particolare θ. Si ha: Eθ pn|S 1 q ě min Eθ pn|S q @S 1 S Definizione 23. il test si dice uniformemente migliore se D S 0 | Eθ pn|S q “ minS Eθ pn|S q @θ. Generalmente non é possibile minimizzare Eθ pn|S q @ θ, ma solamente trovare dei compromessi. 2.4 Ipotesi semplici per un test sequenziale 2.4.1 Efficienza di un test sequenziale Dobbiamo considerare solo due dei valori del parametro θ: θ0 e θ1 . Ciò che vogliamo testare è: # H0 : θ “ θ0 l’ipotesi nulla ovvero quella da confutare (2.2) H1 : θ “ θ1 l’ipotesi alternativa 14 In ogni test sequenziale di ipotesi semplici si associano due valori α e β tali che: Ñ se H0 è vera α è la probabilità di un errore di 1a specie (rifiutare H0 ). Ñ se H0 è falsa β è la probabilità di un errore di 2a specie (accettare H0 ). Proposizione 24. Due test S ed S 1 hanno la stessa potenza se α e β sono uguali a α1 e β1 . Se α ă α1 e β ď β1 oppure α ď α1 e β ă β1 S è detto più potente di S 1 , se, inevce, α ă α1 e β ą β1 oppure α ą α1 e β ă β1 i due test si dicono di potenza non comparabile. Restringiamoci sui test con (α, β) assegnati. Il test migliore sarà quello con Eθ pnq minore. Siano S e S 1 due test di uguale potenza: se Eθ0 pn|S q ď Eθ0 pn|S 1 q e Eθ1 pn|S q ă Eθ1 pn|S 1 q oppure Eθ0 pn|S q ă Eθ0 pn|S 1 q e Eθ1 pn|S q ď Eθ1 pn|S 1 q il test S è da preferire a S 1 . Se D S 0 tale che Eθ0 pn|S 0 q ď Eθ0 pn|S q e Eθ1 pn|S 0 q ď Eθ1 pn|S q @ S di uguale potenza allora S 0 é detto test ottimo. Più specificamente, denotiamo con n0 pα, βq il minimo valore di Eθ0 pn|S q e con n1 pα, βq il minimo valore di Eθ1 pn|S q. Sappiamo che Eθ0 pn|S q ě n0 pα, βq e Eθ1 pn|S q ě n1 pα, βq @ S . Un sequential test è ottimo quando vale l’ugualianza, ma in generale l’esistenza di tale ottimo non può essere provata. Definiamo, dunque, l’efficienza di un test S(α, β) attraverso questi rapporti: n1 pα, βq n0 pα, βq , Eθ0 pn|S q Eθ1 pn|S q Il primo rapporto vale quando H0 è vera e il secondo quando è falsa. Naturalmente l’efficienza giacerà tra 0 e 1. 15 Capitolo 3 SPRT con ipotesi semplici 3.1 Definizione di sequential probability ratio test Sia f px, θq la funzione di densità della variabile aleatoria X, consideriamo il seguente test: # H 0 : θ “ θ0 H 1 : θ “ θ1 Dunque la distribuzione sarà: f px, θ0 q se H0 è vera, mentre f px, θ1 q se H1 vera. La probabilità di ottenere le osservazioni successive x1 , x2 , x3 , ...., xn è: p1m “ f px1 , θ1 q f px2 , θ1 q.... f pxn , θ1 q se H1 é vera p0m “ f px1 , θ0 q f px2 , θ0 q.... f pxn , θ0 q se H0 é vera Il sequential probability ratio test sarà definito così: vengono scelti due valori A e B P R e ad ogni step dell’esperimento (dunque all’m-esimo campionamento) decido se accettare l’ipotesi nulla, rifiutarla o continuare l’esperimento. Se p1m ăA (3.1) Bă p0m continuo ad iterare l’esperimento facendo una nuova osservazione, mentre se: p1m ěA (3.2) p0m il processo terminerà con l’accettazione di H1 e quindi il rifiuto di H0 , infine se: p1m ďB (3.3) p0m il processo terminerà con l’accettazione di H0 ovvero l’ipotesi nulla. Nel paragrafo successivo verrà trattato il problema di come scegliere A e B per un test con errori di primo e seconda specie rispettivamente α e β. 16 3.2 Disuguaglianze utili per la scelta A e B in un SPRT In questa sezione giungeremo a trovare dei metodi per scegliere opportuni A e B per test con erroridi primo e seconda specie α e β attraverso l’utilizzo di alcune semplici disuguaglianze. Definiamo un campione px1 , x2 , .., xn q di tipo 0 se soddisfa: p1m p1m ă A per m “ 1, 2, ..., n ´ 1 e ď B per m “ n p0m p0m e di tipo 1 se soddisfa: p1m p1m ă A per m “ 1, 2, ..., n ´ 1 e ě A per m “ n Bă p0m p0m Bă Sappiamo che un campione di tipo 0 induce all’accettazione dell’ipotesi nulla (p “ p0 ) ed un campione di tipo 1 induce al suo rifiuto (p “ p1 ). Quando ho un campione di tipo 1, dalla disuguaglianza p3.2q si ha: (3.4) p1n ě Ap0n dove p0n rappresenta, in questo caso, la probabilità di accettare H1 quando H0 è vero (errore di prima specie α). p1n è, invece, la probabilità di accettare H1 quando quest’ultimo è vero (1 ´ β con β errore di seconda specie). Giungiamo quindi alla seguente affermazione: 1 ´ β ě Aα ñ Aď 1´β α (3.5) 1´β cioè è un limite superiore per A. α Il discorso è analogo per un campione di tipo 0; fatte le opportune considerazione possiamo asserire che: β β ď Bp1 ´ αq ñ B ě (3.6) 1´α β cioè è un limite superiore per B. 1´α Dalla disuguaglianza p3.5q possiamo ricavare: mentre dalla p3.6q: 1 α ď 1´β A ñ αď 1 A (3.7) β ďB 1´α ñ βďB (3.8) 17 3.3 Scelta di A e B in casi pratici Dopo aver scritto questa serie di disuguaglianze utili alla comprensione riportiamo i risultati di Wald nel libro "Sequential Analysis" nella sezione 3.2. [1]: A“ 1´β “ apα, βq α e B“ β “ bpα, βq 1´α mentre i veri valori sono Apα, βq e Bpα, βq. Ora si vuole giungere alla conclusione che l’utilizzo di apα, βq e bpα, βq, in luogo dei reali valori, garantisce la stessa accuratezza. Naturalmente la scelta di apα, βq e bpα, βq cambierà le probabilità degli errori di prima e seconda specie. Svolgiamo prima una semplice ragionamento: • se A fosse maggiore di Apα, βq e B fosse uguale a Bpα, βq la probabilità di errore di prima specie sarebbe minore, mentre quella di seconda specie lievemente maggiore, • se B fosse minore di Bpα, βq e A fosse uguale a Apα, βq la probabilità di errore di seconda specie sarebbe minore, mentre quella di prima specie lievemente maggiore, • se A fosse maggiore di Apα, βq e B fosse minore di Bpα, βq non è chiaro l’effetto sugli errori. Ponendo α1 e β1 come gli errori di prima e seconda specie per un test con A “ apα, βq e B “ bpα, βq segue dalle disuguaglianze p3.6q e p3.7q: α1 α 1 “ ď 1 ´ β1 1´β apα, βq β β1 ď bpα, βq “ 1 ´ α1 1´α possiamo cosí ricavare due altre importanti relazioni: α α1 ď 1´β (3.9) (3.10) (3.11) β (3.12) 1´α Moltiplichiamo ora la disuguaglianza p3.9q ambo i membri per p1 ´ βqp1 ´ β1 q e la disuguaglianza p3.10q per p1 ´ αqp1 ´ α2 q. Semplificandole e sommandole otteniamo: α1 ` β1 ď α ` β (3.13) β1 ď 18 Poiché solitamente α e β sono valori compresi tra 0.01 e 0.05 e quindi i valori β α e sono rispettivamente molto vicini ad α e β, per p3.11q e p3.12q 1´α 1´β possiamo ignorare di quanto può eccedere α1 rispetto ad α e β1 rispetto a β. Inoltre per la p3.13q almeno uno tra α1 ď α e β1 ď β è vera. Possiamo concludere che a fini pratici l’utilizzo di A “ apα, βq e B “ bpα, βq garantisce almeno la stessa accuratezza e tutela da decisioni sbagliate rispetto ad un test con A “ Apα, βq e B “ Bpα, βq. 3.4 Funzione OC per un sequential probability ratio test In questa sezione si affronta il problema di calcolare la funzione Cpθq in un SPRT. La funzione OC risulta particolarmente importante poiché, sebbene i test che sino ad ora abbiamo esaminato sono ad ipotesi semplici, i reali valori dei parametri da testare possono assumere anche valori diversi da quelli ipotizzati. Ricordiamo che la funzione Cpθq rappresenta al variare di θ la probabilità di terminare il test con l’accettazione dell’ipotesi nulla H0 . Consideriamo l’espressione „ ȷ f px, θ1 q hpθq f px, θ0 q per ogni valore di θ, il valore hpθq , 0 é scelto tale che: „ ż `8 „ ȷ „ f px, θ1 q f px, θ1 q hpθq E f px, θqdx “ 1 f px, θ0 q ´8 f px, θ0 q e se la variabile aleatoria X in considerazione è discreta ÿ „ f px, θ1 q ȷhpθq f px, θqdx “ 1 f px, θ0 q x Possiamo per tanto dire che la funzione „ ȷ f px, θ1 q hpθq ˚ f px, θq “ f px, θq f px, θ0 q (3.14) (3.15) (3.16) è una funzione di densità di una variabile aleatoria. Poichè hpθq , 0 si devono trattare i casi hpθq ă 0 e hpθq ą 0. Iniziamo con hpθq ą 0. 19 Utilizziamo il sequential probability ratio test S ˚ per testare H, ovvero l’ipotesi che f px, θq sia la distribuzione di x contro H˚, l’ipotesi che f ˚ px, θq sia la distribuzione di x. Fino a quando Bhpθq ă f ˚ px1 , θq f ˚ px2 , θq..... f ˚ pxm , θq ă Ahpθq f px1, θq f px2 , θq..... f pxm , θq (3.17) il test S ˚ continua ad iterarsi. Accetterò H se: f ˚ px1 , θq f ˚ px2 , θq..... f ˚ pxm , θq ď Bhpθq f px1, θq f px2 , θq..... f pxm , θq (3.18) f ˚ px1 , θq f ˚ px2 , θq..... f ˚ pxm , θq ě Ahpθq f px1, θq f px2 , θq..... f pxm , θq (3.19) mentre rifiuterò H se: Dal momento che f ˚ px, θq “ f px, θq „ f px, θ1 q f px, θ0 q ȷhpθq sostituendo nella disuguaglianza p3.17q e p3.16q otteniamo Bă f px1 , θq f px2 , θq..... f pxm , θq ăA f px1, θq f px2 , θq..... f pxm , θq (3.20) Così se il test S ˚ conduce all’accettazione di H il test S condurrà all’accettazione di H0 . Ne deduciamo che i due test hanno regioni di accettazione, rifiuto e indifferenza equivalenti. Denotiamo con α1 e β1 gli errori di prima e seconda specie del test S ˚. Desumiamo dalle disuguaglianze p3.5q e p3.6q che Ahpθq ď 1 ´ β1 α1 Bhpθq ě β1 1 ´ α1 e Per le considerazioni effettuate nella sezione precedente, si può approssimare Ahpθq e Bhpθq senza causare errori ai fini pratici. Ahpθq „ 1 ´ β1 α1 Bhpθq „ β1 1 ´ α1 e 20 Da queste ultime due disuguaglianze otteniamo: α1 „ 1 ´ Bhpθq Ahpθq ´ Bhpθq (3.21) e poichè α1 “ 1 ´ Cpθq Ahpθq ´ 1 (3.22) Ahpθq ´ Bhpθq Per la dimostrazione rigorosa ed il calcolo in modo non approssimato della funzione OC si rimanda all’articolo "Sequential Test of Statistical Hypothesis" scritto su "The Annals of Mathematical Statistics" nel 1951 [3]. Cpθq „ 3.5 ASN di un SPRT Come abbiamo osservato nel primo capitolo, il pregio principale del SPRT è quello di utilizzare un numero di ossevazioni generalmente inferiore rispetto ad analoghi test con il numero di campioni fissato. Si vuole ora dare una dimostrazione teorica che provi la suddetta affermazione. Sia n “ numero di osservazioni e Eθ pnq “ valore atteso di n. Quest’ultimo valore rappresenta proprio l’ASN function che dipende dal paramtro θ e l’obiettivo è trovare una buona approssimazione p1m non considerando gli eccessi di oltre la frontiera A e B al termine del prop0m cesso. Il risultato non approssimato è possibile trovarlo nell’ appndice del libro "Sequential Analysis" di Wald. Consideriamo N un intero molto grande tale che Ppn ď Nq è trascurabile. Possiamo scrivere: z1 ` z2 ` ... ` zN “ pz1 ` .. ` zn q ` pzn ` 1 ` .. ` zN q dove: zα “ logp (3.23) f pxα , θ1 q f pxα , θ0 q Se si applica il valore atteso ad entrambi i membri di (3.23) si ottiene: NEpzq “ Epz1 ` ... ` zn q ` Epzn`1 ` ... ` zN q (3.24) Se le variabili zα sono i.i.d. allora segue che: Epzn`1 ` ... ` zn q “ EpN ´ nqEpnqEpzq “ NEpzq ´ EpnqEpzq (3.25) quindi: Epz1 ` ... ` zn q ´ EpnqEpzq “ 0 21 (3.26) Epnq “ Epz1 ` .... ` zn q Epzq (3.27) trascurando gli eccessi oltre A e B si ha: Epz1 ` .... ` zn q „ Cpθq log B ` p1 ´ Cpθqq log A (3.28) sostituendo in (3.27) Eθ pnq „ Cpθq log B ` p1 ´ Cpθqq log A Eθ pzq 22 (3.29) Capitolo 4 Test sequenziali di ipotesi unilaterali Precedentemente ci siamo concentrati sui test di ipotesi semplici H0 contro una singola alternativa H1 . Nelle applicazioni alla realtà, risulta molto più probabile che il parametro da testare possa assumere infiniti altri valori. In questo capitolo affronteremo il tema dei test sequenziali di ipotesi semplici o composte contro alternative unilaterali. 4.1 Test di ipotesi semplici Consideriamo il caso in cui vi sia un solo sconosciuto parametro θ coinvolto nella distribuzione della variabile aleatoria X. Una ipotesi semplice è che θ sia uguale ad un valore specifico θ0 . Ci sono due alternative (decidere di accettare“alternativa 1; rifiutare“alternativa 2 ) e la preferenza per una di esse dipende dal vero valore di θ. Definiamo ω la regione di tutti i valori di θ per i quali l’alternativa 1 è preferita alla 2, dunque l’alternativa 2 è preferita alla 1 per tutti i valori θ al di fuori di ω. Allora il problema della decisione tra le due alternative può essere riformulato come un test di ipotesi Hw . Se Hw è accettato prendimo la decisione 1 e se Hw è rifiutata prendiamo l’alternativa 2. Supponendo che il livello di preferenza per una o l’altra delle alternative vari continuamente con il valore θ, la regione ω non può consistere in un solo valore θ0 . Difatti, se ω contenesse un singolo valore θ0 vorrebbe dire che preferiamo l’alternativa 1 quando θ “ θ0 e la 2 per ogni θ , θ0 non importa quanto θ sia vicino a θ0 . Notiamo che il test di ipotesi semplici ha un senso soltanto quando vi è una discontinuità nella nostra preferenza per le due alternative. Un esempio di discontinuità per la preferenza si ha nel caso in cui si voglia testare, per esempio, la validità di una teoria scientifica riguardante un determinato parametro θ che deve 23 assumere uno specifico valore θ0 . Nelle applicazioni alla realtà è comunque molto raro avere delle discontinuità. 4.2 Test di ipotesi semplici contro alternative in un’unica direzione In questa sezione analiziamo un semplice caso: H0 : θ “ θ0 contro H1 : θ ą θ0 . La zona di preferenza per l’accettazione dell’ipotesi nulla consiste solamente nel valore θ0 . Il livello di preferenza per il rifiuto decrescerà all’aumentare del valore θ nel dominio θ ą θ0 . Possiamo a questo punto trovare un valore θ1 ą θ0 tale che l’accettazione dell’ipotesi è considerata un errore di importanza rilevante. Dunque, la zona di preferenza per il rifiuto si definisce da θ ą θ1 e la zona di indifferenza da θ0 ă θ ă θ1 . Imponiamo dei requisiti sul comportamento dell’ OC function: la probabilità che l’ipotesi sia rifiutata è uguale ad un predetermintato valore α quando θ “ θ0 e uguale a β quando θ ě θ1 . Nella maggior parte dei casi come per esempio per la variabile aleatoria binomiale, di Poisson e normale sequential probability ratio test di potenza pα, βq per testare l’ipotesi θ “ θ0 contro la singola alternativa θ1 che soddisfa tali requisiti, fin quando la probabilità di errore di seconda specie decresce all’aumentare di θ nel dominio θ ą θ0 . Infine il sequential probability ratio test sopra citato è una soddisfacente soluzione al nostro problema. 24 Capitolo 5 SPRT per variabili aleatorie 5.1 SPRT caso Bernoulli 5.1.1 Presentazione del problema Sia X una variabile aleatoria che può assumere solamente i valori 0 e 1. Sia p la probabilità che X assuma il valore 1 e 1´ p quella che assuma il valore 0. Si vuole testare l’ipotesi che p sia minore di un determinato valore p. Nella realtà si ha un riscontro nei problemi attinenti al controllo qualità; ad esempio quando si vuole ispezionare un lotto composto da un numero molto grande di prodotti che possono essere classificati in due categorie (difettosi e non difettosi). Si può assegnare il valore 0 ai prodotti non difettosi ed il valore 1 a quelli difettosi. Il lotto sarà da immettere nel mercato, quindi accettato, se la proporzione p di prodotti difettosi non eccede un determinato valore p. Si accetta, quindi, il lotto se p ď p e si rifiuta quando p ą p. Il problema in questione si può formulare in questo modo: # H0 : p ď p (5.1) H1 : p ą p Poiché una completa ispezione del lotto sarebbe troppo costosa, si accetta di tollerare qualche rischio nel prendere una decisione errata. Se p ą p il rifiuto del lotto è preferito e questa preferenza cresce all’aumentare di p, mentre per p ă p si preferisce accettare il lotto e la preferenza per l’accettazione cresce al diminuire del valore p. Se p “ p si è indifferenti alla scelta da effettuare. Se il valore p è maggiore di p, ma p ´ p è un numero sufficientemente piccolo, l’errore in caso di accettazione non è rilevante; e così pure nel caso di p ă p. Sarà possibile allora specificare due valori p0 e p1 tali che il rifiuto del lotto comporti un errore di rilevante importanza se p ď p0 e l’accettazione comporti errori rilevanti per p ě p1 . per valori di p compresi tra p0 e p1 nessuna decisione è preferibile. 25 La nostra propensione al rischio di errore si misura attraverso i valori: α “ probabilità di rifiutare il lotto quando p ď p0 e β “ probabiltà di accettare il lotto quando p ě p1 . I valori p0 , p1 , α, β sono determinati in base a considerazione di natura non statistica. La scelta di p0 e p1 dipende dall’importanza dell’errore causato da una scelta errata al variare del vero valore p, mentre la scelta di α e β dipende, come detto in precedenza, sia dai costi in termine di tempo dell’ispezione che dal livello di accuratezza che si vuole ottenere. L’applicazione al controllo qualità con una variabile aleatoria qui presentato è un caso semplice. Per capire meglio quali possano essere le applicazioni dell’analisi statistica al campo del controllo qualità si veda il libro "Economic Control of Quality Manufactured Product. American Society for Quality (ASQ)" di Shewhart Walter. 5.1.2 Implementazione del SPRT caso Bernoulli Ipotizziamo di aver trovato i nostri valori p0 , p1 , α, β in base alle considerazioni prima effettuate. Un test che soddisfa le condizioni: Ppx P C; p ď p0 q ď α e Ppx P C c ; p ě p1 q ď β è il sequential probability ratio test di potenza pα, βq che testa le ipotesi # H0 : p “ p0 (5.2) H1 : p “ p1 poiché la funzione OC di un sequential probability ratio test è monotona decrescente. Il test è definito così: Sia xi il risultato della i-esima ispezione e p la proporzione di prodotti difettosi nel lotto allora la funzione di verosimiglianza f px, pq di un campione di taglia m (che chiamiamo per comodità pm ) è data da: p sm p1 ´ pqm´sm dove sm “ Quando p “ p1 la probabilità (23) è: m ÿ (5.3) xi i“1 p1m “ p sm p1 ´ p1 qm´sm (5.4) p0m “ p sm p1 ´ p0 qm´sm (5.5) mentre quando p “ p0 si ha: Il logaritmo del rapporto di p0m e p1m equivale a: logp p1m p1 1 ´ p1 q “ sm log ` pm ´ sm q log p0m p0 1 ´ p0 26 (5.6) dunque il test continua a campionare fino a quando: log p1m 1´β β ă logp q ă log 1´α p0m α (5.7) Si rifiuterà l’ipotesi nulla quando: logp p1m 1´β q ě log p0m α (5.8) logp p1m β q ď log p0m 1´α (5.9) mentre si accetterà se Vogliamo ora riscrivere le tre disuguaglianze precedenti in modo più comodo: 1 ´ p0 1 ´ p0 β 1´β log log log 1 ´ p1 1 ´ p1 1´α α `m ă sm ă `m p1 p1 p1 p1 1 ´ p1 1 ´ p1 1 ´ p1 1 ´ p1 log log log log ´ log ´ log ´ log ´ log p0 1 ´ p0 p0 1 ´ p0 p0 1 ´ p0 p0 1 ´ p0 log e per praticità chiamiamo 1 ´ p0 1´β log 1 ´ p1 α `m dm “ p1 1 ´ p1 p1 1 ´ p1 log ´ log log ´ log p0 1 ´ p0 p0 1 ´ p0 (5.10) 1 ´ p0 β log 1 ´ p1 1´α am “ `m p1 1 ´ p1 p1 1 ´ p1 log ´ log log ´ log p0 1 ´ p0 p0 1 ´ p0 (5.11) log e log In questo modo è possibile tracciare un grafico che spiega l’andamento del test: i valori am sono punti di una retta sopra la quale ci troviamo nella regione di accettazione e i punti rm definiscono una retta sotto la quale si ha la regione di rifiuto. I valori compresi tra le due rette sono appartenenti alla regione di indifferenza. Fino a quando non si ha una intersezione con una delle rette disegnate a partire dai punti di dm o am il test continua a procedere. 27 14 12 10 Prodotti difettosi 8 Zona di rifiuto 6 4 2 0 -2 -4 zona di accettazione 0 5 10 15 20 25 30 Numero di iterate Figura 5.1: Il grafico mostra lo svolgimento del SPRT caso bernoulli di ipotesi H0 : p “ 0.2 contro H1 : p “ 0.5 con α e β errori di prima e seconda specie equivalenti rispettivamente a 0.01 e 0.05. La linea gialla rapprensenta il numero di prodotti difettosi trovati sino a quel momento. Al di sotto della linea blu abbiamo la regione di accettazione, al di sopra di quella rossa la regione di rifiuto, mentre per i valori compresi tra queste due linee si ha la regione di indifferenza. Il test continua procedere finchè la linea gialla non interseca una delle altre due linee. 28 5.1.3 Funzione OC caso Bernoulli Per completezza diamo ora una buona approsimazione della funzione OC, la quale descrive la probabilità di accettare l’ipotesi nulla al variare del parametro da testare. In questo modo si confermano le considerazioni fatte nel capitolo p3q. Per definizione sappiamo che 1 ´ α “ Cpp0 q e β “ Cpp1 q. Dalla disuguaglianza p3.15q deriva p p1 hpθq 1 ´ p1 hpθq ` p1 ´ pq “1 p0 1 ´ p0 (5.12) Considerando h un parametro risolviamo la precedente equazione rispetto a p: 1 ´ p1 h 1 ´ p0 p“ p1 h 1 ´ p1 h ´ p0 1 ´ p0 1´ (5.13) Facendo variare i valori h possiamo numericamente trovare i valori p. In aggiunta dalla p3.22q abbiamo: 1 ´ βh ´1 α Cppq „ 1 ´ βh β h ´ α 1´α pβq p1 ´ βq eB“ poichè A “ α 1´α Possiamo, giunti a questo punto, scrivere i valori rp, Cppqs 29 5.2 Caso variabile Normale 5.2.1 Presentazione del problema Sia X una variabile aleatoria distribuita come una normale con media µ sconosciuta e varianza σ2 “ 1. Si vuole testare l’ipotesi che µ sia minore di uno specifico valore µ1 . Anche questo tipo di problema, come visto nel caso di una distribuzione bernouliana, trova applicazione nel campo del controllo qualità. Supponiamo per esempio di avere un lotto composto da un grande numero di prodotti i quali devono essere sottoposti ad una ispezione prima di essere immessi nel mercato. Ipotizziamo di effettuare delle osservazioni su una caratteristica del prodotto (altezza, peso, durezza del materiale, ecc...). Il valore di tali osservazioni è distribuito come una normale con media sconosciuta e varianza conosciuta. Si considera migliore il prodotto man mano che µ diminuisce. In generale si può trovare un valore µ1 per cui preferiamo accettare il lotto per µ ă µ1 e rifiutarlo per µ ą µ1 . Per µ “ µ1 è indifferente la scelta da prendere. Come visto in precedenza, nel dominio µ ă µ1 , la preferenza per l’accettazione cresce al diminuire del vero valore µ, mentre nel dominio µ ą µ1 , la preferenza per il rifiuto cresce all’aumentare di µ. E’ possibile trovare due valori µ0 e µ1 con µ0 ă µ ă µ1 tali che l’accettazione del lotto se µ ě µ1 e il rifiuto quando µ ď µ0 sono considerati errori di rilevante importanza. Per valori compresi tra µ0 e µ1 ci troviamo nella regione di indifferenza. Si decide un determinato valore α per cui la probabilità di un errore di prima specie deve essere minore di esso, e un determinato valore β per cui la probabilità di un errore di seconda specie deve essere minore di quest’ultimo. 5.2.2 Implementazione del SPRT caso normale Un test che soddisfa le condizioni prima esposte è il sequential probability ratio test di potenza pα, βq che testa le ipotesi: # H0 : µ “ µ0 (5.14) H1 : µ “ µ1 Il test è definito così: sia xi l’ispezione i-esima, la densità di probabilità del campione xi è: pxi ´ µq2 ´ 1 2 f pxi , µq “ ? e (5.15) 2π 30 Quando µ “ µ0 e µ “ µ1 abbiamo: ´ 1 f pxi , µ0 q “ ? e 2π pxi ´ µ0 q2 pxi ´ µ1 q2 ´ 1 2 2 e f pxi , µ1 q “ ? e 2π (5.16) quindi passando al logaritmo del rapporto delle funzioni di densità congiunta (27) otteniamo: f pxi , u1 q 1 “ pu1 ´ u0 qxi ` pu20 ´ u21 q zi “ log 2 f pxi , uo q e m ÿ m p1m “ z1 ` z2 ` .... ` zm “ pu1 ´ u0 q xi ` pu20 ´ u21 q (5.17) log p0m 2 i“1 Il Sequential Probability Ratio Test di ipotesi semplici procede in questo modo: m ř 1´β m ñ rifiuta l’ipotesi nulla, se xi ` pu20 ´ u21 q ě 2 α i“1 m ř m β xi ` pu20 ´ u21 q ď se ñ accetta l’ipotesi nulla, 2 1´α i“1 altrimenti continua a campionare. 5.2.3 Considerazioni sul numero di iterate necessarie al termine del test Supponiamo di avere un vettore y “ py1 , y2 , ..., yn q generato da variabili aleatorie Yi „ Normp1, 1q indipendenti e identicamente distribuite. Si vuole testare in modo sequenziale H0 : µ “ 0 contro H1 : µ “ 1. Naturalmente il test, dopo aver esaminato un numero m non fissato di campioni, dovrà condurre al rifiuto dell’ipotesi nulla con una probabilità di errore di prima specie pari ad α. Il numero m di campioni estratti, necessari al test per rifiutare l’ipotesi nulla, è una variabile aleatoria. Essa dipenderà sia dalle dalle relizzazioni del campione che dagli errori di prima e seconda specie. In primo luogo si vuole determinare la distribuzione dei tempi necessari al test per rifiutare l’ipotesi nulla con α e β fissati: per farlo iteriamo il test un numero sufficientemente grande di volte in modo da ottenere un vettore contenente tutti i tempi di rifiuto. Fissando α “ 10´4 e β “ 5 ¨ 10´4 l’istogramma di questo vettore (quello rappresentato in Figura 4) dà alcune importanti informazioni sulla distibizione del numero di iterate necessarie al test per prendere una decisione. Fissando β, al diminuire di α la distribuzione delle frequenze assolute aumenta in media. man mano aumenta di valore e inoltre presenta una sempre più forte densità nella zona vicino allo zero. Nel caso in cui, invece diminuissimo solo β la distribuzione rimane invariata o si sposta solo leggermente verso destra. 31 600 Frequenza assoluta 500 400 300 200 100 0 0 10 20 30 40 50 60 70 80 90 Dimensione del campione Figura 5.2: Distribuzione dei campioni(tempi) necessari al rifiuto dell’ipotesi nulla in un SPRT sulla normale H0 : µ “ 0 contro H1 : µ “ 1 dove H1 è vera. 5.2.4 Confronto SPRT con test di Neyman Pearson caso Normale Il test di Neyman-Pearson di ampiezza α ha, come visto nel Lemma 9, la regione Lpθ0 ; px1 , x2 , .., xn qq ď ku dove la critica così costituita: C “ tpx1 , x2 , .., xn q | Lpθ1 ; px1 , x2 , .., xn qq funzione di verosimiglianza nel caso della normale è Lpθq “ ˆ 1 ? 2π ˙n e ´ n ř i“1 pxi ´ µq2 2 Prendendo sempre in considerazione il logaritmo del rapporto delle funzioni di verosimiglianza si vuole giungere, attraverso qualche passaggio, alla determinazione di k. n ÿ Lpθ0 ; px1 , x2 , .., xn qq n logp q “ pu0 ´ u1 q xi ` pu21 ´ u20 q ď logpkq 2 Lpθ1 ; px1 , x2 , .., xn qq i“1 n logpkq ´ pµ21 ´ µ20 q 2 xi ď “ k1 µ ´ µ 0 1 i“1 n ÿ 32 Ora, poiché la Ppx P Cq “ α, scriviamo: Pp n ř i“1 Xi ď k1 q “ α. In aggiunta conosciamo la distribuzione della sommatoria di n variabili aleatorie normali X1 pµ1 , σ21 q, X2 pµ2 , σ22 q, ..., Xn pµn , σ2n q “ Xpµ1 ` µ2 ` ... ` µn , σ21 ` σ22 ` ... ` σ2n q dunque per la somma di n variabili aleatorie normali con varianza unitaria, indin ř pendenti e identicamente distribuite abbiamo che Xi „ Xpnµ, nq. Dobbiamo ora esplicitare la funzione di ripartizione: i“1 n ÿ k1 ´ µ0 1 qq “ α Pp Xi ď k1 ; µ “ µ0 q “ p1 ` er f p ? 2 2n i“1 k1 “ ? 2n er f ´1 p2α ´ 1q ` µ0 n k1 pµ0 ´µ1 q` pµ20 ´µ21 q 2 k“e Vogliamo scrivere anche il valore β -errore di seconda specie- per fare dei confronti con il sequential test. Ricordando che PpX P C c ; µ “ µ1 q “ β si trova: n ÿ 1 k 1 ´ µ1 qqs “ β Pp Xi ą k1 ; µ “ µ1 q “ 1 ´ r p1 ` er f p ? 2 2n i“1 Si considera ora il test trattato nella sezione precedente con errori di prima e seconda specie equivalenti rispettivamente a 10´4 e 5 ¨ 10´4 . Si studia la distribuzione del numero di iterate necessarie al termine del SPRT per poi confrontarlo con il test di Neyman-Pearson. Come visto prima il numero di campioni necessari al Sequential Test è una variabile aleatoria; in questo specifico caso il numero di campioni necessari al termine del test è, con probabilità maggiore del 77% è minore di 30. Determinando un test di Neyman Pearson con errori dello stesso valore, il numero di campioni necessari è 51. Il risparmio è notevole (il 50% in media). Si può verificare che al diminuire degli errori il risparmio in percentuale del numero di campioni esaminati cresce, anche se molto lentamente. Per una dimostrazione teorica dei risultati qui riportati si consulti il testo di Abraham Wald-Sequential Analysis-Dover Publications (2004) (Dover Phoenix Editions) [1] nella sezione p3.6q. 33 Capitolo 6 Applicazioni del SPRT al moto browniano 6.1 Alcuni richiami sul moto browniano Si introduce ora la definizione di moto browniano al fine di effettuare alcuni test su di esso. Definizione 25. Un processo stocastico tXptq, t P T u è una famiglia parametrizzata di variabili aleatorie definite su uno spazio campionario Ω a valori in uno spazio misurabile S detto spazio degli stati del processo. Se T P N si parla di processo a tempo discreto, se, invece, T P R si parla di processo a tempo continuo. Per quanto riguarda lo spazio degli stati: quando S P N il processo è discreto, mentre quando S P R il processo si dice continuo. Definizione 26. Un processo stocastico tXptq, t ą 0u è detto Moto Browniano se: 1. Xp0q “ 0 2. Xptq ha incrementi indipendenti e stazionari 3. Xptq „ Normp0, σ2 tq Dalla proprietà di stazionarietà e dal fatto che conosciamo la distribuzione di Xptq si può dedurre che: Xptq ´ Xpsq „ Normp0, σ2 pt ´ sqq Definizione 27. Un moto browniano si dice moto browniano standard se σ2 “ 1 e si indica con Bptq, t ě 0. Definizione 28. Si ha un moto browniano con drift quando Xptq „ Normpµt, σ2 tq che corrisponde a dire: Xptq “ µt ` σBptq 34 6.2 Presentazione del problema In questa sezione si pone il problema di riuscire a capire se il moto browniano Xptq „ Normpµt, σ2 tq con σ2 conosciuto ha un drift che non eccede un determinato valore µ1 . Il problema sorge da alcune applicazioni, ad esempio quelle finanziarie, in cui abbiamo dati che ci vengono forniti in tempo reale e continuamente (o quasi) e le osservazioni vengono fatte ad intervalli di tempo anche irregolari. All’aumentare del drift si preferisce sempre di più comprare il titolo che si sta ispezionando. Generalmente è possibile trovare un µ1 per il quale se µ ą µ1 si preferisce entrare e se µ ă µ1 si preferisce uscire o non entrare. E’ necessario ricordare che gli incrementi tra due tempi t e s si distribuiscono come una normale con media µpt ´ sq e varianza σ2 pt ´ sq. Denotiamo con Xi gli incrementi stazionari e indipendenti che hanno come intertempo h “ 1. La sommatoria di questi intertempi con i che va da 1 a n con i P N` genera un moto browniano Xpnq che in questo caso è a tempo discreto. L’affermazione è vera per i seguenti motivi: siano Xi con i “ 1, 2, .., n variabili aleatorie normali di media µi e varianza σ2i allora n ÿ i“1 Xi „ Normpµ1 ` µ2 ` .. ` µn , σ21 ` σ22 ` .. ` σ2n q se le Xi sono indipendenti e identicamente distrubuite si ottiene: n ÿ i“1 Xi „ Normpnµ, nσ2 q Il problema in questione si riconduce semplicemente al caso del SPRT nel caso normale. Difatti, se si campiona ad intervalli di tempo regolari, per testare delle ipotesi sul drift del moto browniano basta effettuare un sequential probability ratio test con ipotesi sulla media degli incrementi. Se, però, gli intervalli sono irregolari, la funzione di densità congiunta di f ppx1 , x2 , .., xn q; µq dipenderà anche dagli intertempi h tra una osservazione e l’altra. 6.3 Individuazione del cambiamento nel drift In questa sezione si pone il problema di riuscire a capire quando un moto browniano Xptq „ Normpµt, σ2 tq con σ2 noto cambia il suo drift. Supponiamo che un moto browniano che si sta realizzando in tempo reale cambi il suo drift ad un certo istante θ. Quindi @t ă θ Xptq „ Normp0, tq, mentre @t ě θ Xpt ´ θq „ Normpµpt ´ θq, t ´ θq. Si vuole trovare un metodo che riesca repentinamente a 35 SPRT sul moto browniano 6 5 B(t) 4 Accetto in questo punto H 0 3 2 1 0 0 5 10 15 numero di iterate 20 25 Figura 6.1: test di ipotesi semplici sul drift del moto browniano che conduce all’accettazione dell’ipotesi nulla µ “ 0. capire quando ciò accade. Ovviamente ci sono numerosi esempi applicativi in cui si ha la necessità di affrontare questo tipo di problema. Nel libro "Sequential Analysis: Hypothesis Testing and Changepoint Detection" scritto da Alexander Tartakovsky, Igor Nikiforov e Michele Basseville[5]. In ambito medico, supponiamo di testare l’effetto di una determinata cura o un determinato farmaco su uno specifico valore del paziente (ad esempio il livello di tossicità). Poniamo che il livello di tossicità si comporti come un moto browniano standard Xptq con Xp0q “ φ dove φ rappresenta il livello di tossicità del paziente prima della somministrazione. Possiamo ipotizzare che, prima di immettere il prodotto sul mercato, si operino, per un tempo fissato preliminarmente, controlli sulle variazioni del valore conseguenti alla somministrazione del farmaco/cura. Questa fase viene chiamata "fase di sperimentazione", il cui scopo è fornire una valutazione sulla sicurezza e/o tollerabilità del medicinale. E’ maggiormente desiderabile un farmaco che abbassi o lasci invariato il valore specifico testato, mentre il farmaco non dovrà passare alle fasi successive di sperimentazione, se la somministrazione comporta un evidente innalzamento del valore. Il test si riconduce sempre ad un SPRT caso normale -con varianza conosciuta, il quale testa la media degli incrementi (con campionamenti effettuati ad intervalli regolari). Anche in questo caso è possibile trovare due valori µ0 e µ1 , tali che µ0 ă µ1 e che: per valori superiori a µ1 si commette un errore rilevante permettendogli di accedere alle fasi successive, per valori inferiori a µ0 l’errore nello stop 36 SPRT sul moto browniano 25 Rifiuto in questo punto H 0 20 15 Accetto B(t) 10 Accetto Accetto 5 0 Accetto -5 -10 0 50 100 150 200 250 numero di iterate Figura 6.2: applicazione del metodo sequenziale di ipotesi H0 : µ “ 0 contro H1 : µ “ 0.5 nel caso in cui dal tempo 0 sino a τ si hanno incrementi Xi „ Normp0, 1q, mentre per tempi maggiori di τ gli incrementi Xi „ Normp0.5, 1q. della produzione sarebbe importante. Bisogna anche considerare un altro aspetto: la somministrazione del farmaco potrebbe non presentare effetti collaterali immediati e dunque è necessario monitorare il paziente per un periodo medio-lungo. Il test terminerà se ci troviamo nella zona di rifiuto, mentre continuerà a campionare se ci troviamo nella regione di indifferenza o in quella di accettazione. Quando il test conduce all’accettazione dell’ipotesi nulla non fermiamo l’ispezione poiché l’effetto collaterale potrebbe presentarsi in qualsiasi momento della fase di sperimentazione. Una risposta al nostro problema è un sequential probability ratio test caso normale di ipotesi semplici H0 : µ “ µ0 contro H1 : µ “ µ1 con un numero di iterate massimo (corrispondente al tempo di sperimentazione) che, però riparte da capo ogni volta che viene accettata l’ipotesi nulla. In sostanza iteriamo il SPRT caso normale sino al rifiuto dell’ipotesi nulla. In Figura 6.3 possiamo vedere il comportamento del test su una particolare realizzazione di un moto browniano. 37 6.4 Tempi di individuazione Facciamo ora un esempio numerico per verificare l’efficacia del metodo. Il test che prendiamo in esame è quello di ipotesi: # H0 : µ “ 0 (6.1) H1 : µ “ 0.5 con α “ 0.010 e β “ 0.001. Questo test, come in precedenza spiegato, continuerà ad iterarsi sino al rifiuto dell’ipotesi nulla o all’arrivo del numero massimo di iterate prefissato. Il moto browniano testato sappiamo a priori che è un moto browniano standard sino τ “ 200 e poi acquista un drift µ “ 0.5. Dopo averlo implementato in macchina, si effettua questo stesso test 10000 volte -ovviamente le realizzazioni del browniano sono differenti ogni volta- al fine di fare qualche osservazione sui tempi -o numero di iterate- che intercorrono tra l’istante τ “ 200 e l’istante t in cui viene rifiutata l’ipotesi nulla. Il grafico riportato in Figura 6.4 ci mostra le frquenze assolute degli t ´ τ. I valori maggiori sono quelli che occorrono nel 98% dei casi. Questo indica che la probabilità di falso allarme rimane molto limitata. Osserviamo inoltre che il 95% dei dati si trovano nell’intervallo tra 0 e 50. Questo dimostra la bontà del test. Tuttavia, possono sorgere alcune problematiche. Le frequenze assolute per le classi appartenenti ai valori negativi di t ´ τ sono costanti; i valori negativi rappresentano gli errori del test quando accetta l’ipotesi che ci sia un drift quando in realtà non c’è. Se effettuiamo la stessa procedura ma prendiamo un τ (punto di cambiamento) molto più grande le frequenze per tutti i valori negativi rimarranno costanti, ma in questo caso ci saranno molte più classi negative. La probabilità di effettuare una scelta sbagliata aumenta man mano che allontaniamo l’istante τ dall’istante in cui inizia il procedimento. Questo problema deriva essenzialmente dall’errore di prima specie: ogniqualvolta il test comincia vi è una probabilità „ α di rifutare l’ipotesi nulla commettendo un errore; così iterando il test questa probabilità si moltiplica per il numero di test necessari ad arrivare sino a τ. Potrebbe essere utile utilizzare un procedimento nel quale il valore α è molto più piccolo. Questo cambiamento genererebbe sicuramente un forte abbassamento nella probabilità di un errore nel rifiutare, ma comporterebbe un innalzamento della media degli t ´ τ. Possiamo concludere che la suddetta procedura può essere utilizzata soltanto quando il periodo di monitoraggio è sufficientemente breve. 38 1500 1000 500 0 -200 -150 -100 -50 0 50 100 150 200 250 Figura 6.3: Istogramma degli t ´ τ ottenuti mediante l’iterazione per 10000 volte del metodo sequenziale di ipotesi H0 : µ “ 0 contro H1 : µ “ 0.5 con α “ β “ 0.01 nel caso in cui dal tempo 0 sino a τ si hanno incrementi Xi „ Normp0, 1q, mentre per tempi maggiori di τ gli incrementi Xi „ Normp0.5, 1q. 6.5 Analisi del metodo Il metodo da noi utilizzato presenta un problema principale: la sottostima del drift dopo l’istante τ di cambiamento. In questa sezione si vuole trovare un metodo per attenuare questa problematica. Verrà proposto anche un altro approfondimento relativo alla dipendenza dei parametri dati nel test. 6.5.1 Sottostima Per capire meglio la problematica osserviamo la Figura 6.5 e la Figura 6.6. Nel primo caso abbiamo una situazione ideale: ovvero l’ipotesi nulla è stata accettata appena prima che il drift mutasse e questo ha reso possibile un’individuazione veloce. Nel secondo caso l’ultima accettazione avviene molto prima e dunque questo rallenta la decisione della procedura. Questo avviene perché, dopo aver accettato l’ipotesi nulla, nella sommatoria degli incrementi x1 ` x2 ` ... ` xθ ` ...xn quelli che hanno media µ “ 0 avranno un certo peso, il quale condurrà ad una sottostima del drift e si avrà bisogno di un maggior numero di campioni per ri39 SPRT sul moto browniano 50 Rifiuto in questo punto H 40 30 Accetto B(t) 20 Accetto 10 Accetto Accetto Accetto Accetto 0 -10 -20 0 50 100 150 200 250 300 350 400 450 numero di iterate Figura 6.4: Illustrazione delle regioni di rifiuto,accettazione e indifferenza applicate al metodo sequenziale di ipotesi H0 : µ “ 0 contro H1 : µ “ 0.5 nel caso in cui dal tempo 0 sino a τ si hanno incrementi Xi „ Normp0, 1q, mentre per tempi maggiori di τ gli incrementi Xi „ Normp0.5, 1q. SPRT sul moto browniano 20 Rifiuto in questo punto H 0 15 10 5 B(t) 0 -5 -10 -15 -20 -25 0 50 100 150 200 250 numero di iterate Figura 6.5: Illustrazione analoga alla Figura 6.5. In questo caso le realizzazioni hanno condotto ad un rifiuto meno repentino. 40 fiutare l’ipotesi nulla rispetto ad un test che prende in considerazione řn solo i dati da τ in poi. Una possibile soluzione si ottiene sostituendo sm “ i“1 xi con una media pesata degli incrementi xi con pesi pi opportunamente scelti in modo che più il dato è vicino al tempo presente e maggiore sarà il suo peso moltiplicato per le iterate. Dunque x1 p1 ` x2 p2 ` ...xn pn sm “ n p1 ` p2 ` ... ` pn dove p1 ă p2 ă .. ă pn . Tuttavia non è facile creare un buon metodo che scelga i pesi in modo opportuno. Decidiamo allora di considerare solo gli ultimi k incrementi e ogni volta che si effettua uno step (quindi una nuova osservazione) si cancella il più "vecchio". Il test procederà, dunque, in questo modo: consideriamo il campione x1 , x2 , ..xn dopo aver effettuano k iterate senza effettuare una decisione l’agoritmo considererà solamente gli ultimi k incrementi. Se ci trovassimo allo step n´esimo senza ancora nella regione di indifferenza considereremmo solamente xn´k , ....., xn In questo modo il nostro test avrà, dopo aver effettuato le prime k iterate senza prendere una decisione, una numerosità campionaria prefissata. Osservazione sulla scelta di k Per scegliere k bisogna pensare a come attenuare il problema di sottostima. Una possibilità per la scelta di k è questa: si valuta Eµ pnq del sequential test Q˚ di ipotesi # H0 : µ “ µ0 con µ1 ą µ0 H1 : µ “ µ1 Eµ pnq rappresenta il valore atteso del numero di campioni necessari al test Q˚ per effetuare una decisione al variare del vero valore µ. A questo punto imponiamo Eµ1 pnq “ k in modo che, con alte probabilità, quando rifiutiamo l’ipotesi nulla le ultime k ispezioni appartengano ad incrementi generati da una variabile aleatoria normale con media maggiore o uguale a µ1 . 6.5.2 Dipendenza dai parametri Sia µ1 il drift del moto prima del cambiamento, si vuole costruire un test il quale: fino a quando µ ´ µ1 ă ε (con ε sufficientemente piccolo) con alte probabilità si 41 2000 1800 1800 1600 1600 1600 1400 1400 1400 1200 1000 800 Frequenza assoluta 2000 1800 Frequenza assoluta Frequenza assoluta 2000 1200 1000 800 1200 1000 800 600 600 600 400 400 400 200 200 200 0 0 0 5 10 15 20 25 30 Dimensione del campione 35 40 45 50 0 5 10 15 20 25 30 Dimensione del campione 35 40 45 50 0 0 5 10 15 20 25 30 35 40 45 50 Dimensione del campione Figura 6.6: I tre grafici illustrano la frequenza assoluta dei tempi di rifiuto di un SPRT di ipotesi H0 : µ “ µ0 contro H1 : µ “ µ1 con α “ β “ 0.01 dove il vero valore è µ “ 1. Il primo grafico a sinistra è il caso µ1 “ 1, il secondo è µ1 “ 0.5 ed il terzo µ1 “ 0.2 rimane nella regione di indifferenza e solo quando si altera il drift, il test condurrà ad una decisione. Il processo prima presentato riscontrava il problema di intersecare troppe volte la retta che delimitava la regione di accettazione. Ogni volta che avveniva un’intersezione le regioni venivano modificate come illustrato in Figura 6.5. Questo problema potrebbe essere attenuato se riuscissimo a non accettare l’ipotesi nulla, anche se vera, troppo facilmente. In questo modo sarà possibile iterare meno volte il Sequential test e non accumulare errori. Il nostro proposito è dunque di ragionare sulle variazioni delle regioni di accettazione, rifiuto e indifferenza. Sia S il nostro SPRT caso normale di ipotesi semplici µ0 e µ1 tali che 1 ą µ1 ą µ0 e sia µ1 ´ µ0 “ ϵ. Al diminuire di ϵ la regione di indifferenza si amplia e inoltre i coefficienti angolari delle due rette che delimitano la regione di indifferenza diminuiscono. Effettuiamo degli esempi, che vengono riportati in Figura 6.7 al fine di verificare la precedente osservazione. Usiamo ora la procedura con µ1 minore di quello utilizzato prima per testare un moto (che noi conosciamo) equivalente a quello creato nella sezione precedente, nonostante si ha una preferenza per l’accettazione per tutti i valori ă 0 ed una preferenza per il rifiuto per valori ą 0.5. Prendiamo un caso particolare: # H0 : µ “ 0 (6.2) H1 : µ “ 0.2 con α “ 0.01 e β “ 0.01. Il test, condurrà al rifiuto con probabilità P “ p1 ´ αq se il vero valore è 0.2, mentre per valori più grandi di µ la probabilità aumenterà (diminuirà αpθq). Vediamo in Figura come le regioni siano cambiate. Il cambiamento applicato si traduce anche in tempi maggiori necessari al test per fornire una scelta e di conseguenza in un cambiamento sostanziale nella di42 40 30 20 X(t) 10 0 -10 -20 -30 -40 0 50 100 150 200 250 300 numero di iterate Figura 6.7: applicazione del metodo sequenziale di ipotesi H0 : µ “ 0 contro H1 : µ “ 0.2 nel caso in cui dal tempo 0 sino a τ si hanno incrementi Xi „ Normp0, 1q, mentre per tempi maggiori di τ gli incrementi Xi „ Normp0.5, 1q. In blu e arancione si possono notate le linee che delimitano rispettivamente le regioni di rifiuto e accettazione. stribuzione dei tempi (numero di iterate) τ ´ t. La Figura 6.9 dimostra come sia cambiato il grafico delle frequenze assolute. Questo metodo ha il vantaggio di non accumulare errori all’aumentare di τ, ma man mano che l’istante di cambiamento è più lontano rispetto all’istante di partenza dell’ispezione le frequenze di t ´ τ si sposteranno verso destra: cioè per il problema di sottostima il test sarà più lento ad effettuare una decisione. Anche in questo caso si può utilizzare il metodo prensentato nella sezione 6.10.1. 43 1000 900 800 frequenza assoluta 700 600 500 400 300 200 100 0 -200 -150 -100 -50 0 50 100 150 200 250 300 tempo t-theta Figura 6.8: Istogramma rappresentante i tempi t ´ τ trovati dall’iterazione per 10000 volte del metodo sequenziale di ipotesi H0 : µ “ 0 contro H1 : µ “ 0.2 nel caso in cui dal tempo 0 sino a τ si hanno incrementi Xi „ Normp0, 1q, mentre per tempi maggiori di τ gli incrementi Xi „ Normp0.5, 1q. In blu e arancione si possono notate le linee che delimitano rispettivamente le regioni di rifiuto e accettazione. 44 Capitolo 7 Conclusioni In questa tesi sono stati esposti i pilastri fondamentali della teoria statistica sequenziale usando come riferimento principale il libro "Sequential Analysis" di Abraham Wald. Si è fatto uso di questa teoria per alcune applicazioni. L’applicazione del Sequential Probability Ratio Test al caso di una variabile aleatoria bernoulliana riveste un ruolo fondamentale nei problemi attinenti il controllo qualità. L’applicazione al caso normale risulta, invece, essenziale per applicare il metodo sequenziale ad un moto browniano. La questione di fondo che anima questa tesi è, difatti, legata ad un problema di detection: cioè l’individuazione del cambiamento del drift di un moto browniano. Applicando la teoria dei test sequenziali al moto browniano abbiamo ottenuto dei metodi per la detection che hanno tempi di individuazione inversamente proporzionali alle probabilità di individuare un falso allarme. Sono, dunque, molto soddisfacenti nel caso in cui l’obbiettivo sia solo l’individuazione corretta e non il tempo necessario ad essa. L’utilizzo di α molto piccolo è una soluzione efficiente in questo caso particolare. Gli algoritmi presentati vogliono essere una rappresentazione parziale dei possibili metodi di soluzione, ma forniscono comunque alcuni importanti spunti per l’approccio al problema in questione. Nell’ultima parte della tesi si è accennato ad alcune soluzioni che potrebbero rendere l’algoritmo più efficiente. Tuttavia bisognerebbe continuare lo studio delle dipendenze dei parametri dati in un SPRT al fine di capire come essi influiscano sulla suddivisione dello spazio campionario nelle tre regioni (accettazione, rifiuto e indifferenza). Possibili sviluppi futuri potranno riguardare anche la scelta di quali sono le ispezioni che si devono considerare al fine di ottimizzare il metodo. Tale scelta dipenderà non solo da considerazioni preliminari, ma anche dalle realizzazioni del campione. 45 Capitolo 8 Appendice 8.1 Codice matlab: Sequential Probability Ratio Test caso Bernoulli Viene riportato qui il codice implementato su Matlab dell’algoritmo presentato nella sezione 5.2. function[n,h]=SPRT_bernoulli(p0,p1,x,alpha,beta) % Input:p0 é il valore da testare come ipotesi nulla % p1 è il valore da testare come ipotesi alternativa % x è il lotto di prodotti molto grande % alpha è l’errore di prima specie % beta è l’errore di seconda specie % % Output: n vale 0 se accetto l’ipotesi nulla e 1 se la rifiuto % h è il numero delle iterate n=99; h=1; c=0; for iter=1:1000 if x(h)==1 d=log(p1/p0); end if x(h)==0 d=log((1-p1)/(1-p0)); end c=c+d; if c >log( (1-beta)/alpha ) n=1; 46 break end if c < log(beta/(1-alpha)) n=0; break end end L’output sarà n “ 1 se rifiuto l’ipotesi nulla, altrimenti n “ 0. 8.2 Codice matlab: sequential probability ratio test caso Normale Riportiamo qui il codice implementato su Matlab dell’algoritmo presentato nella sezione 6.2. function[n,h]=SPRT_normale(u0,u1,y,alpha,beta) % Input: u0 è il valore da testare come ipotesi nulla % u1 è il valore da testare come ipotesi alternativa % y è il lotto di prodotti distribuito come una normale % alpha è l’errore di prima specie % beta èl’errore di seconda specie % Output: n vale 0 se accetto l’ipotesi nulla e 1 se la rifiuto % h è il numero delle iterate n=99; h=1; c=0; for iter=1:1000 d= ((u1-u0)*y(h))+ 0.5*((u0)2 -(u1)2 ) c=c+d; if c>log( (1-beta)/alpha ) n=1; break end if c < log(beta/(1-alpha)) n=0; break end h=h+1; end end 47 8.3 Neyman test Qui di seguito il codice implementato in Matlab dell’algoritmo visto in sezione 6.4 function [ beta ] = Neymantest( u0, u1, alpha, n ) x=randn(n,1)+1; L=0; k0= sqrt(2*n)* erfinv(2*alpha-1) + n*u0; k1= exp(k0*(u0-u1)+ (n/2)*(u12 -u02 )); for i=1:n L=L+(u0-u1)*x(i)+ (1/2)*(u12 -u02 ); end if L>k1 H=0; fprintf(Accetto ipotesi nulla) else H=1; fprintf(Rifiuto ipotesi nulla ipotesi nulla) end beta= 1- ( (1/2)*(1+erf( (k0-u1)/(sqrt(2*n)) ))); end 8.4 Test sul Moto Browniano Riportiamo ora un codice che racchiude tutti i test e tutte le modificazioni viste nella sezione 7. function[h,n,t,aM,rM,dM]=SprtBrownianoCompost(k,u0,u1,alpha, beta,drift,tau) aM=linspace(0,10,101); rM=linspace(0,10,101); dM=linspace(0,10,101); diff=0; zz=linspace(-500,100,1001); te=linspace(0,10,10001); cc=1; for iter=1:10000 s=0; 48 n=99; q=0; h=1; c=0; d=linspace(0,10,10001); cont=linspace(0,10,501); B=linspace(0,10,10001); k=1; y=randn(10001,1); j=tau; for i=1:10000-tau j=j+1; y(j)=y(j)+drift; end for iter=1:10000 k=k+1; B(k)= B(k-1)+y(k); end for iter= 1:10000 for i = 1:length(d) h=h+1; dM(h)=dM(h-1)+y(h-1); if h>k dM(h)=dM(h)-y(h-k); end aM(h)=(1/(u1-u0))*log(beta/(1-alpha))+(i)*((u0+u1)/2)+diff; if i>k aM(h)=(1/(u1-u0))*log(beta/(1-alpha))+(k)*((u0+u1)/2)+diff; end rM(h)=(1/(u1-u0))*log((1-beta)/alpha)+(i)*((u0+u1)/2)+diff; if i>k rM(h)=(1/(u1-u0))*log((1-beta)/alpha)+(k)*((u0+u1)/2)+diff; end if dM(h)>=rM(h); n=1; break end if dM(h)<=aM(h); n=0; q=q+1; break 49 end end s=s+1; cont(s)=h; c=0; diff=B(h) te(cc)=h-tau; if n==1 break end end cc=cc+1; end hist(te,50) end 50 Bibliografia [1] Abraham Wald, (1948). Sequential Analysis. Dover Phoenix Editions. [2] David Siegmund, (1985). Sequential Analysis. Test and confidence. Springer series in statistics. [3] Abraham Wald, (1951). Sequential Test of Statistical Hypothesis. The Annals of Mathematical Statistics. [4] R.L. Berger,(2001). Statistical Inference, Duxbury Press [5] Alexander Tartakovsky, Igor Nikiforov, Michèle Basseville ,(2015).Sequential Analysis: Hypothesis Testing and Changepoint Detection, Chapman and Hall Book [6] Shewhart Walter, (1931). Economic Control of Quality Manufactured Product. American Society for Quality (ASQ) 51