Processi stocastici

Capitolo 1

Finora abbiamo sempre assunto i nostri dati campionari i.i.d.. In particolare, questa assunzione implica (in ordine crescente di importanza): indipendenza dell’errore, omoschedasticità e normalità. Le ultime due, per il teorema del limite centrale, possono essere considerate valide al crescere del numero di dati; la prima l’abbiamo sempre assunta come vera. Quest’assunzione, d’ora in poi, verrà messa in discussione.

Dipendenza temporale nei dati

Supponiamo adesso di avere dei dati che abbiano una componente di dipendenza nei loro errori: in particolare facciamo riferimento a una dipendenza temporale formalizzata come una serie di osservazioni che si succedono ordinatamente nel tempo.

Lo studio di fenomeni in cui i dati osservati sono in successione ordinata richiede formulazioni probabilistiche differenti da quelle viste finora, che spesso falliscono proprio per la dipendenza dal tempo tra osservazioni successive. Distinguiamo pertanto due nuove categorie:

• La prima è legata ai dati e tipicamente dà luogo al cosiddetto survival study, in cui i dati aleatori riguardano il tempo di uno o più eventi puntuali o il numero di arrivi in un intervallo temporale.
• Il secondo si riferisce ad osservazioni di variabili casuali che si evolvono nel tempo, dando luogo alle cosiddette time series.

Processi di arrivo

Approfondiamo la prima categoria, ovvero i processi di arrivo. Ammettiamo che diversi arrivi possano verificarsi nel tempo uno dopo l'altro. In questo caso si tratta di un fenomeno che non cessa di esistere dopo il primo arrivo, come un fallimento paralizzante di un elemento, ma un processo che può vedere arrivi multipli.

Sia N(t) la variabile conteggio che esprime il numero di arrivi nell’intervallo (0, t]; questa è non-decrescente con t, così come lo è il suo valore atteso M(t). Se sono possibili arrivi simultanei, la funzione non sarà continua, tuttavia rimane continua a destra. La sua derivata m(t) si definisce arrival rate (o rate of occurence of failure ROCOF) ed è il tasso istantaneo di cambio nel valore atteso degli arrivi.

• Se il k-esimo arrivo arriva dopo il tempo t (T > t), allora al tempo t si ha N(t) < k: ciò significa che è possibile ragionare in due modi: per fissato t ci si può chiedere quanti arrivi ci sono stati oppure fissare k e chiederci questo k-esimo arrivo dove si colloca.
• Se T ≤ t allora N(t) ≥ k.
• Se t < T < t, allora N(t₁) < k ≤ N(t₂).

In questo caso invece se facciamo le differenze temporali, o prendiamo 2 valori del numero di arrivi, se il k-esimo arrivo si verifica nell’istante tra t₁ e t₂ significa che il nostro k sarà minore o uguale N(t₂) ma maggiore di N(t₁).

C’è una corrispondenza biunivoca tra tempi e numero di arrivi sulla distribuzione congiunta: questo implica che i gradi di libertà congiunti di T₁, T₂, …, T_n, per ogni n, sono in corrispondenza uno ad uno con quelli di N(T₁), N(T₂), …, N(T_n).

Processi stocastici e intensità

Un processo stocastico si dice definito con incrementi stazionari se risulta che dipende unicamente dal tempo di interarrivo s e non dalla sua posizione assoluta t, per ogni k. Un processo stocastico si dice definito con incrementi indipendenti se risulta (non gli eventi, gli incrementi!).

La funzione intensità i.f. è definita come probabilità che in un intervallo infinitesimo avvenga almeno un arrivo: lim = 0.

Sottolineiamo la differenza con il tasso di rischio h.r. h(t): questo misura la probabilità istantanea che l’unico possibile evento si verifichi in un intervallo infinitesimo, ammettendo che non si sia verificato prima. A differenza del tasso di rischio definito nel processo di Poisson, qui gli arrivi simultanei sono permessi.

In particolare è una funzione a gradini costante pari a 0 fin quando non arriva il primo evento all’istante t₁, poi costante a 1 fin quando non arriva il secondo evento all’istante t₂ e così via.

Attenzione perché la nostra variabile N è una funzione continua a destra; infatti, nel momento in cui si verifica un nuovo evento, la nostra funzione fa il salto quindi è discontinua a sinistra e continua a destra (come si vede dai segni della disequazione), oppure è costante.

Comunque dimostrabile che, ammettendo che gli arrivi simultanei siano impossibili (nello stesso intervallo), si dimostra che i(t) e m(t) sono equivalenti. Si osservino i seguenti passaggi: Δ, lim Δ → 0, … .

Dunque la i.f. misura la frequenza attesa degli arrivi. Sotto le seguenti condizioni si ottiene un processo di Poisson (P.P.) di densità di probabilità conosciuta:

• Il processo è a incrementi indipendenti.
• Non esistono eventi simultanei ed è definita una i.f. i(t)=dI(t)/dt con valore atteso M(t)=I(t) ed R(t)=exp{–I(t)}.

Se N(t) è definito in (a, b], si ha:

I processi dove può aversi più di un arrivo sono quelli in cui l’evento non è fatale e lascia l’affidabilità del sistema intatta e si parla di minimal repair. Possiamo campionare fintantoché non si raggiunge un certo numero di arrivi n, o per un dato periodo di osservazione. Supponiamo che sia T l’istante in cui si è raggiunto un numero n di arrivi. Allora T=T_n è una variabile casuale corrispondente all’istante dell’ultimo arrivo, mentre n non è fissato.

Densità di probabilità e processi di Poisson

La d.f. del primo T e secondo arrivo dato il primo T₂|T₁ sono:

Considerando la d.f. dell’ultimo arrivo f_n(t|t₁, t₂, …, t_n-1), la densità di probabilità congiunta f(t₁, t₂, …, t_n) è:

Si dimostra che la p.d. dell’n-esimo arrivo entro t risulta un P.P.

Possiamo ottenere la d.f. dell’n-esimo tempo di arrivo assumendo un processo di degradazione che segue la distribuzione di Weibull. Questa variabile casuale è chiamata Gamma generalizzata, e per n=1 è uguale ad una Weibull, mentre per α=1 ad una Gamma. Per n=α=1 si ottiene un processo di Poisson omogeneo. La d.f. e la p.d. P(n) dell’n-esimo arrivo al tempo t risultano:

Gli istanti di arrivo di un P.P., dato l’istante dell’ultimo arrivo T_n, sono n – 1 variabili casuali così distribuite:

Supponiamo adesso di campionare per un periodo di osservazione t. Allora adesso anche N(t) è una v.c. Effettuando i dovuti passaggi si trova che:

Lo studio di questi casi particolari è necessario perché possono aversi diverse regole d’arresto. Il processo osservato potrebbe andare avanti all’infinito: la regola d’arresto è il modo in cui fermo l’osservazione, ovvero se ad un certo tempo fissato oppure ad un n-esimo arrivo.

Vi è un modo per considerare anche processi di rinnovo. In questo caso ipotizziamo però che dopo ogni arrivo, vi sia un elemento disponibile pronto a subentrare, per un massimo di k elementi. Ciò significa che dopo ogni arrivo l’affidabilità ricomincia da capo, come al tempo 0.

In questo caso la vita utile del sistema è data dalla somma delle individuali T_k = X₁ + X₂ + … + X_k e vale:

Con questa calcoliamo la f. di densità congiunta, ovvero la probabilità condizionata che si abbia il primo arrivo a t₁, il secondo a t₂, l’n-esimo a t_n.

Con questa calcoliamo la f. di densità congiunta condizionata a n dato t.

X₁, X₂, …, X_k sono i tempi di interarrivo (X₂ = quanto abbiamo aspettato tra il primo e il secondo arrivo).

La d.f. di T_k è la convoluzione delle X cioè l’intervallo di tempo prima che si verifichi il k-esimo arrivo che è dato dalla somma delle variabili aleatorie di interarrivo che ora sono indipendenti (in quanto relative a processi diversi). Risulta E(T_k)= E(X₁)+ E(X₂)+…+ E(X_k), var(T_k)= var(X₁)+ var(X₂)+…+ var(X_k).

Quando k diventa abbastanza grande, è possibile applicare il teorema limite centrale. Se supponiamo poi che le X siano anche identiche con media μ e varianza σ², si ha E(T_k)=kμ, var(T_k)=kσ² e dunque possiamo scrivere:

\[\frac{k}{\sqrt{k}}\Phi(x) = T(k)\]

Con la prima si considera T(k) come v.c., mentre con la seconda la v.c. è k, ovvero dove naturalmente \( z - 1 \). Quadrando, si ottiene una disequazione si 2° grado in N(t), con soluzione:

\[2μ - 4μ^2 \leq 2μN(t) \]

Se vogliamo Prob{N(t)≥k}= Prob{T_k ≤t}, allora deve scegliersi l’estremo superiore, viceversa l’estremo superiore.

Inferenza su sistemi riparabili

Vediamo ora come fare inferenza su sistemi riparabili nel caso di processi di Poisson omogenei (HPP) e non (NHPP) secondo una legge esponenziale.

Nel caso di HPP si ha che MLE=m = λ (reciproco della media aritmetica) e che 2λt - 1, dal quale si stima facilmente l’intervallo di stima di λ.

Se ragioniamo fissando il tempo t, il numero di arrivi N è una v.c. distribuita come una Poisson con media e varianza pari a μ(t)=λt ed MLE=N/t. Naturalmente il tempo totale è dato dalla somma di tutti questi tempi ed essendo somme di arrivi governati da una variabile casuale esponenziale con tasso di arrivo costante, si distribuisce come una Gamma. Legando una Gamma generica alla gamma con λ=1 otteniamo λt = y= ½ z e la chi-quadrato è quella particolare Gamma che ha λ= 1/2 e ν= 2α (in questo caso α è il numero di eventi che si sono sommati). Il numero di arrivi è una variabile di Poisson, il valore atteso è pari a λt, pertanto l’MLE di λ è N/t e la distribuzione del campione per N>10 può approssimarsi ad una Gaussiana.

Nel caso di NHPP che segue una legge esponenziale (Weibull), un sistema riparabile osservato sino all’n-esimo arrivo ha una siffatta funzione di verosimiglianza (calcolata a partire dalla funzione di densità):

Se i(t_k) è una costante troviamo il processo di Poisson omogeneo, se è generica allora abbiamo il processo di Poisson non omogeneo. Per quest’ultimo caso dobbiamo assumere che la forma della distribuzione e l’unica quantità non nota siano i parametri di questa distribuzione.

Se consideriamo le t₁, t₂ come variabili casuali i parametri sono α e λ. Se questi parametri sono noti allora guardiamo la funzione come una funzione di densità, t₁, t₂ sono valori osservati mentre dobbiamo fare inferenza su α e λ.

Il problema delle verosimiglianze viene risolto eguagliando a zero e derivando rispetto ai parametri:

Per ottenere gli intervalli di confidenza per β e, di conseguenza, per α si noti che 2n λ / β quindi 2n λ / β < (intervallo di confidenza che stacca a sinistra e a destra).

Per quanto riguarda θ si fa riferimento a tavole della statistica fondamentale W = (λ/θ), ottenute per simulazione. Tuttavia gli intervalli di confidenza per θ sono spesso più ampi di quelli così stimati, la stima di α non è stabile e inoltre θ non ha una diretta interpretazione se non quando α=1. Di solito questo θ non viene stimato; del resto, si è più interessati al primo parametro perché dice se il sistema è stabile o se il processo si intensifica o meno.

A tal proposito un’altra importante stima è quella della funzione di intensità a t, che per un numero n di arrivi fissati può essere così stimato attraverso MLE:

Come vediamo, per n fissato, non dipende da λ; gli intervalli di confidenza sono ricavabili dai principali software statistici.

Sistemi riparabili e tempo di osservazione fisso

Passiamo al caso in cui si studia un sistema riparabile osservando per un tempo fissato t. Se α > 1 significa che l’intensità aumenta con il passare del tempo e la frequenza di arrivo va aumentando, ovvero l’intertempo è sempre minore; viceversa se è 0<α<1.

Si noti che la sommatoria nella stima di α si arresta a n-1 perché se al posto di t si mettesse t_k il rapporto fa 1 e, poiché il logaritmo di 1 fa zero, l’ultimo addendo si annulla. Segnaliamo inoltre che si distribuisce come una chi-quadrato con 2(n-1) gradi di libertà (e non 2n perché l’ultimo addendo è nullo per quanto detto precedentemente).

Nel caso di NHPP con legge esponenziale (Weibull), la funzione di densità che si trova è la stessa del caso di NHPP con n fissato, sostituendo t con t_n ed n con N (la v.c. che esprime il numero di arrivi registrati a t_n). Così facendo le formule degli stimatori sono le medesime, con la differenza che λ n / β.

Se si collezionano dati da sistemi multipli, è consigliabile effettuare test per saggiare l’ipotesi che i processi nei diversi sistemi si comportino allo stesso modo. Sia allora t il tempo fissato relativo al j-esimo sistema o, equivalentemente, l’istante dell’ultimo arrivo. Sia N = n₁ + … + n_j il totale delle osservazioni. Poiché i sistemi sono indipendenti, la verosimiglianza del sistema complessivo si può ottenere come prodotto di quella dei singoli. Pertanto, se i processi seguono una P. e hanno una i.f. esponenziale, si ottiene la seguente logverosimiglianza:

Il fatto che le due funzioni di verosimiglianza siano così simili non è scontato ma scaturisce dalla scelta della variabile casuale Weibull. Poiché la funzione di verosimiglianza è simile, succedono due particolarità: uguale per tutti.

Risolvendo otteniamo le seguenti stime, che sono non lineari tranne nel caso in cui t è fissato e:

Analizziamo: questa somma non si ferma più a n-1 perché prima si metteva t_n a denominatore, cioè il tempo dell’ultimo arrivo, ora invece si mette t, cioè il tempo di arresto. Non è detto che il tempo dell’ultimo arrivo coincida con il tempo di arresto ed è questo il motivo per cui l’ultimo addendo della somma non è univocamente nullo.

Una semplificazione si ha facendo riferimento alla verosimiglianza condizionata, ovvero quando si dispone del valore della verosimiglianza. In questo caso, per processi con t fissato o con n fissato si ha rispettivamente:

• M=∑_i n_i
• M=∑_i (n_i – 1)

La stima di θ non cambia (basta mettere M al posto di N), mentre per β si ha:

Per testare che due sistemi abbiano la stessa β si può utilizzare un test, dividendo le due statistiche per i loro rispettivi g.d.l., che sotto l’ipotesi nulla H₀: β₁ = β₂ si distribuisce come una F-Fisher: F= λ e quindi è possibile effettuare il test di Fisher.

Se invece si vuole saggiare se più sistemi, ai quali si può attribuire un andamento di tipo Weibull, hanno lo stesso β, si utilizza il test di Bartlett basato sul rapporto di verosimiglianza R. Essendo β* la media aritmetica dei β, possiamo calcolare tale rapporto come segue:

Con 2LR/a con regione di rifiuto a destra:

Se può ammettersi che tutti i sistemi hanno lo stesso α, la funzione di verosimiglianza si semplifica, ottenendo:

Si può fare inferenza anche attraverso metodi grafici. Ad esempio, si possono riportare in ordinata N(t_i) (cumulata del numero di arrivi) e in ascissa il tempo globale t_i, a partire da 0. In un siffatto grafico, una curvatura evidente verso l’alto sarà sintomo di un deterioramento, viceversa vi è un miglioramento del sistema. Se la curvatura assume un simile andamento, allora il processo può modellarsi con una v.c. Weibull; se invece non ci sono evidenti curvature si può ammettere che il processo sia di rinnovo con tasso di arrivo costante.

Il grafico di Duane riporta invece il cumulative failure rate N(t_k)/t_k rispetto a t_k. In questo caso un andamento crescente significa un incremento della i.f. nel tempo, viceversa un decremento. Assumendo una Weibull i.f.:

La stima dei due parametri può ottenersi tramite interpolazione (e con l’OLS) sul grafico di Duane. Se l’ipotesi di linearità è accettabile, se la pendenza è 0 si può assumere un modello HPP (α=1).

Test di Laplace

Un altro test per HPP è il test di Laplace che è mirato a vedere se l’intensità aumenta o diminuisce. Supponendo di campionare sino ad un numero n di arrivi, dato il tempo dell’ultimo arrivo T = t_n, i primi n – 1 tempi di n arrivi sono statistiche ordinate secondo una distribuzione uniforme nell’intervallo (0, t_n), e dunque hanno media t_n/2 e varianza t_n²/12. Per n>20, possiamo approssimarla ad una distribuzione normale con valore atteso n/2 e deviazione std. t_n/√12 e calcolare il valore L:

Se risulta L > z significa che...