Post on 07-Oct-2020
2Simulazione
Con il termine simulazione si intende la riproduzione del comportamento di un
sistema. In generale, si parla di simulazione sia nel caso in cui viene utilizzato
un modello concreto, sia nel caso in cui viene utilizzato un modello astratto che
riproduce la realta mediante l’uso del computer. Un esempio di modello concreto
e il modello in scala di una nave che viene poi posto in un’apposita vasca per
effettuare prove simulate allo scopo di stimare opportune misure di prestazione. E
chiaro che esistono leggi teoriche della fisica dalle quali ottenere informazioni sulle
prestazioni della nave, ma le analisi di queste leggi e spesso troppo complicata,
per essere effettuata; naturalmente, e anche impraticabile (o quanto meno non
conveniente) la costruzione reale della nave e la prova diretta in mare.
All’interno della Ricerca Operativa, la simulazione utilizza modelli astratti che
vengono costruiti al fine di “replicare” il funzionamento di un sistema. Essa gioca
un ruolo molto importante soprattutto nel progettare un sistema stocastico e nel
definirne le procedure operative: il funzionamento di un sistema e “simulato” uti-
lizzando distribuzioni di probabilita per generare casualmente eventi del sistema
e dal sistema simulato si ottengono osservazioni statistiche sulle prestazioni dello
stesso. Naturalmente affinche cio possa essere realizzato e necessario costruire un
modello di simulazione, che permetta di descrivere le operazioni di un sistema e
come esse devono essere simulate.
Gli aspetti rilevanti che fanno della simulazione uno strumento largamente uti-
lizzato sono legati al fatto che essa permette di
• rappresentare sistemi reali anche complessi tenendo conto anche delle sor-
genti di incertezza;
• riprodurre il comportamento di un sistema in riferimento a situazioni che
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
124 SIMULAZIONE
non sono sperimentabili direttamente.
D’altra parte deve essere sempre tenuto sempre ben presente il fatto che
• la simulazione fornisce indicazioni sul comportamento del sistema, ma non
“risposte” esatte;
• l’analisi dell’output di una simulazione potrebbe essere complessa e potreb-
be essere difficile individuare quale puo essere la configurazione migliore;
• l’implementazione di un modello di simulazione potrebbe essere laboriosa
ed inoltre potrebbero essere necessari elevati tempi di calcolo per effettuare
una simulazione significativa.
2.1 GENERALITA SUI MODELLI DI SIMULAZIONE
Come abbiamo gia osservato, per simulare il comportamento di un sistema e
necessario costruire un modello di simulazione. Il modello dovra essere suffi-
cientementre complesso da rispondere alle esigenze dal caso, ma deve comunque
rimanere il piu semplice possibile. Devono inoltre essere chiari i limiti di utilizzo
del modello stesso.
2.1.1 Elementi di un modello di simulazione
Vediamo ora gli elementi che costituiscono un modello di simulazione.
• Variabili di stato
Innanzitutto ricordiamo che un sistema e descritto in ogni istante di tempo
da un insieme di variabili che prendono nome di variabili di stato. Quindi,
ad esempio, in riferimento ad un sistema a coda, e una variabile di stato
il numero degli utenti presenti nel sistema in un certo istante di tempo.
Ricordiamo, inoltre, che esistono sistemi discreti in cui le variabili cambiano
istantaneamente in corrispondenza di precisi istanti di tempo che sono finiti
oppure appartenenti ad un insieme numerabile e sistemi continui in cui le
variabili variano con continuita rispetto al tempo. Si osservi fin d’ora che la
scelta di un modello continuo o discreto da utilizzare non e necessariamente
obbligata dalla tipologia del sistema; si puo infatti decidere, ad esempio, di
costruire un modello discreto per un sistema continuo, a seconda dello studio
che si vuole effettuare. Un esempio tipico e il caso in cui nel rappresentare
una linea ferroviaria, la posizione del treno puo essere descritta da una
variabile reale che fornisce la distanza dalla stazione di origine, oppure da
variabili binarie che descrivono lo stato libero–occupato di ciascuna delle
sezioni di blocco in cui e divisa la linea.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
GENERALITA SUI MODELLI DI SIMULAZIONE 125
• Eventi
Si definisce evento un qualsiasi accadimento istantaneo che fa cambiare il
valore di almeno una delle variabili di stato. L’arrivo di un utente ad un
sistema a coda e un evento, cosı come il completamento di un servizio.
Esistono eventi esterni al sistema (eventi esogeni) ed eventi interni (eventi
endogeni). Ad esempio, l’inizio del servizio ad un utente che e in coda in
un sistema a coda e un evento endogeno, perche interno al sistema; l’arrivo
di un utente ad un sistema a coda e un evento esogeno.
• Entita e attributi
Le entita sono singoli elementi del sistema che devono essere definiti. Un
esempio di entita e un utente presso un sistema a coda, oppure puo essere
un servente. Nel primo caso l’entita fluisce all’interno del sistema e si parla
di entita dinamica, nel secondo caso si parla di entita statica.
Le entita possono essere caratterizzate da attributi che forniscono un valore
di un dato assegnato all’entita stessa. Ad esempio, in un sistema a coda
monoservente dove le entita sono il servente e gli utenti, un attributo di
un’entita “utente” potrebbe essere il suo tempo di arrivo al sistema, mentre
il servente e caratterizzato dall’attributo “status” che puo assumere valore
di “libero” o “occupato”. E chiaro che alcuni attributi possono essere di
interesse in alcuni casi e non in altri.
Le entita possono essere raggruppate in classi che sono insiemi di entita dello
stesso tipo, ovvero si possono raggruppare le entita in base ad attributi.
Se, ad esempio, consideriamo persone di sesso maschile e femminile come
utenti di un sistema a coda, essendo le entita le persone, esse possono essere
raggruppate in dua classi in base all’attributo “sesso”.
• Risorse
Le risorse sono elementi del sistema che forniscono un servizio alle entita.
Un’entita puo richiedere una o piu unita di risorsa e se questa non e di-
sponibile l’entita dovra mettersi, ad esempio, in una coda in attesa che si
renda disponibile, oppure intraprendere un’altra azione. Se invece la risorsa
e disponibile, essa viene “catturata” dall’entita, “trattenuta” per il tempo
necessario e poi “rilasciata”. Un esempio di risorsa potrebbe essere data
da un operaio che sovrintende il funzionamento di una macchina che non
puo funzionare senza l’operaio stesso; quando e richiesto l’utilizzo di questa
macchina, se la risora “operaio” e disponibile allora l’esecuzione del lavoro
e effettuata altrimenti si attende che ci sia risorsa (operaio) disponibile.
L’operaio verra “trattenuto” per la durata dell’esecuzione del lavoro e poi
“rilasciato”. Si osservi che, in generale, un elemento del modello potrebbe
essere considerato parimenti un’entita o una risorsa. Questo, ovviamente,
dipende da come si e scelto di costruire un modello.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
126 SIMULAZIONE
• Attivita e ritardi
Un’attivita e un’operazione la cui durata e nota a priori all’inizio dell’ese-
cuzione dell’attivita stessa. Tale durata puo essere una costante, un valore
aleatorio generato da una distribuzione di probabilita, oppure data in input
o calcolata in base ad altri eventi che accadono nel sistema. Un esempio e
dato dal tempo di servizio in un sistema a coda.
Un ritardo e un periodo di tempo di durata indefinita che e determinata
dalle condizioni stesse del sistema. Il tempo che un’entita trascorre presso
una coda prima che si liberi una risorsa della quale necessita e un ritardo.
2.1.2 Classificazione dei modelli si simulazione
I modelli di simulazione si possono classificare in base a diversi criteri; una prima
distinzione gia vista e tra
• modelli continui, in cui le variabili variano con continuita;
• modelli discreti, in cui il valore delle variabili cambia in ben definiti istanti
di tempo.
Un’altra distinzione e tra:
• modelli statici, che rappresentano un sistema in un particolare istante di
tempo;
• modelli dinamici, che rappresentano un sistema in evoluzione nel tempo.
Infine, si possono distinguere
• modelli deterministici, che non contengono componenti probabilistici;
• modelli stocastici, che presentano elementi soggetti ad aleatorieta.
In questa trattazione considereremo modelli di simulazione discreti, dinamici,
stocastici che vengono comunemente chiamati modelli di simulazione ad eventi
discreti. Molte applicazioni sono ben rappresentate da modelli di questo tipo
ed inoltre approssimando variazioni continue con variazioni discrete e possibile
utilizzare modelli ad eventi discreti anche per approssimare il comportamento di
sistemi continui semplificando quindi molto l’analisi.
2.1.3 Simulazione ad eventi discreti
Nella simulazione ad eventi discreti il sistema e rappresentato, nella sua evoluzione
nel tempo, con variabili che cambiano instantaneamente il loro valore in ben
definiti istanti di tempo appartenenti ad un insieme numerabile. Questi istanti
sono quelli nei quali accadono gli eventi. E chiaro che, essendo questi modelli
di natura dinamica, e necessario registrare, ovvero tenere memoria, del tempo
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
GENERALITA SUI MODELLI DI SIMULAZIONE 127
(simulato) che procede. In particolare sara necessario definire un meccanismo di
avanzamento del tempo per far procedere il tempo simulato da un valore ad un
altro. La variabile che in un modello di simulazione fornisce il valore corrente del “Simulation
clock”tempo simulato si chiama “simulation clock”, ed esistono due modi per definire
il suo avanzamento:
• avanzamento del tempo al prossimo evento,
• avanzamento del tempo ad incrementi prefissati.
Il primo e quello piu diffuso ed e quello a cui faremo riferimento. In questo caso
il “simulation clock” e inizializzato a zero e viene avanzato al tempo dell’accadi-
mento del primo degli eventi futuri; poi il sistema viene aggiornato tenendo conto
dell’evento che e accaduto, si aggiornano i tempi degli eventi futuri e si itera il
procedimento. A differenza dell’avanzamento ad incrementi prefissati, i periodi
di inattivita non vengono considerati.
Un esempio puo essere visto considerando un sistema di code in cui gli eventi sono
l’arrivo di un cliente, la conclusione di un servizio; entrambi sono eventi perche
provocano il cambiamento di valore di qualche variabile di stato. Il meccanismo
di avanzamento del tempo segue in questo caso l’accadere di questi due eventi
nell’ordine cronologico in cui essi si verificano.
Un esempio di simulazione ad eventi discreti
Vediamo, ora, un semplice esempio di come si realizza un simulazione ad eventi
discreti. Consideriamo a tale scopo un sistema a coda costituito da una coda
e da un singolo servente e supponiamo che i tempi di interarrivo siano unifor-
memente distribuiti tra 1 e 3 minuti e che anche i tempi di servizio siano uni-
formemente distribuiti tra 0.5 e 2 minuti. Vediamo, ora, come si puo effettuare
una simulazione di questo sistema. Poiche si tratta di un sistema regolato da
due processi stocastici (gli arrivi e i servizi) per generare gli eventi e necessario
generare osservazioni casuali dalle due distribuzioni di probabilita che regolano
i due processi (come questo puo essere effettuato sara oggetto di considerazioni
successive nel paragrafo 2.3). Supponiamo di avere a disposizione le due liste
che forniscono, rispettivamente i tempi di interarrivo generati casualmente dalla
distribuzione corrispondente e i tempi di servizio anch’essi generati casualmente
dalla distribuzione corrispondente:
Tempi di interarrivo Tempi di servizio
1.9 1.7
1.3 1.8
1.1 1.5
1.0 0.9...
...
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
128 SIMULAZIONE
Supponendo che al tempo t = 0 nessun utente e presente nel sistema. Osservando
i valori campionati riportati nelle due liste, si ricava facilmente la successione degli
eventi:
Tempo t Eventi
1.9 arriva un utente inizia il servizio
3.2 arriva un utente e si pone in coda
3.6 finisce un servizio e il primo utente in coda inizia il servizio
4.3 arriva un utente e si pone in coda
5.3 arriva un utente e si pone in coda
5.4 finisce un servizio e il primo utente in coda inizia il servizio...
...
Limitando questa semplice simulazione al tempo t = 5.4 (in modo che due utenti
sono entrati e hanno completato il servizio), possiamo calcolare, ad esempio, il
tempo medio di permanenza nel sistema: il primo utente rimane nel sistema
1.7 minuti, il secondo 2.2 minuti e quindi il valore medio e 1.95. Questa stima,
ovviamente non ha alcun senso perche ottenuta dalla particolare sequenza di
numeri casuali delle due liste. Quindi, se l’esempio da un lato vuole mettere
evidenza il meccanismo di una simulazione ad eventi discreti, dall’altro mette fin
d’ora in evidenza un errore che si potrebbe commettere nel reputare affidabili
i risultati di una sola esecuzione e che ha avuto una durata arbitraria. D’altra
parte c’e anche da tener presente che se siamo interessati a valutare misure di
prestazioni del sistema a regime, ovvero quando sono state raggiunte condizioni di
stazionarieta, sara necessario non prendere in considerazione il sistema durante il
periodo iniziale di transitorio. Queste problematiche rappresentano un elemento
chiave di ogni simulazione e saranno considerate in dettaglio nel seguito.
2.1.4 Schema dello studio di un problema basato sulla simulazione
In questo paragrafo riportiamo uno schema che descrive la successione delle varie
fasi che caratterizzano uno studio basato sulla simulazione.
1. Analisi del problema
Consiste nel comprendere il problema cercando di capire quali sono gli scopi
dello studio e di identificare quali sono le componenti essenziali e quali sono
le misure di prestazione che interessano. Naturalmente, se una versione
del sistema e gia operativa, si deve osservare tale sistema per dedurne le
caratteristiche fondamentali.
2. Formulazione del modello di simulazione
Poiche stiamo trattando sistemi stocastici, per formulare un modello di si-
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
GENERALITA SUI MODELLI DI SIMULAZIONE 129
mulazione e necessario conoscere le distribuzioni di probabilita delle quan-
tita di interesse. Infatti, per generare vari scenari rappresentativi di come
un sistema funziona, e essenziale che una simulazione generi osservazioni
casuali da queste distribuzioni. Ad esempio, nei sistemi a coda e necessaria
la distribuzione dei tempi di interarrivo e i tempi di servizio; nella gestione
delle scorte e necessaria la distribuzione della richiesta dei prodotti e la
distribuzione del tempo tra un’ordine e il ricevimento della merce; nella ge-
stione dei sistemi di produzione con macchine che occasionalmente possono
guastarsi, sara necessario conoscere la distribuzione del tempo fino a che
una macchina si guasta e la distribuzione dei tempi di riparazione. General-
mente e possibile solo stimare queste distribuzioni derivandole, ad esempio,
dall’osservazione di sistemi simili gia esistenti. Se dall’analisi dei dati si
vede che la forma di questa distribuzione approssima una distribuzione ti-
po standard, si puo utilizzare la distribuzione teorica standard effettuando
un test statistico per verificare se i dati possono essere rappresentati bene
mediante quella distribuzione di probabilita. Se non esistono sistemi simili
dai quali ottenere dati osservabili si deve far ricorso ad altre fonti di infor-
mazioni: specifiche delle macchine, manuali di istruzioni delle stesse, studi
sperimentali, etc.
La costruzione di un modello di simulazione e un procedimento complesso.
In particolare, facendo riferimento alla simulazione ad eventi discreti, la
costruzione di un modello prevede le seguenti fasi:
(a) Definizione delle variabili di stato.
(b) Identificazione dei valori che possono essere assunti dalle variabili di
stato.
(c) Identificazione dei possibili eventi che fanno cambiare lo stato del
sistema.
(d) Realizzazione di una misura del tempo simulato, “simulation clock”,
che registra lo scorrimento del tempo simulato.
(e) Realizzazione di un metodo per generare casualmente gli eventi.
(f) Identificazione delle transizioni di stato generate dagli eventi.
3. Analisi del modello di simulazione
Nella fase di analisi del modello deve essere verificata l’accuratezza del mo-
dello realizzato con diverse modalita. Di solito cio viene fatto attraverso
un’analisi concettuale del modello che puo essere effettuata insieme agli
esperti del settore applicativo in modo da evidenziare eventuali errori e/o
omissioni.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
130 SIMULAZIONE
4. Scelta del software e costruzione di un programma
Dopo aver costruito il modello, esso deve essere tradotto in un programma.
A tale scopo e possibile utilizzare diversi strumenti.
• Linguaggi “general purpose”.
Linguaggi come JAVA, C++, FORTRAN, etc. Erano molto utilizzati
alla nascita della simulazione ma richiedono molto tempo di program-
mazione e quindi si preferisce, in genere, utilizzare linguaggi specifici
per la simulazione.
• Linguaggi di simulazione generali.
Forniscono molte caratteristiche necessarie per realizzare un modello
di simulazione riducendo cosı il tempo di realizzazione; esempi sono
SIMAN, MODSIM, GPSS, SIMSCRIPT, etc. Anche se meno flessibili dei
linguaggi “general purpose” sono il modo piu naturale per realizzare
un modello di simulazione.
• Simulatori.
Sono packages per la simulazione orientati alle applicazioni. Esistono
numerosi pacchetti software di tipo interattivo per la simulazione co-
me ARENA, SIMIO, WITNESS, EXTEND, MICRO SAINT. Alcuni sono
abbastanza generali anche se dedicati a specifici tipi di sistemi come
impianti industriali, sistemi di comunicazione, altri invece sono molto
specifici come, ad esempio, nel caso di simulatori di centrali nucleari
o di simulatori della fisiologia cardiovascolare. I simulatori permetto-
no di costruire un programma di simulazione utilizzando menu grafici
senza bisogno di programmare. Sono abbastanza facili da imparare e
un inconveniente che molti di essi hanno e di essere limitati a model-
lare quei sistemi previsti dalle loro caratteristiche standard. In ogni
caso alcuni simulatori prevedono la possibilita di incorporare routi-
nes scritte in un linguaggio general purpose per trattare elementi non
standard. Spesso hanno anche capacita di animazione per mostrare
la simulazione in azione e questo permette di illustrare facilmente la
simulazione anche a persone non esperte.
• Fogli elettronici (spreadsheets).
Quando si hanno problemi di piccole dimensioni si possono anche uti-
lizzare fogli elettronici, come ad esempio Excel, per avere un’idea del
funzionamento di un sistema.
5. Validazione del modello di simulazione
Nella fase successiva e necessario verificare se il modello che e stato realizza-
to fornisce risultati validi per il sistema in esame. Piu in particolare si deve
verificare se le misure di prestazione del sistema reale sono bene approssi-
mate dalle misure generate dal modello di simulazione. Cio e molto difficile
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
GENERALITA SUI MODELLI DI SIMULAZIONE 131
da effettuare, specialmente in fase di progettazione quando il sistema reale
non esiste.
6. Progettazione della simulazione
Prima di passare all’esecuzione della simulazione e necessario decidere co-
me condurre la simulazione. Spesso una simulazione e un processo che
evolve durante la sua realizzazione e dove i risultati iniziali aiutano a con-
durre la simulazione verso configurazioni piu complesse. Ci sono inoltre
problematiche di tipo statistico:
• la determinazione della lunghezza del transitorio del sistema prima di
raggiungere condizioni di stazionarieta, momento dal quale si inizia
a raccogliere dati se si vogliono misure di prestazione del sistema a
regime;
• la determinazione della lunghezza della simulazione (durata) dopo che
il sistema ha raggiunto l’equilibrio. Infatti, si deve sempre tener pre-
sente che la simulazione non produce valori esatti delle misure di pre-
stazione di un sistema in quanto ogni singola simulazione puo essere
vista come un “esperimento statistico” che genera osservazioni stati-
stiche sulle prestazioni del sistema. Queste osservazioni sono poi uti-
lizzate per produrre stime delle misure di prestazione e naturalmente
aumentando la durata della simulazione puo aumentare la precisione
di queste stime;
• la determinazione del numero delle repliche di una simulazione. In-
fatti, sara necessario effettuare un certo numero di run indipendenti
affinche si raggiunga l’accuratezza desiderata nella stima delle misure
di interesse.
7. Esecuzione della simulazione e analisi dei risultati
L’output della simulazione fornisce stime statistiche delle misure di presta-
zione di un sistema. Un punto fondamentale e che ogni misura sia accompa-
gnata dall’“intervallo di confidenza” all’interno del quale essa puo variare.
Questi risultati potrebbero evidenziare subito una configurazione del si-
stema migliore delle altre, ma piu spesso verranno identificate piu di una
configurazione candidata ad essere la migliore. In questo caso potrebbero
essere necessarie ulteriori indagini per confrontare queste configurazioni.
8. Presentazione delle conclusioni
In conclusione, e necessario redigere una relazione ed una presentazione
che riassuma lo studio effettuato, come e stato condotto e includendo la
documentazione necessaria. Includere nella presentazione un’animazione di
una simulazione e di solito molto efficace.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
132 SIMULAZIONE
2.1.5 Applicazioni tipiche della simulazione
La simulazione e uno strumento molto flessibile: puo essere utilizzata per studiare
la maggior parte dei sistemi esistenti. E impossibile enumerare tutte le aree
specifiche in cui la simulazione puo essere utilizzata. Come esempi, riportiamo,
di seguito, solo alcune importanti tipiche categorie di applicazioni in cui si usa la
simulazione.
• Progettazione e definizione delle procedure operative di un sistema di ser-
vizio.
• Gestione di sistemi di scorte.
• Progetto e definizione delle procedure operative di sistemi di produzione.
• Progetto e funzionamento del sistemi di distribuzione.
• Analisi dei rischi finanziari.
• Gestione dei progetti.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
ANALISI DELL’INPUT 133
2.2 ANALISI DELL’INPUT
Per effettuare una simulazione di un sistema che presenta elementi stocastici
e necessario specificare le distribuzioni di probabilita che regolano i processi che
caratterizzano il sistema stesso. Ad esempio, nei sistemi a coda devono essere note
le distribuzioni di probabilita dei tempi di interarrivo e dei tempi di servizio. Se e
possibile raccogliere dati reali (osservazioni) sulle variabili aleatorie di interesse,
essi possono essere utilizzati per determinare queste distribuzioni facendo uso di
tecniche di inferenza statistica. Questa prima fase viene di solito chiamata analisi
dell’input e sara oggetto di studio di questo paragrafo.
Una volta stabilite tali distribuzioni, la simulazione procede generando valori
casuali da queste distribuzioni, ovvero, durante ogni esecuzione, la simulazione
genera osservazioni casuali di variabili aleatorie distribuite secondo particolari
distribuzioni di probabilita. Successivamente, vengono utilizzate tecniche stati-
stiche per interpretare i risultati ottenuti da una simulazione (analisi dell’output).
Questi ultimi due aspetti saranno trattati rispettivamente nei successivi paragrafi
2.3 e 2.4.
2.2.1 Generalita
In generale, nello studio di un fenomeno riguardante un insieme di elementi (po-
polazione) che presenta caratteristiche aleatorie, molto spesso si dispone solo di
informazioni su una parte di essi (campione) e si vogliono dedurre proprieta ge-
nerali riguardanti l’intera popolazione. L’inferenza statistica si occupa di questa
problematica e riveste un importante strumento di analisi. Esula dallo scopo di
queste note una trattazione degli elementi dell’inferenza statistica che verranno
assunti come noti. Una breve rassegna di alcuni importanti risultati e riportata
nell’Appendice C. Richiamiamo di seguito solo alcune importanti osservazioni e
le ben note definizioni di media campionaria e varianza campionaria
Solitamente viene fatta l’assunzione che esiste una distribuzione di probabilita
della popolazione nel senso che se da essa vengono estratti casualmente alcu-
ni elementi, ad essi sono associate variabili aleatorie indipendenti identicamente
distribuite secondo tale distribuzione. In questo senso, un insieme di variabili
aleatorie X1, . . . , Xn di variabili aleatorie indipendenti tutte con la stessa distri-
buzione si dice campione di questa distribuzione. L’interesse principale risiede
nella possibilita di dedurre caratteristiche della distribuzione non nota sulla base
dei dati a disposizione. Naturalmente ci sono casi in cui della distribuzione della
popolazione non si conosce nulla (se non il fatto che essa e discreta o continua),
mentre in altri casi la distribuzione e nota ma non sono noti alcuni suoi parametri.
Formalmente, sia dato un campione X1, . . . , Xn estratto da una popolazione,
ovvero le Xi sono variabili aleatorie indipendenti identicamente distribuite, e sia
µ e σ2 rispettivamente la loro media e la loro varianza (ovvero la media e la
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
134 SIMULAZIONE
varianza della popolazione). La media campionaria e data daMedia cam-
pionaria
Xn =1
n
n∑i=1
Xi.
Xn e una variabile aleatoria funzione delle Xi e si verifica facilmente che risulta
E(Xn) = µ e V ar(Xn) =σ2
n.
La varianza campionaria e data daVarianza
campiona-
ria s2n =
1
n− 1
n∑i=1
(Xi − Xn
)2e si verifica facilmente che risulta E(s2
n) = σ2.
2.2.2 Scelta delle distribuzioni di input
Se e possibile raccogliere dati reali sulle variabili aleatorie di interesse, essi possono
essere utilizzati per determinare una distribuzione di probabilita che meglio li
rappresenta secondo tre metodi:
1. i dati sono usati direttamente nella simulazione (“trace driver simulation”).
2. i dati sono raccolti per generare una distribuzione empirica, ovvero per
definire una funzione di distribuzione empirica che verra usata per produrre
l’input della simulazione;
3. i dati raccolti sono utilizzati per definire una distribuzione teorica. Vengono
utilizzate tecniche statistiche per analizzare se una distribuzione teorica tra
quelle note sia adatta a rappresentare i dati, effettuando i test di ipotesi
per verificare la rappresentativita della distribuzione ipotizzata (problema
del “fitting”).
Il primo approccio ha senso solamente quando si possono raccogliere grandi quan-
tita di dati rappresentativi del funzionamento del sistema; ha l’ovvio difetto di
rappresentare il “passato” ed e usato raramente; puo essere utile per effettua-
re una validazione del modello, ovvero per confrontare il modello con il sistema
reale, ma non permette un’analisi previsionale.
Il secondo approccio elimina questo inconveniente poiche, almeno per distribuzio-
ni continue, puo essere ottenuto ogni valore compreso tra il minimo e il massimo
osservati.
Se si puo determinare una distribuzione teorica che si adatta bene ai dati, il terzo
approccio e quello preferibile. I motivi per cui una distribuzione teorica in genere
e preferibile a una empirica sono i seguenti:
• le distribuzioni empiriche possono avere irregolarita (specialmente se i dati
sono scarsi) mentre le distribuzioni teoriche sono piu “smooth”, nel senso
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
ANALISI DELL’INPUT 135
che tendono a regolarizzare i dati e rappresentano un comportamento ge-
nerale;
• le distribuzioni empiriche non permettono di generare valori al di fuori del
range di valori osservati, mentre le misure di prestazione possono, a volte,
dipendere anche da eventi “eccezionali” che corrispondono a valori fuori da
tale range;
• le distribuzioni teoriche sono un modo compatto di rappresentare un in-
sieme di valori, mentre in una distribuzione empirica, se ci sono n dati
disponibili, si ha bisogno di 2n valori per rappresentarla: il dato e le corri-
spondenti probabilita cumulative (si hanno quindi grandi quantita di dati
da memorizzare);
• le distribuzioni teoriche si possono variare piu facilmente. Ad esempio se la
distribuzione esponenziale degli arrivi di un sistema di code ha media pari a
1/λ = 5, per effettuare una diminuzione del 20% sara sufficiente considerare
1/λ = 4.9.
Tuttavia esistono situazioni in cui nessuna distribuzione teorica si adatta ai dati
osservati e allora in questo caso si deve usare una distribuzione empirica.
Un difetto dell’uso di distribuzioni teoriche sta nel fatto che esse possono generare
anche valori molto grandi (anche se con probabilita molto piccole), quando nella
pratica questi non vengono mai assunti realmente.
2.2.3 Distribuzioni empiriche
Supponiamo di disporre di n osservazioni X1, . . . , Xn di una variabile aleatoria e
di voler costruire, a partire da esse, una distribuzione continua. Supponiamo di
aver ordinato le Xi per valori crescenti e sia X(i) l’i-esima osservazione in ordine
crescente, ovvero risulti
X(1) ≤ X(2) ≤ . . . ≤ X(n).
Si puo costruire la distribuzione empirica come una distribuzione continua lineare
a tratti, cosı definita:
F (x) =
0 se x < X(1)
i− 1
n− 1+
x−X(i)
(n− 1)(X(i+1) −X(i))se X(i) ≤ x < X(i+1),
i = 1, . . . , n− 1
1 se X(n) ≤ x
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
136 SIMULAZIONE
Si osservi che per ogni i vale F (X(i)) =i− 1
n− 1che e approssimativamente (per n
grande) la proporzione delle Xi che sono minori di X(i).
Esempio 2.2.1 Disponendo dei seguenti valori osservati: 1, 0.4, 4, 2, 2.5, 3.6, 3 costruire
il grafico della distribuzione empirica. Dopo aver ordinato le osservazioni si ottiene il
grafico della F (x) riportato nella Figura 2.2.1.
1/6
1/3
1/2
2/3
5/6
1
1
2
3 4
0.4
2.5
3.6
Figura 2.2.1 Grafico della distribuzione empirica dell’Esempio 2.2.1
Come abbiamo gia osservato, uno svantaggio nell’utilizzare una distribuzione em-
pirica e che le variabili aleatorie generate da essa durante un’esecuzione di una
simulazione non possono essere mai piu piccole di X(1) o piu grandi di X(n).
Analogamente si possono costruire distribuzioni empiriche per distribuzioni di-
screte; infatti, e sufficiente per ogni x definire p(x) come “proporzione” delle Xi
che sono uguali ad x.
2.2.4 Distribuzioni teoriche
• Distribuzioni Continue
Le distribuzioni teoriche continue alle quali si puo fare riferimento nella co-
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
ANALISI DELL’INPUT 137
struzione di un modello di simulazione sono molte. Quelle piu comunemen-
te utilizzate sono la distribuzione uniforme, la distribuzione esponenziale,
la distribuzione gamma, la distribuzione normale, la distribuzione lognor-
male, la distribuzione di Weibull, la distribuzione beta, la distribuzione
triangolare.
In realta spesso si tratta di famiglie di distribuzioni in quanto sono presenti
uno o piu parametri che possono essere classificati in:
· parametro di posizione: specifica un punto del range della distribuzione
e una sua variazione provoca solamente una traslazione;
· parametro di scala: specifica l’unita di misura dei valori,
· parametro di forma: specifica l’andamento della distribuzione.
• Distribuzioni Discrete
Le distribuzioni teoriche discrete che vengono di solito utilizzate come in-
put di una simulazione sono: la distribuzione uniforme, la distribuzione
di Bernoulli, la distribuzione binomiale, la distribuzione geometrica, la
distribuzione di Poisson, la distribuzione binomiale negativa.
Per una descrizione dettagliata di ogni singola distribuzione di probabilita e del-
le caratteristiche specifiche, si rimanda ad un qualsiasi testo di Calcolo delle
Probabilita.
2.2.5 Scelta di una distribuzione teorica
Determinare quale distribuzione teorica e adatta a rappresentare dei dati e un
problema complesso. Un modo efficiente per effettuare questa scelta puo essere
cosı schematizzato:
Dopo una preliminare verifica dell’indipendenza delle osservazioni, si cerca di in-
dividuare una o piu famiglie di distribuzioni candidate, stimando poi nella fase
successiva i parametri di queste distribuzioni. A questo punto e necessario effet-
tuare una verifica della rappresentativita dei dati reali da parte della distribuzioni
scelta. Se tale verifica non ha successo, e necessario individuare una differente
famiglia di distribuzioni ed eventualmente iterare il processo fino a che la verifica
e soddisfatta.
Indipendenza delle osservazioni
Preliminarmente e necessario verificare l’indipendenza delle osservazioni in quan-
to questa e un’assunzione essenziale per l’utilizzo di tecniche statistiche quali la
stima della massima verosimiglianza o il test chi–quadro. Un primo strumen-
to di analisi e basato su una tecnica grafica. Siano X1, . . . , Xn le osservazioni
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
138 SIMULAZIONE
Verif ica dell段ndipendenza
delle osservazioni
Individuazione di una
famiglia di distribuzioni
Stima dei parametri
della distribuzione
Verif ica della rappresentat ivit�
elencate cosı come sono state osservate nel tempo; un modo possibile per ave-
re un’idea informale sull’indipendenza consiste nel valutare la correlazione fra
diverse osservazioni. Sia
ρj =
n−j∑i=1
(Xi − Xn)(Xi+j − Xn)
(n− j)s2n
la stima del coefficiente di correlazione ρj di Xi e Xi+j , ovvero di due osservazioni
distanti j. Se le osservazioni sono indipendenti allora il coefficiente di correlazione
e nullo, cioe ρj = 0 per ogni j = 1, . . . , n− 1. Tuttavia poiche ρj e una stima di
ρj , anche nel caso di osservazioni indipendenti ρj potrebbe essere non nullo. Ci
si aspetta, comunque che esso sia prossimo a zero, e quindi possiamo dire che se
ρj e diverso da zero in maniera significativa, allora le Xi non sono indipendenti.
Ci sono due modi grafici per verificare informalmente se le Xi sono indipendenti:
il grafico di ρj al variare di j e il diagramma di dispersione delle osservazioni
X1, . . . , Xn, ovvero le coppie (Xi, Xi+1) con i = 1, 2, . . . , n − 1. In caso di os-
servazioni indipendenti i punti dovrebbero risultare distribuiti casualmente sul
piano, altrimenti, in presenza di correlazioni, essi saranno concentrati intorno a
rette.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
ANALISI DELL’INPUT 139
Individuazione di una famiglia di distribuzioni
Una volta verificata l’indipendenza delle osservazioni, il passo successivo e quello
di individuare una distribuzione da scegliere come input della simulazione che
sia rappresentativa della variabile aleatoria in ingresso alla simulazione. In una
prima fase vorremmo individuare una famiglia generale senza occuparci, per ora,
dei suoi parametri. In alcuni casi, quando esiste una conoscenza “a priori” del
fenomeno che la variabile aleatoria rappresenta, essa puo essere utilizzata per
ottenere la distribuzione. Cio e fatto su base teorica e non richiede osservazioni.
Ad esempio, se supponiamo che dei clienti arrivano ad un sistema di servizio uno
alla volta e che il numero dei clienti che arrivano in intervalli disgiunti e indi-
pendente, allora ci sono motivi teorici per assumere che i tempi di interarrivo
siano variabili aleatorie indipendenti identicamente distribuite secondo la distri-
buzione esponenziale. Oppure, puo anche accadere che la conoscenza “a priori”
permetta solo di escludere alcune distribuzioni. Tuttavia, nella pratica spesso
queste informazioni “a priori” non sono disponibili, o comunque non sufficienti.
Quello che si fa piu frequentemente e ricorrere a due strumenti di analisi stati-
stica: le statistiche riassuntive delle osservazioni e i grafici dell’andamento delle
osservazioni.
• Statistiche riassuntive
Dalle osservazioni e possibile ricavare stime di parametri dalle quali cercare
di individuare una famiglia di distribuzioni che meglio realizza il fitting
dei dati. I parametri che di solito vengono presi in considerazione sono i
seguenti:
– l’intervallo [X(1), X(n)] che ha per estremi il piu piccolo e il piu grande
valore osservati e che approssima il range della distribuzione;
– la stima della media µ data Xn =1
n
n∑i=1
Xi;
– la stima della mediana data da{X(n+1)/2 se n e dispari
[X(n/2) +X((n/2)+1)]/2 se n e pari;
– la stima della varianza σ2 data da
s2n =
n∑i=1
(Xi − Xn)2
n− 1;
– stime di misure di variabilita:
· nel caso continuo, la stima del rapporto cv =√σ2/µ
(coefficiente di variazione) data da cv =√s2n/Xn
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
140 SIMULAZIONE
· nel caso discreto, la stima del rapporto τ = σ2/µ
data da τ = s2n/Xn
– la stima del grado di asimmetria ν =E[(X − µ)3]
(σ2)3/2
data da ν =
n∑i=1
[(Xi − Xn)3]/n
(s2n)3/2
.
Un confronto tra media e mediana puo farci capire se considerare la distri-
buzione simmetrica o no; questo perche nel caso di distribuzioni continue
simmetriche, media e mediana coincidono, come ad esempio nel caso della
distribuzione normale (nel caso di distribuzioni discrete questo e vero solo
se il numero dei valori distinti che possono essere assunti e pari, altrimenti
sono solo approssimativamente uguali). E importante tener presente che si
hanno solo le stime dei parametri, pertanto anche nel caso di distribuzioni
continue simmetriche, stima di media e mediana possono non essere esat-
tamente uguali.
Per quanto riguarda la misura di variabilita data nel caso continuo dal
rapporto cv, si ha che cv = 1 per la distribuzione esponenziale indipenden-
temente dal parametro; nel caso discreto si considera invece il rapporto τ e
si ha che τ = 1 per la distribuzione di Poisson e τ < 1 per la distribuzione
binomiale.
Il grado di asimmetria ν vale zero per distribuzioni simmetriche mentre se
ν > 0 la distribuzione ha asimmetria verso destra (ν = 2 per la distribuzione
esponenziale); se ν < 0 la distribuzione ha asimmetria verso sinistra.
• Uso di grafici (tecnica dell’istogramma)
Per distribuzioni continue e molto utile costruire un istogramma di valori
assunti dalla variabile aleatoria sulla base delle osservazioni che si hanno
a disposizione. Un istogramma puo essere considerato una “stima grafica”
del grafico della densita di probabilita corrispondente alla distribuzione dei
dati X1, . . . , Xn. Le funzioni densita di probabilita tipiche hanno forme
che in molti casi sono ben riconoscibili e quindi un istogramma puo fornire
utili indicazioni. Ricordiamo che per costruire un istogramma si suddivide
l’intervallo formato dal minimo e dal massimo dei valori assunti dai dati, in
k intervalli disgiunti adiacenti
[b0, b1), [b1, b2), . . . , [bk−1, bk)
di uguale ampiezza ∆b = bi − bi−1. Per j = 1, 2, . . . , k si definisce hj il
numero delle osservazioni che cadono nel j-esimo intervallo diviso il numero
totale delle osservazioni, ovvero la proporzione delle Xi contenute nel j-
esimo intervallo [bj−1, bj). Si definisce la funzione
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
ANALISI DELL’INPUT 141
h(x) =
0 se x < b0
hj se bj−1 ≤ x < bj , j = 1, 2, . . . , k
0 se x ≥ bk.
Il grafico di h(x) e costante a tratti e puo fornire una buona indicazione
sul tipo di distribuzione che ha la variabile aleatoria in questione, confron-
tandolo con i grafici delle densita di probabilita ignorando, per il momento,
posizione e scala, ma considerando solo la forma.
Mostriamo ora le motivazioni che sono alla base del fatto che la forma di
h(x) dovrebbe “somigliare” alla densita di probabilita f dei dati. A questo
scopo, sia X una variabile aleatoria con densita di probabilita data da f .
Allora per ogni j fissato (j = 1, 2, . . . , k), applicando il teorema della media
si ha
P (bj−1 ≤ X ≤ bj) =
∫ bj
bj−1
f(x)dx = ∆b f(y)
per qualche y ∈ (bj−1, bj). D’altra parte hj approssima P (bj−1 ≤ X ≤ bj)
che e il valore di h(y) perche h(x) = hj per ogni x ∈ [bj−i, bj); si ha dunque:
h(y) = hj ≈ ∆bf(y). Quindi, h(y) e approssimativamente proporzionale a
f(y), ovvero h e f hanno approssimativamente forma simile.
Una difficolta e data dall’assenza di criteri generali per scegliere k. C’e una
regola detta regola di Sturges che suggerisce di scegliere k = b1 + log2 nc. Regola di
SturgesTale regola non e pero utile in generale e si raccomanda piuttosto di provare
differenti valori di ∆b e scegliere il piu piccolo che fornisce un istogramma
“smooth”. La scelta di ∆b e problematica in genere; infatti se si sceglie
troppo piccolo, l’istogramma sara molto poco uniforme (frastagliato); se si
sceglie troppo grande l’istogramma puo avere una forma a blocchi e il vero
andamento della densita che stiamo cercando sara mascherata in quanto i
dati sono sovraggregati.
Nel caso di variabili discrete si puo ugualmente rappresentare la distribuzio-
ne di probabilita usando un istogramma; per ogni valore xj che puo essere
assunto dai dati, sia hj = nj/n dove nj e il numero delle occorrenze di tali
valori, ovvero hj e la proporzione delle Xi che sono uguali a xj . Si tracciano
poi le barre verticali di altezza hj in corrispondenza di xj .
Stima dei parametri
Dopo aver individuato una o piu famiglie di distribuzioni candidate a rappre-
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
142 SIMULAZIONE
sentare i dati osservati, e necessario determinare i parametri di queste distri-
buzioni in modo che siano completamente definite e utilizzabili in ingresso ad
una simulazione. Gli stessi dati utilizzati per individuare le famiglie di distri-
buzioni sono utilizzati per stimare i loro parametri. La stima dei parametri e
un argomento vasto che si assume noto. Alcuni cenni riguardanti lo stimatore
di massima verosimiglianza (Maximum Likelihood Estimator) sono riportati nel
paragrafo C.3 dell’Appendice C.
Verifica della rappresentativita della distribuzione di probabilita
Dopo aver individuato una o piu distribuzioni di probabilita candidate e i relativi
parametri, si devono esaminare queste distribuzioni di probabilita per verificare
se esse rappresentano bene i dati osservati. A tale scopo, si utilizzano di solito
procedure euristiche basate su confronti grafici e test statistici.
• Procedure grafiche
Per distribuzioni continue, si confronta l’istogramma dei dati con il grafico
della densita di probabilita della distribuzione di probabilita ipotizzata,
oppure, per distribuzioni discrete, si confronta l’istogramma con la funzione
p(x) della distribuzione ipotizzata.
Un altro possibile confronto e tra il grafico della distribuzione empirica e il
grafico della funzione di distribuzione della distribuzione ipotizzata.
• Test statistici
Come e noto, i test delle ipotesi possono essere utilizzati per verificare se
le osservazioni X1, . . . , Xn sono un campione indipendente di una partico-
lare distribuzione di probabilita con funzione di distribuzione F . I due test
piu comuni sono i test Chi–quadro e Kolmogorov–Smirnov e sono adat-
ti al caso che stiamo esaminando anche se essi presentano le loro limita-
zioni intrinseche. Un’esposizione sintetica di questi due test e riportata
nell’Appendice C.
2.2.6 Scelta delle distribuzioni di input in assenza di dati
In alcuni casi, nella pratica, puo accadere che non sia possibile raccogliere dati
sul funzionamente del sistema che si vuole studiare perche esso e ancora in fase
di progettazione e quindi non ancora esistente. In questi casi non sono quindi
disponibili dati da utilizzare per selezionare una distribuzione di input ad una
simulazione e quindi non sono applicabili le tecniche viste fino ad ora. Senza
entrare nei dettagli, osserviamo solamente che sara necessario far ricorso a proce-
dure euristiche che si basano sulla natura del sistema, sul ricorso a persone esperte
di sistemi della tipologia di interesse, sulle limitazioni fisiche o convenzionali del
processo in esame.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
GENERAZIONE DI OSSERVAZIONI CASUALI 143
2.3 GENERAZIONE DI OSSERVAZIONI CASUALI
Una volta determinate le distribuzioni di input, la simulazione dovra generare du-
rante ogni esecuzione osservazioni casuali di variabili aleatorie distribuite secondo
particolari distribuzioni di probabilita. Ad esempio, nel simulare un sistema di
code M/M/1 si avra bisogno di generare i tempi di interarrivo in accordo alla di-
stribuzione esponenziale e cosı anche per i tempi di servizio. Un esempio banale
di cio e stato mostrato alla pagina 127, dove nell’effettuare una semplice esem-
plificazione di una simulazione abbiamo avuto bisogno delle due liste di valori
generati casualmente dalle distribuzioni corrispondenti.
Ogni metodo per generare osservazioni casuali da distribuzioni fissate utilizza va-
riabili indipendenti identicamente distribuite secondo la distribuzione uniforme in
[0, 1) (che indichiamo con U(0, 1)), nel senso che costruisce le osservazioni casuali
desiderate a partire da numeri casuali generati uniformemente in [0, 1) attraverso
opportune trasformazioni. Quindi preliminarmente analizziamo brevemente nel
prossimo paragrafo la generazione di numeri casuali con distribuzione uniforme
e nel paragrafo successivo studieremo come generare osservazioni casuali secondo
distribuzioni fissate a partire dalla distribuzione U(0, 1).
2.3.1 Generazione di numeri pseudocasuali con distribuzione uniforme
In generale, per generare successioni di numeri casuali si potrebbero utilizza-
re metodi quali il lancio dei dadi, l’estrazione da urne, ma ovviamente questi
metodi non sono utilizzabili nella simulazione dove e necessario generare lunghe
successioni di numeri casuali in tempi molto brevi. Attualmente per generare suc-
cessioni di numeri casuali si ricorre all’uso del calcolatore e in realta quello che si
fa e la generazione deterministica di successioni di numeri aventi proprieta stati-
stiche che approssimano molto bene quelle di successioni veramente casuali e che
ad un’analisi statistica risultano indistinguibili da successioni di numeri casuali. Numeri
pseudo-
casuali
I numeri determinati con questa procedura si chiamano numeri pseudocasuali.
Per generare numeri pseudocasuali con distribuzione uniforme esistono diversi
metodi; i piu utilizzati sono i generatori congruenziali lineari. In questi metodi Generatori
con-
gruenziali
lineari
una successione di numeri interi Zi viene definita dalla seguente formula ricorsiva
Zi+1 = aZi + c ( mod m )
dove a si chiama moltiplicatore e c viene detto incremento1. Il termine Z0 si
chiama seme. Il generatore si dice moltiplicativo se c = 0. Quindi vengono
generati al piu m numeri interi Zi distinti con 0 ≤ Zi ≤ m − 1. Per ottenere
numeri casuali Ui in [0, 1) e sufficiente definire Ui = Zi/m.
1La notazione “ mod m ” indica la congruenza modulo m, ovvero il resto della divisione per m
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
144 SIMULAZIONE
Esempio 2.3.1 Vediamo un esempio di generatore moltiplicativo (c = 0). Prendiamo
a = 3, Z0 = 3 e m = 7. Si ottiene
Z1 = 9 ( mod 7 ) = 2
Z2 = 6 ( mod 7 ) = 6
Z3 = 18 ( mod 7 ) = 4
Z4 = 12 ( mod 7 ) = 5
Z5 = 15 ( mod 7 ) = 1
Z6 = 3 ( mod 7 ) = 3
Z7 = 9 ( mod 7 ) = 2
Naturalmente la successione e periodica di periodo al piu pari ad m. Se un
generatore ha periodo pari ad m, ovvero il periodo massimo, si dice che il gene-Periodo
pieno ratore ha periodo pieno e, in questo caso, ogni scelta del seme Z0 portera alla
generazione dell’intero ciclo di valori da 0 a m− 1. Se invece un generatore non
ha periodo pieno allora la lunghezza del ciclo puo dipendere dal particolare valore
del seme Z0. E importante avere periodo lungo o, ancora meglio, periodo pieno in
modo che vengono generati tutti gli interi tra 0 e m−1 ed inoltre essi appariranno
esattamente una volta ogni ciclo e questo contribuisce all’uniformita delle Ui.
Il teorema che segue riporta le condizioni che devono essere soddisfatte dai
parametri m, a e c affinche il generatore abbia periodo pieno.
Teorema 2.3.1 Un generatore congruenziale lineare ha periodo pieno se e
solo se sono soddisfatte le seguenti condizioni:
i) m e c sono primi tra loro;
ii) se q e un numero primo che divide m, allora q divide a− 1;
iii) se 4 divide m, allora 4 divide a− 1.
E obiezione comune a tutti i generatori pseudo–random il fatto che le Zi non
sono realmente casuali, ma ogni Zi e completamente determinata dai quattro
parametri m, a, c e Z0. Tuttavia, un’attenta scelta di questi parametri induce
le Zi ad assumere valori tali che le Ui corrispondenti siano indipendenti identica-
mente distribuite secondo la distribuzione uniforme in [0, 1). Un’altra obiezione
riguarda il fatto che, ovviamente ogni numero reale nell’intervallo [0, 1) dovra
avere la stessa probabilita di essere generato, mentre le Ui assumono solamente
valori razionali. A questo inconveniente si puo ovviare scegliendo m molto grande
(m ≥ 109) in modo che i numeri generati Ui costituiscono un sottoinsieme denso
dell’intervallo [0, 1).
Nei linguaggi di programmazione sono di solito disponibili generatori di numeri
pseudocasuali ed esistono metodi statistici per valutarne la qualita.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
GENERAZIONE DI OSSERVAZIONI CASUALI 145
Tra tutti i generatori di numeri pseudocasuali basati su metodi congruenziali
moltiplicativi, nel caso di computer a 32 bit, il piu utilizzato e il generatore di
Learmouth–Lewis dato da c = 0, a = 75 e m = 231 − 1.
2.3.2 Generazione di osservazioni casuali da una distribuzione di probabilita
Esaminiamo ora come, a partire da numeri casuali uniformemente distribuiti in
[0, 1) sia possibile generare osservazioni da una fissata distribuzione di probabilita.
Per questo scopo sono state introdotte molte tecniche. Ci limitiamo nel seguito
a considerarne due tra le piu utilizzate: il metodo della trasformazione inversa e
il metodo dell’accettazione–reiezione.
Metodo della trasformazione inversa
E un metodo per generare osservazioni da una distribuzione di probabilita. Il
metodo si basa sul seguente risultato teorico [Ross, 2003a] valido nel caso di
distribuzioni continue:
Proposizione 2.3.2 Sia U una variabile aleatoria uniforme in [0, 1). Allora
per ogni funzione di distribuzione continua F , la variabile aleatoria
X = F−1(U)
ha funzione di distribuzione F .
Dimostrazione: Sia FX la funzione di distribuzione della variabile aleatoria X.
Quindi, per ogni y risulta
FX(y) = P (X ≤ y) = P(F−1(U) ≤ y
).
Poiche F (x) e una funzione di distribuzione, essa e monotona crescente e quindi
si ha che F−1(U) ≤ y se e solo se U ≤ F (y). Quindi, poiche U ha distribuzione
uniforme, si ha
FX(y) = P (U ≤ F (y)) = F (y),
ovvero F e la funzione di distribuzione della X.
Quindi, sulla base di questo risultato, data una distribuzione di probabilita
con funzione di distribuzione F , a partire dalla distribuzione uniforme in [0, 1)
possiamo costruire una variabile aleatoria la cui funzione di distribuzione e F .
Osservazione 2.3.3 Nel caso in cui la funzione F non e strettamente monotona,
non e possibile definire la sua inversa; in questi casi, si puo utilizzare la funzione
pseudoinversa data da
F−1(x) = inf {y ∈ IR | x ≤ F (y)} .
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
146 SIMULAZIONE
Esempio 2.3.4 Supponiamo di voler costruire una successione di numeri pseudocasuali
come osservazioni dalla distribuzione esponenziale, ovvero con funzione di distribuzione
F (x) = 1 − e−λx. Innanzitutto determiniamo F−1: da u = F (x) = 1 − e−λx si ricava
x = −1/λ ln(1− u), ovvero
F−1(u) = − 1
λln(1− u).
Quindi se U e una variabile aleatoria uniformemente distribuita in [0, 1),
X = F−1(U) = − 1
λln(1− U) (2.3.1)
e una variabile aleatoria con distribuzione esponenziale con media 1/λ. Quindi, da-
ta una successione di numeri pseudocasuali con distribuzione uniforme in [0, 1), dal-
la (2.3.1) possiamo ottenere una successione di numeri pseudocasuali con distribuzione
esponenziale.
Esempio 2.3.5 Utilizzando quanto ricavato nel precedente Esempio 2.3.4 si puo otte-
nere la generazione di osservazioni casuali dalla distribuzione di Erlang. Infatti sappiamo
che la somma di k variabili aleatorie indipendenti identicamente distribuite secondo la
distribuzione esponenziale, ciascuna con media 1/(kµ) ha distribuzione di Erlang di para-
metro k e media 1/µ. Quindi avendo una successione di numeri uniformemente distribuiti
in [0, 1), u1, . . . , uk, le osservazioni dalla distribuzione di Erlang possono essere ottenute
da
x =
k∑i=1
ln(1− ui)−kµ
che e equivalente a
x = − 1
kµln
[k∏i=1
(1− ui)
].
Osservazione 2.3.6 E importante osservare che se una variabile aleatoria U ha
distribuzione uniforme in [0, 1), anche 1 − U ha distribuzione uniforme in [0, 1)
e quindi nella (2.3.1) si puo sostituire nell’argomento del logaritmo (1 − U) con
U . Tuttavia, come vedremo nel seguito nel paragrafo 2.5.1 sulle tecniche per la
riduzione della varianza, questo cambiamento potrebbe indurre un cambiamento
nella correlazione delle variabili X generate.
Il metodo della trasformazione inversa puo essere esteso ed utilizzato anche nel
caso di distribuzioni discrete, ovvero quando si assume che la variabile X sia una
variabile aleatoria discreta. In questo caso, naturalmente si ha
F (x) = P (X ≤ x) =∑xi≤x
p(xi),
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
GENERAZIONE DI OSSERVAZIONI CASUALI 147
dove p(xi) = P (X = xi).
Supponiamo quindi che X assuma i valori x1, x2, . . . e supponiamo che essi siano
ordinati, ovvero x1 < x2 < · · · . Data una variabile U uniformemente distribuita
in [0, 1) si definisce la variabile X nel seguente modo: si determina il piu piccolo
intero positivo k tale che U ≤ F (xk) e si pone X = xk. Dobbiamo ora dimostra-
re che effettivamente la X cosı generata e quella desiderata, ovvero che risulta
P (X = xi) = p(xi) per ogni i. Infatti si ha:
• per i = 1 risulta X = x1 se e solo se U ≤ F (x1), ma F (x1) = p(x1) perche
le xi sono ordinate. Ora, poiche la U e uniformemente distribuita in [0, 1),
si ha P (X = x1) = P (U ≤ F (x1)) = F (x1) = p(x1)
• per i ≥ 2 risulta X = xi se e solo se F (xi−1) < U ≤ F (xi) per come e scelto
i. Inoltre, poche la U e uniformemente distribuita in [0, 1) si ha
P (X = xi) = P (F (xi−1) < U ≤ F (xi)) = F (xi)− F (xi−1) = p(xi)
Il metodo della trasformazione inversa nel caso discreto ha una giustificazione
molto intuitiva: si divide l’intervallo [0, 1) in sottointervalli contigui di ampiezza
p(x1), p(x2), . . . e si assegna X a seconda del fatto che questi intervalli contengano
la U che e stata generata.
Esercizio 2.3.7 Si vuole generare una successione di osservazioni casuali da una distri-
buzione triangolare avente funzione di distribuzione
F (x) =
0 x ≤ 0x2
20 < x ≤ 1
1− (2− x)2
21 < x ≤ 2
1 x > 2
con il metodo della trasformazione inversa. Generare i primi quattro termini della suc-
cessione utilizzando come generatore di numeri pseudocasuali in [0, 1) un generatore
congruenziale lineare moltiplicativo con Z0 = 3, m = 5 e a = 3.
Esercizio 2.3.8 Si vuole generare una successione di osservazioni casuali da una distri-
buzione di probabilita avente per funzione di distribuzione
F (x) = exp
{− exp
(−x− µ
σ
)}con il metodo della trasformazione inversa. Generare i primi quattro termini della suc-
cessione (ponendo σ = 2 e µ = 1) utilizzando come generatore dei numero pseudocasuali
in [0, 1) un generatore congruenziale lineare moltiplicativo con Z0 = 5, m = 7 e a = 3.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
148 SIMULAZIONE
Esercizio 2.3.9 Si vuole generare una successione di osservazioni casuali dalla distribu-
zione di probabilita di Cauchy avente per funzione di distribuzione
F (x) =1
2+
1
πarctan
x− st
−∞ < x <∞
con il metodo della trasformazione inversa. Generare i primi quattro termini della suc-
cessione (ponendo s = 0 e t = 1) utilizzando come generatore dei numero pseudocasuali
in [0, 1) un generatore congruenziale lineare moltiplicativo con Z0 = 5, m = 7 e a = 3.
Metodo dell’accettazione–reiezione
Il metodo della trasformazione inversa e basato sul calcolo della trasformazione
inversa F−1 che non sempre puo essere calcolata o comunque non in maniera
efficiente. Per questa ragione sono stati sviluppati altri metodi fra i quali il
metodo che esaminiamo in questo paragrafo detto “acceptance–rejection” o anche
“metodo del rigetto”.
Consideriamo il caso continuo e supponiamo di voler generare osservazioni ca-
suali da una distribuzione di probabilita avente funzione di distribuzione F e
densita di probabilita f (il caso discreto si tratta in maniera del tutto analoga).
Supponiamo di disporre di un metodo per generare osservazioni casuali da una
variabile aleatoria Y avente per densita di probabilita una funzione g(x). Il me-
todo accettazione–reiezione utilizza queste osservazioni per generare osservazioni
casuali dalla distribuzione di probabilita avente per densita di probabilita la fun-
zione f(x). In particolare, si generano osservazioni casuali della variabile aleatoria
Y e poi si accettano o si rifiutano queste osservazioni come osservazioni casuali
delle distribuzione con densita data dalla f con una probabilita proporzionale al
rapporto f(Y )/g(Y ). Piu in dettaglio, sia c una costante tale che
f(y)
g(y)≤ c, per ogni y.
Il metodo si puo schematizzare nel seguente modo:
Passo 1: si genera un osservazione casuale della variabile aleatoria Y ;
Passo 2: si genera un numero casuale U dalla distribuzione uniforme in [0, 1)
(indipendente da Y );
Passo 3: se U ≤ f(Y )
cg(Y )allora si pone X = Y e STOP
altrimenti si torna al Passo 1.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
GENERAZIONE DI OSSERVAZIONI CASUALI 149
Le generazione viene quindi effettuata mediante un algoritmo iterativo che ad
ogni passo genera un coppia di osservazioni casuali (Y, U) e si arresta quando e
soddisfatta la diseguaglianza U ≤ f(Y )
cg(Y ), accettando il valore di Y come osserva-
zione casuale della X. Alla fine dell’algoritmo si avra quindi che la variabile X ha
la distribuzione condizionata di Y dato l’evento “accettazione”, ovvero l’evento{U ≤ f(Y )
cg(Y )
}.
Si puo dimostrare che la variabile aleatoria X le cui osservazioni casuali sono
generate con il metodo accettazione–reiezione ha la densita di probabilita voluta,
ovvero vale la seguente proposizione.
Proposizione 2.3.10 La variabile aleatoria X generata con il metodo
accettazione–reiezione ha densita di probabilita f .
Dimostrazione: Osserviamo innanzitutto che, poiche la posizione X = Y si ha
quando avviene l’“accettazione”, allora risulta
P (X ≤ x) = P
(Y ≤ x/U ≤ f(Y )
cg(Y )
).
Inoltre, per l’indipendenza di Y da U si ha
P
(U ≤ f(Y )
cg(Y )/Y = y
)= P
(U ≤ f(y)
cg(y)
)=
f(y)
cg(y)(2.3.2)
Utilizzando la (2.3.2) si ottiene
P
(U ≤ f(Y )
cg(Y )
)=
∫ ∞−∞
P
(U ≤ f(Y )
cg(Y )/Y = y
)g(y)dy =
∫ ∞−∞
f(y)
cg(y)g(y)dy =
1
c
perche f e una densita di probabilita. Utilizzando le relazioni gia ottenute,
possiamo quindi calcolare
P (X ≤ x) = P
(Y ≤ x/U ≤ f(Y )
cg(Y )
)=P(Y ≤ x , U ≤ f(Y )
cg(Y )
)P(U ≤ f(Y )
cg(Y )
) =
= c
∫ x
−∞P
(U ≤ f(Y )
cg(Y )/Y = y
)g(y)dy
= c
∫ x
−∞
f(y)
cg(y)g(y)dy =
∫ x
−∞f(y)dy = F (x).
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
150 SIMULAZIONE
Esempio 2.3.11 Applichiamo il metodo dell’accettazione–reiezione per generare osser-
vazioni casuali da una variabile aleatoria avente densita di probabilita f(x) = 20x(1−x)3,
0 < x < 1 (si tratta della distribuzione beta con parametri 2 e 4). Innanzitutto dobbiamo
scegliere la funzione g(x): in questo una scelta di una funzione molto semplice e g(x) = 1,
0 < x < 1, ovvero si sceglie la distribuzione uniforme. Ora si deve determinare una co-
stante c tale che risultif(x)/g(x) ≤ c per ogni x. Naturalmente, in questo caso si puo
determinare facilmente c massimizzando la funzione f(x)/g(x); poiche la derivata prima
di f(x)/g(x) si annulla in x? = 1/4 e la derivata seconda in questo punto e negativa, il
punto x? e un punto di massimo per la funzione f(x)/g(x) e risulta
f(x)
g(x)≤ 135
64= c,
e quindi il test diventa U ≤ 256/27x(1− x)3.
Quindi il metodo in questo caso puo essere cosı realizzato:
Passo 1: si genera un’osservazione casuale U1 e dalla distribuzione uniforme in [0, 1);
Passo 2: si genera un’osservazione casuale U2 dalla distribuzione uniforme in [0, 1);
Passo 3: se U2 ≤256
27U1(1− U1)3 allora si pone X = U1 e STOP, altrimenti si torna al
Passo 1.
Esistono inoltre tecniche speciali per generare osservazioni da variabili aleatorie
distribuite secondo distribuzioni particolari per le quali si rimanda alla letteratura
dedicata.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
PROGETTAZIONE DI UNA SIMULAZIONE E ANALISI DELL’OUTPUT 151
2.4 PROGETTAZIONE DI UNA SIMULAZIONE E ANALISI DELL’OUTPUT
In questo paragrafo analizziamo uno degli aspetti fondamentali di una simula-
zione, ovvero la progettazione della simulazione stessa e l’analisi statistica dei
risultati. L’errore tipico che si commette nell’effettuare una simulazione e quello
di eseguire un singolo “run” di lunghezza arbitraria e di prendere i risultati che
si ottengono come stime delle caratteristiche del vero modello. Poiche in una
simulazione si utilizzano osservazioni casuali da una o piu distribuzioni di proba-
bilita, queste stime sono solamente realizzazioni particolari di variabili aleatorie
che possono anche avere varianza molto grande. Questo implica che, in un parti-
colare run di una simulazione, queste stime possono anche differire di molto dalle
corrispondenti reali caratteristiche del modello.
Affinche i risultati di uno studio effettuato attraverso la simulazione abbiano senso
e necessario l’uso di tecniche statistiche per progettare e analizzare gli esperimen-
ti di una simulazione. Nel fare cio si manifesta subito il problema derivante dal
fatto che i processi di output della simulazione sono, in generale, autocorrelati
e non stazionari e questo rende inapplicabili direttamente le tecniche statisti-
che classiche che sono invece basate su osservazioni indipendenti, identicamente
distribuite.
Un altro problema esistente nell’ottenere stime precise delle caratteristiche del
vero modello e dato dal tempo di calcolo necessario per raccogliere la quan-
tita necessaria di dati di output per effettuare un’analisi statistica efficace, che
potrebbe essere anche molto elevato.
2.4.1 Dati di output
Siano Y1, Y2, . . . i dati di output di un singolo run di una simulazione; ciascuna
Yi puo essere vista come una variabile aleatoria e quindi la collezione di variabili
aleatorie {Yi, i = 1, 2, . . .} e un processo stocastico. Ad esempio, in un sistema a
coda, le Yi possono rappresentare il tempo di attesa in coda dell’i-esimo utente.
In generale, le variabile aleatorie Yi potrebbero non essere ne indipendenti, ne
identicamente distribuite e quindi nell’analisi di questi dati non possono essere
applicati direttamente i metodi di analisi statistica. Per ovviare a questo in-
conveniente, si effettuano piu repliche della simulazione, ciascuna di lunghezza
m e si basa l’analisi sulle varie repliche. Formalmente, siano y11, y12, . . . y1m la
realizzazione delle variabili aleatorie Y1, . . . , Ym ottenute con la prima replica.
Nella seconda replica si avranno differenti realizzazioni delle variabili aleatorie
Y1, . . . , Ym; siano esse y21, y22, . . . y2m. In generale, supponiamo di aver effettuato
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
152 SIMULAZIONE
n repliche indipendenti di lunghezza m, ottenendo cosı
y11 y12 · · · y1m
y21 y12 · · · y2m...
......
yn1 yn2 · · · ynm
(2.4.1)
Le realizzazioni di una stessa replica non sono indipendenti, identicamente di-
stribuite, ma se per ogni i = 1, . . . ,m consideriamo le osservazioni y1i, y2i, . . . yni,
ovvero l’i-esima colonna della (2.4.1) allora esse costituiscono osservazioni in-
dipendenti, identicamente distribuite della variabile aleatoria Yi. Quindi l’a-
nalisi statistica e applicabile alle osservazioni y1i, y2i, . . . yni, per ogni fissato
i = 1, . . . , n.
2.4.2 Transitorio e stato stazionario
Un altro problema importante nella progettazione e nell’analisi dell’output di
una simulazione e dato dalla scelta delle condizioni iniziali e dal fatto che in un
sistema si deve distinguere lo stato stazionario e il transitorio.
Sia {Yi, i = 1, 2, . . .} il processo stocastico di output e, per ogni i = 1, 2, . . ., sia
Fi(y | I) = P (Yi ≤ y | I) ,
dove y ∈ IR e I rappresenta le condizioni iniziali, ovvero le condizioni alle quali
inizia la simulazione al tempo 0; Fi(y | I) e quindi la probabilita condizionata
che l’evento {Yi ≤ y} accada, date le condizioni iniziali I. Ad esempio, in un
impianto di produzione, I potrebbe rappresentare il numero dei job presenti al
tempo 0, oppure in un sistema di trasporto potrebbe indicare la posizione dei
veicoli al tempo 0.
La Fi(y | I) e la distribuzione transitoria del processo di output al tempo (di-Distribuzio-
ne transi-
toria
screto) i con condizioni iniziali I. In generale, la Fi(y | I) e diversa per differenti
valori di i e per ogni insieme di condizioni iniziali I.
Per y e I fissati, la successione {Fi(y | I)}i puo essere convergente per i → ∞.
Se per ogni y e I esiste il limite
limi→∞
Fi(y | I) = F (y), (2.4.2)
allora F (y) e la distribuzione stazionaria del processo di output {Yi, i = 1, 2, . . .}.Distribuzio-
ne stazio-
naria
Lo stato stazionario si raggiunge, in teoria, al limite per i→∞, ma, nella pratica,
molto spesso esiste un indice temporale finito k tale che le distribuzioni da quel
punto in poi rimarrano approssimativamente coincidenti. Quando questo accade,
si assume che il sistema e nello stato stazionario a partire dal tempo k.
Osservazione 2.4.1 Ovviamente il raggiungimento dello stato stazionario non
equivale a dire che le variabili aleatorie Yk, Yk+1, . . . assumeranno tutte lo stesso
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
PROGETTAZIONE DI UNA SIMULAZIONE E ANALISI DELL’OUTPUT 153
valore in una particolare replica, ma solamente che hanno (approssimativamente)
la stessa distribuzione
Osservazione 2.4.2 La distribuzione F (y) non dipende dalle condizioni iniziali
I, ma puo dipendere dalle condizioni iniziali la rapidita di convergenza delle
distribuzioni Fi(y | I) a F (y).
2.4.3 Tipi di simulazioni
Una distinzione fondamentale deve essere compiuta tra simulazioni per le quali
esiste un evento “naturale” il cui accadimento indica la conclusione della simula-
zione stessa e simulazioni in cui questo non accade. Ad esempio, se si e interessati
a simulare un sistema di code che rappresenta un ufficio postale che apre alle 8.30
e chiude alle 14.00 e finisce di operare dopo che e stato servito l’ultimo cliente
entrato prima della chiusura, e chiaramente presente l’evento dato dal fatto che
il sistema e vuoto dopo le ore 14.00 che identifica la fine dell’operativita del siste-
ma stesso e quindi della simulazione corrispondente. Un’altro esempio potrebbe
essere dato da un’industria navale che ha un contratto per produrre 30 navi in
24 mesi; si vorrano simulare diverse configurazioni in cui puo operare l’industria
per studiare quale permetta di evadere l’ordine entro la data fissata; in questo
caso l’evento e la fabbricazione di 30 navi. In altri casi, invece un tale evento
non esiste; si pensi, ad esempio, ad un sistema di produzione continuo in cui si
vogliono valutare misure di prestazione a regime.
Formalmente distinguiamo i due casi parlando, nel primo caso, di simulazione con Simulazione
con o senza
termina-
zione
terminazione (finita), quando esiste un evento naturale che specifica la lunghezza
di ciascuna replica. In questo caso le condizioni iniziali influenzano le misure di
prestazione del sistema e verra naturalmente effettuata un’analisi del transitorio.
Nel secondo caso si parla di simulazione senza terminazione nel senso che non
esiste un evento che specifica la lunghezza di una replica; in questo caso si e
di solito interessati al comportamento a regime del sistema, ovvero all’analisi
dello stato stazionario e questo puo comportare che per stimare una misura di
prestazione possano essere necessari run molto lunghi della simulazione.
Si deve anche notare come in alcuni casi reali lo stato stazionario potrebbe non
essere raggiunto perche le caratteristiche del sistema cambiano nel tempo; tutta-
via, il corrispondente modello di simulazione potrebbe raggiungerlo, poiche nel
costruire il modello si e assunto che le caratteristiche del modello non cambia-
no nel tempo. In questi casi, quando ci sono cambiamenti delle caratteristiche
del sistema e necessario inserire questi cambiamenti nel modello ed effettuare di
nuovo lo studio.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
154 SIMULAZIONE
2.4.4 Simulazioni con terminazione: analisi del transitorio
Supponiamo di effettuare n repliche indipendenti di una simulazione con ter-
minazione dove ciascuna replica ha le medesime condizioni iniziali e finisce con
l’evento che determina il termine della replica. Per semplicita, assumiamo che ci
sia una singola misura di prestazione da considerare rappresentata dalla variabile
aleatoria X e sia Xi la realizzazione di questa variabile aleatoria ottenuta nella
i-esima replica. Per quanto visto nel paragrafo 2.4.1, le Xi sono variabili aleatorie
indipendenti, identicamente distribuite e questo permette di determinare una sti-
ma della media µ = E(X) e il relativo intervallo di confidenza dalle X1, . . . , Xn.
Infatti si puo definire lo stimatore (corretto) della media µ dato da
Xn =
n∑i=1
Xi
n
e lo stimatore (corretto) della varianza dato da
s2n =
n∑i=1
(Xi − Xn)2
n− 1.
Inoltre, dato α ∈ (0, 1), un intervallo di confidenza approssimato al 100(1− α)%
e dato da
Xn ± tn−1,1−α2
√s2n
n,
dove tn−1,1−α2
e il punto critico (1− α2 ) della distribuzione t student a n−1 gradi
di liberta. Indicando con
δ(n, α) = tn−1,1−α2
√s2n
n(2.4.3)
la semiampiezza dell’intervallo di confidenza, si puo scrivere l’intervallo di confi-
denza nella forma [Xn − δ(n, α) , Xn + δ(n, α)
].
Se n fosse fissato “a priori” (fixed–sample–size), non ci sarebbe nessuna possibilita
di controllare la precisione della stima Xn, ovvero l’ampiezza dell’intervallo di
confidenza. Quindi e necessario considerare la possibilita di decidere il numero
delle repliche da effettuare in modo da raggiungere una precisione desiderata.
Naturalmente esistono due modi di misurare l’errore che si commette:
errore assoluto∣∣Xn − µ
∣∣errore relativo
∣∣Xn − µ∣∣
|µ|
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
PROGETTAZIONE DI UNA SIMULAZIONE E ANALISI DELL’OUTPUT 155
Iniziamo ad esaminare il caso dell’errore assoluto. In genere si richiede di effet-
tuare un numero di repliche n in modo tale che risulti Errore
assolutoP(|Xn − µ| ≤ β
)≥ 1− α, α ∈ (0, 1), β > 0. (2.4.4)
La (2.4.4) afferma che l’errore assoluto commesso e pari al piu a β con probabi-
lita pari ad almeno 1 − α. Quindi, se, ad esempio, costruiamo 100 intervalli di
confidenza al 90% utilizzando il criterio di arresto sul numero delle repliche che
abbiamo adottato, ci aspettiamo che l’errore assoluto∣∣Xn − µ
∣∣ sia al piu pari a
β in circa 90 dei 100 casi, mentre nei rimanti 10 casi l’errore assoluto potrebbe
essere maggiore di β. Esistono regole che determinano il numero minimo n?a di
repliche da effettuare affinche la (2.4.4) sia soddisfatta; una delle piu semplici e
la seguente:
n?a = min {n ≥ 2 | δ(n, α) ≤ β} .
Per quanto riguarda l’errore relativo, si richiede che il numero di repliche n da Errore
relativoeffettuare sia tale che risulti
P
(|Xn − µ||µ|
≤ γ)≥ 1− α, α ∈ (0, 1), γ > 0. (2.4.5)
Anche in questo caso esistono regole empiriche che forniscono il valore del numero
minimo n?r di repliche da effettuare affinche la (2.4.5) sia soddisfatta; una di queste
e la seguente:
n?r = min
{n ≥ 10 | δ(n, α)
|Xn|≤ γ
1 + γ
}.
Il termine γ/(1 + γ) deriva dal fatto che nel calcolare l’errore relativo stiamo
utilizzando |Xn|, ovvero una stima di |µ| e non il valore esatto |µ|. Infatti,
supponiamo di aver effettuato un numero di repliche sufficienti a garantire che
risultiδ(n, α)
|Xn|≤ L. (2.4.6)
Allora si ha
1− α ≈ P(|Xn − µ| ≤ δ(n, α)
)= P
(|Xn − µ||Xn|
≤ δ(n, α)
|Xn|
)≤
≤ P
(|Xn − µ||Xn|
≤ L)
= P(|Xn − µ| ≤ L|Xn − µ+ µ|
)≤
≤ P(
(1− L)|Xn − µ| ≤ L|µ|)
= P
(|Xn − µ||µ|
≤ L
1− L
).
Quindi, affiche valga la (2.4.5) si deve avere L/(1−L) = γ, ovvero L = γ/(1+γ).
Per applicare nella pratica quanto ora discusso sono state definite due strategie
generali:
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
156 SIMULAZIONE
• procedura a due fasi: in una prima fase vengono effettuate n0 repliche sulla
base delle quali si calcola s2n0
e la semiampiezza dell’intervallo di confidenza
δ(n0, α). Poi, eventualmente, si effettuano altre repliche fino al raggiungi-
mento della precisione desiderata non ricalcolando la stima della varianza
all’aumentare delle repliche;
• procedura iterativa: si aumenta iterativamente il numero delle repliche fino
a che non si raggiunge il livello di accuratezza desiderato, ricalcolando la
stima della varianza ad ogni nuova replica.
Procedura a due fasi
Con questa procedura e possibible determinare il numero delle repliche necessarie
per stimare il valore atteso µ = E(X) con una precisione fissata. Nella prima fase
vengono effettuate un numero fissato n0 di repliche e sulla base di queste si calcola
la stima della varianza s2n0
e l’intervallo di confidenza per µ. A questo punto, nella
seconda fase, assumendo che la stima della varianza non cambi significativamente
all’aumentare del numero delle repliche, si effettuano, eventualmente, ulteriori
repliche fino ad ottenere la precisione voluta. Considerando l’errore assoluto, il
numero totale delle repliche da effettuare per ottenere un errore assoluto pari al
piu a β e dato da
n?a = min
{i ≥ n0 | ti−1,1−α
2
√s2n0
i≤ β
},
ovvero si incrementa il numero delle repliche aggiuntive i di 1 fino a quando si
ottiene la disuguaglianza verificata. Ovviamente, se n?a > n0, allora si dovranno
effettuare n?a − n0 repliche aggiuntive.
Analogamente, considerando l’errore relativo, supponendo di aver costruito un
intervallo di confidenza per µ basato su un numero fissato n0 di repliche, allora
si puo verificare che il numero totale delle repliche da effettuare per ottenere un
errore relativo pari a γ e
n?r = min
i ≥ n0 |ti−1,1−α
2
√s2n0
i|Xi|
≤ γ
1 + γ
.
Anche in questo caso, se tale valore n?r e maggiore del numero delle repliche
gia effettuate n0, sara necessario effettuare n∗r − n0 repliche aggiuntive fino al
soddisfacimento della disuguaglianza.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
PROGETTAZIONE DI UNA SIMULAZIONE E ANALISI DELL’OUTPUT 157
Procedura iterativa
Quanto visto nel paragrafo precedente presenta l’inconveniente di utilizzare la
stima della varianza basata sulle n0 repliche fissate nel calcolo del numero di
repliche necessarie per ottenere una precisione desiderata. Tuttavia tali stime
potrebbe essere imprecise e quindi questo potrebbe portare ad una scelta del nu-
mero delle repliche troppo grande con notevole spreco di tempo di calcolo, oppure
tale numero potrebbe essere troppo piccolo e quindi, di fatto, non si otterrebbe la
precisione voluta. Infatti, se s2n0
e una sottostima della varianza, allora il numero
di repliche necessarie determinato dalla procedura e piu piccolo del necessario,
mentre se s2n0
e una sovrastima della varianza, allora tale numero e piu grande del
necessario e quindi vengono effettuate repliche inutili. Questo inconveniente puo
essere superato utilizzando una procedura iterativa che ha lo scopo di determi-
nare una stima di µ con errore relativo pari a γ (o assoluto pari a β) e intervallo
di confidenza del 100(1 − α)%. La differenza fondamentale rispetto alla proce-
dura a due fasi vista nel paragrafo precedente sta nel fatto che, ad ogni replica
aggiuntiva eventualmente effettuata, viene ricalcolata la stima della varianza che
invece nel caso precedente rimaneva fissata a s2n0
. Uno schema algoritmico che
permette di ottenere una stima di µ con errore relativo non superiore a γ ovvero
che permette di soddisfare la (2.4.5) e il seguente:
Passo 0: si effettuano n0 repliche della simulazione e si pone n = n0;
Passo 1: si calcolano Xn e δ(n, α) da X1, . . . , Xn;
Passo 2: seδ(n, α)
|Xn|≤ γ
1 + γsi usa Xn come stima di µ e STOP;
altrimenti si effettua un’ulteriore replica, si pone n = n + 1 e si va al
Passo 1.
Da una sperimentazione diretta su un grande numero di modelli e di distribuzioni
di probabilita per le quali sono noti i valori di µ e stato verificato che la procedura
ora descritta fornisce una buona copertura dell’intervallo di confidenza al 90%
utilizzando n0 ≥ 10 e γ ≤ 0.15.
Naturalmente, la stessa procedura puo essere utilizzata considerando l’errore asso-
luto al posto dell’errore relativo; tuttavia, il fatto che l’errore assoluto sia piccolo,
ovviamente e estremamente dipendente dal modello e molto sensibile alla scelta
di β e quindi, in generale, si preferisce una procedura basata sull’errore relativo.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
158 SIMULAZIONE
2.4.5 Simulazioni senza terminazione: analisi dello stato stazionario
Sia {Yj , j = 1, 2, . . .} il processo di output di una simulazione senza terminazione,
e supponiamo che esista il limite (2.4.2). Supponiamo di voler stimare la media
µ = limj→∞
E(Yj) = E(Y )
della variabile aleatoria di interesse Y . Naturalmente, datem osservazione Y1, . . . Ym,
una stima di µ si puo ottenere da
Ym =
m∑j=1
Yj
m. (2.4.7)
Tuttavia, a causa del transitorio iniziale dipendente dalla scelta delle condizioni
iniziali, Ym non e uno stimatore corretto in quanto e possibile che E(Ym) 6= µ per
ogni m.
Questo problema e detto problema del transitorio iniziale o problema dello star-Problema
dello
“startup”
tup. Per superare questo problema si utilizza una tecnica di cancellazione dei dati
iniziali detta anche “warming up” del modello che consiste nel non considerare
nella stima le prime osservazioni che sono quelle piu influenzate dalle condizioni
iniziali. Quindi, invece di utilizzare la stima Ym data dalla (2.4.7), si considera
Y(m,`) =
m∑j=`+1
Yj
m− `,
dove quindi ` e il numero delle osservazioni che vengono scartate. Ovviamente
c’e il problema di come scegliere `, ovvero il periodo di warming up, in modo che
E(Y(m,`)) ≈ µ. Se il valore di ` e troppo piccolo c’e il rischio che la stima risenta
delle condizioni iniziali, mentre se e troppo grande si potrebbe avere uno spreco
di tempo di calcolo.
Sono state proposte diverse regole per determinare `, nessuna delle quali parti-
colarmente soddisfacente. C’e, invece, una procedura basata su un’analisi grafica
che e efficace. Poiche si vuole determinare un valore ` tale che per i > ` risulti
E(Yi) ≈ µ, graficamente questo si traduce nel determinare quando la curva E(Yi),
i = 1, 2, . . . si “stabilizza” intorno al valore µ. Il problema e che, molto spesso,
si ha un andamento influenzato da un’alta varianza e quindi molto oscillante.
Per superare questo inconveniente, e stata introdotta una procedura nota comeProcedura
di Welch procedura di Welch che, prima di effettuare l’analisi grafica prevede un “tratta-
mento” dei dati per ridurre la varianza. Uno schema della procedura di Welch e
il seguente:
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
PROGETTAZIONE DI UNA SIMULAZIONE E ANALISI DELL’OUTPUT 159
Passo 1: Si effettuano n repliche ciascuna di lunghezza m:
Y11 Y12 · · · Y1m
Y21 Y12 · · · Y2m...
......
Yn1 Yn2 · · · Ynm.
Passo 2: Si costruisce la successione Y1, Y2, . . . Ym dove
Yj =
n∑i=1
Yij
n.
Passo 3: Dato k ≤ bm/4c, si sostituisce ciascun termine della successione
Y1, Y2, . . . Ym con
Yj(k) =
k∑h=−k
Yj+h
2k + 1, j = k + 1, . . . ,m− k
j−1∑h=−(j−1)
Yj+h
2j − 1, j = 1, . . . , k
Passo 4: Si sceglie quel valore di ` oltre il quale la successione{Yj(k)
}appare
giunta a convergenza.
Al Passo 2 si effettua la media fra i dati corrispondenti nelle diverse repliche,
costruendo cosı le Yj . Poiche risulta E(Yj) = E(Yj) e V ar(Yj) = V ar(Yj)/n,
si ha che la media e rimasta invariata, mentre la varianza e stata ridotta di un
fattore 1/n.
Al Passo 3, si sostituisce ciascuna Yj ottenuta al passo precedente con la media
tra essa e le Yj immediatamente precedenti e successive ovvero, se j non e troppo
vicino all’inizio della successione, Yj(k) non e altro che la media di 2k + 1 valori
delle Yj centrate alla j-esima. Il valore k che determina quanti termine inclu-
dere nelle media si chiama “time window”. Il valore Yj(k) e chiamato “moving
average”.
Ad esempio, per k = 2 si ha:
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
160 SIMULAZIONE
Y1(2) = Y1
Y2(2) =1
3
(Y1 + Y2 + Y3
)Y3(2) =
1
5
(Y1 + Y2 + Y3 + Y4 + Y5
)Y4(2) =
1
5
(Y2 + Y3 + Y4 + Y5 + Y6
)...
Ym−2(2) =1
5
(Ym−4 + Ym−3 + Ym−2 + Ym−1 + Ym
).
Nella pratica, le linee guida che conviene seguire per la determinazione della
lunghezza del transitorio sono le seguenti: si disegna il grafico delle Yj(k) per
diversi valori di k e si sceglie il piu piccolo valore di k (se esiste) per il quale
il grafico e “sufficientemente smooth”; si utilizza questo grafico per determinare
la lunghezza ` del periodo di warming up. Per quanto riguarda il valore del
time window k, esso deve essere sufficientemente grande da rendere regolare il
grafico delle Yj(k), ma non cosı grande da non permettere l’individuazione del
transitorio.
Per concludere, ribadiamo il fatto che e molto importante scgeliere in modo
appropriato i tre parametri fondamentali di una simulazione, m, n e `:
• per quanto riguarda la lunghezza delle repliche m, essa dovra essere suffi-
cientemente grande tanto da essere piu grande del valore che ci si aspetta
per ` e tale da permettere nella simulazione un numero elevato di occorrenze
di tutti gli eventi, anche quelli poco probabili;
• per quanto riguarda il numero delle repliche n, e opportuno iniziare con
valori di n pari a 5 o 10 per poi aumentare il valore, se necessario;
• la lunghezza ` del periodo di warm-up deve essere opportunamente deter-
minata, eventualmente utilizzando un trattamento dei dati di output per
ridurre la loro varianza.
Stima della media basata sull’approccio “repliche–cancellazioni”
Supponiamo di aver effettuato n repliche di una simulazione, ciascuna di lunghez-
zam (repliche pilota) e di aver determinato il valore della lunghezza del transitorio
`. L’analisi statistica da effettuare per il calcolo della media e dell’intervallo di
confidenza e del tutto analoga a quella gia vista nel paragrafo 2.4.4 nel caso di
simulazioni con terminazione, con l’unica differenza che consiste nell’utilizzare,
in ciascuna replica, solamente le osservazioni successive alla `-esima.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
PROGETTAZIONE DI UNA SIMULAZIONE E ANALISI DELL’OUTPUT 161
Sia n′ il numero delle repliche ed m′ la lunghezza di ogni replica, assumendo che
m′ � ` e dove ` e stato determinato con la procedura di Welch. Per ogni replica
i si definisce una variabile aleatoria
Xi =
m′∑k=`+1
Yik
m′ − `, i = 1, . . . n′.
Per l’indipendenza si ha che le Xi sono indipendenti, identicamente distribuite
con E(Xi) ≈ µ e, analogamente al caso dell’analisi del transitorio, si puo definire
Xn′ =
n′∑i=1
Xi
n′
che e uno stimatore (approssimativamente corretto) della media µ; un intervallo
di confidenza approssimato del 100(1− α)% per µ e
Xn′ ± tn′−1,1−α2
√s2n′
n′.
Nella maggior parte dei casi, per il basso costo del tempo di calcolo, non costitui-
sce un problema l’uso di n repliche pilota per determinare ` e poi l’uso delle sole
ultime m′ − ` osservazioni provenienti da n′ repliche differenti da quelle pilota.
Naturalmente continuano a valere tutte le considerazioni fatte nel paragrafo 2.4.4
per quanto riguarda la possibilita di ottenere una precisione desiderata.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
162 SIMULAZIONE
2.5 TECNICHE PER LA RIDUZIONE DELLA VARIANZA
Come abbiamo avuto modo piu volte di osservare, l’output di una simulazione
necessita di un’analisi statistica approfondita per ottenere una buona interpre-
tazione dei risultati. Alcune volte, simulazioni di grandi dimensioni richiedono
molto tempo di calcolo e nonostante questo, la precisione dei risultati che si ot-
tengono potrebbe non essere del tutto soddisfacente. Naturalmente la qualita
delle stime e tanto migliore quanto piu e piccola la varianza delle variabili alea-
torie di output. Questo significa che, se si riesce a ridurre tale varianza senza
alterare il valore atteso, si puo ottenere una precisione migliore, ovvero interval-
li di confidenza piu piccoli a parita di carico computazionale della simulazione,
oppure si riesce a raggiungere una precisione desiderata con un minor numero di
repliche e/o con repliche di minore lunghezza. A tale scopo, si possono utilizzare
le cosiddette tecniche per la riduzione della varianza. Ovviamente e generalmen-
te impossibile conoscere “a priori” la misura della riduzione della varianza che
dovra essere ottenuta. Esistono diverse tecniche di riduzione della varianza; nel
seguito ne consideriamo due tra le piu utilizzate: l’uso di variabili antitetiche e
il condizionamento.
2.5.1 Variabili antitetiche
L’uso di variabili antitetiche si basa sulla seguente osservazione: sia X una varia-
bile di output di una simulazione e supponiamo di voler stimare il valore atteso
θ = E(X). A tale scopo, supponiamo di aver generato due variabili aleatorie
X1, X2 identicamente distribuite con media θ. Consideriamo la nuova variabile
aleatoria (X1 +X2)/2; ovviamente essa ha la stessa media θ, e la sua varianza e
data da
V ar
(X1 +X2
2
)=
1
4
(V ar(X1) + V ar(X2) + 2Cov(X1, X2)
).
Ora, se X1 e X2 fossero negativamente correlate (ovvero Cov(X1, X2) < 0),
utilizzando (X1 +X2)/2 si avrebbe una riduzione della varianza.
Sulla base di questa semplice osservazione, la tecnica delle variabili antitetiche
consiste nell’effettuare coppie di repliche della simulazione e utilizzare la media
delle coppie delle osservazioni. La correlazione negativa e indotta utilizzando
numeri casuali complementari nella generazione delle osservazioni casuali delle
distribuzione di probabilita.
Come abbiamo visto nel paragrafo 2.3.2, durante una simulazione, per generare
osservazioni casuali secondo una distribuzione di probabilita si utilizza una suc-
cessione di numeri casuali Uk uniformemente distribuiti nell’intervallo [0, 1). La
tecnica delle variabili antitetiche prevede di utilizzare nella prima replica della
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
TECNICHE PER LA RIDUZIONE DELLA VARIANZA 163
coppia numeri casuali Uk per generare le osservazioni di una certa distribuzio-
ne di probabilita, mentre nella seconda replica, per generare osservazioni dalla
medesima distribuzione di probabilita si utilizzano numeri casuali 1−Uk (comple-
mentari). Questo e possibile, in quanto come gia visto nell’Osservazione 2.3.6, se
Uk e distribuita uniformemente in [0, 1), anche 1−Uk e distribuita uniformemente
in [0, 1).
Formalmente, si ha
X1 = h(U1, . . . , Um)
X2 = h(1− U1, . . . , 1− Um),
con Uk variabile aleatoria con distribuzione uniforme in [0, 1) e h un’opportuna
funzione. Si puo dimostrare [Ross, 2003a], che se h e una funzione monoto-
na, allora le X1 e le X2 sono correlate negativamente. Infatti, vale la seguente
proposizione.
Proposizione 2.5.1 Se h(x1, . . . , xm) e una funzione monotona nel-
le x1, . . . , xm, allora, se U1, . . . , Um e un insieme di numeri casuali
uniformemente distribuiti in [0, 1) allora
Cov (h(U1, . . . , Um), h(1− U1, . . . , 1− Um)) ≤ 0.
Si procede, quindi, in questo modo: si effettuano n coppie di repliche di una
simulazione ottenendo le osservazioni(X
(1)1 , X
(2)1
),(X
(1)2 , X
(2)2
), . . . ,
(X(1)n , X(2)
n
)dove la X
(1)j e ottenuta dalla prima replica della j−esima coppia di repliche
utilizzando Uj e dove la X(2)j e ottenuta dalla seconda replica della j−esima
coppia di repliche utilizzando 1− Uj .Quindi, per ogni j = 1, . . . , n, si hanno coppie di osservazioni
(X
(1)j , X
(2)j
)e
quindi si puo definire
Xj =X
(1)j +X
(2)j
2, j = 1, . . . n,
e calcolare la stima della media
Xn =
n∑j=1
Xj
n.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
164 SIMULAZIONE
Si ha
V ar(Xn) =V ar(Xj)
n=
1
nV ar
(X
(1)j +X
(2)j
2
)=
=V ar
(X
(1)j
)+ V ar
(X
(2)j
)+ 2Cov
(X
(1)j , X
(2)j
)4n
.
Quindi, inducendo una correlazione negativa tra X(1)j e X
(2)j , la varianza V ar(Xn)
viene ridotta. Naturalmente e importante che le repliche siano sincronizzate,
ovvero se nel generare il j−esimo elemento della prima replica si utilizza Uj ,
allora nel generare il corrispondente elemento j−esimo della seconda replica deve
essere utilizzato 1− Uj .
2.5.2 Condizionamento
Un’altra tecnica di riduzione della varianza e basata sul condizionamento. Questa
tecnica e basata su alcuni noti risultati di Calcolo delle probabilita riguardanti
il valore atteso condizionato e la varianza condizionata. In particolare, si fa
riferimento ai risultati riportati nella seguente proposizione.
Proposizione 2.5.2 Siano X e Y due variabili aleatorie. Allora valgono le
seguenti relazioni:
i) E(X) = E (E(X/Y ))
ii) V ar(X) = E (V ar(X/Y )) + V ar (E(X/Y )).
Dimostrazione: [Questa dimostrazione non fa parte del programma di esame]
Dimostriamo queste relazioni nel caso in cui le variabili X e Y siano discrete (il
caso di variabili continue si tratta in maniera analoga).
Dato un qualsiasi valore y di Y si puo calcolare E(X | Y = y) e risulta
E(X | Y = y) =∑x
xP (X = x | Y = y) =∑x
xP (X = x, Y = y)
P (Y = y),
da cui
E (E(X | Y )) =∑y
E(X | Y = y)P (Y = y)
=∑y
∑x
xP (X = x, Y = y)
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
TECNICHE PER LA RIDUZIONE DELLA VARIANZA 165
=∑x
x∑y
P (X = x, Y = y)
=∑x
xP (X = x) = E(X)
e questo dimostra la i). Per dimostrare la ii) si consideri che vale
V ar(X/Y ) = E(
(X − E(X/Y ))2 /Y)
= E(X2/Y )− (E(X/Y ))2 . (2.5.1)
Prendendo il valore atteso ad entrambi i membri della (2.5.1) si ha
E (V ar(X/Y )) = E(E(X2/Y )
)− E
((E(X/Y ))2
)= E(X2)− E
((E(X/Y ))2
). (2.5.2)
Inoltre vale
V ar(E(X/Y )) = E((E(X/Y ))2
)− (E(E(X/Y )))2
= E((E(X/Y ))2
)− (E(X))2 . (2.5.3)
Sommando la (2.5.2) e la (2.5.3) si ha
E (V ar(X/Y )) + V ar(E(X/Y )) = E(X2)− (E(X))2,
ovvero la ii).
Dalla ii) segue facilmente che poiche la varianza e sempre non negativa, cia-
scun termine della ii) deve essere maggiore o uguale a zero e quindi risulta
E (V ar(X/Y )) ≥ 0, da cui si ha
V ar(X) ≥ V ar(E(X/Y )). (2.5.4)
Questo ci permette di affermare che se X e una variabile di output di una si-
mulazione e supponiamo di voler stimare il valore atteso θ = E(X), per la i), si
ha
θ = E(X) = E (E(X/Y )) ,
ovvero E(X/Y ) e uno stimatore corretto della media e, per la (2.5.4), si ha che
la sua varianza e inferiore. Quindi puo essere piu conveniente utilizzare E(X/Y )
come stima di θ.
In sintesi, lo spirito di questa tecnica e quello di sfruttare qualche proprieta
speciale di un modello di simulazione (ove esista), sostituendo una stima di una
quantita con il suo valore analitico esatto. Rimuovendo questa fonte di variabilita
si dovrebbe avere una variabile aleatoria di output piu “stabile”, anche se non c’e
alcuna garanzia che questo sia verificato.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
166 SIMULAZIONE
2.6 L’APPROCCIO SIMULAZIONE OTTIMIZZAZIONE
(“SIMULATION–BASED OPTIMIZATION”)
La Simulazione e l’Ottimizzazione possono interagire e integrarsi quando si vuole
determinare qual e la combinazione di dati di input di una simulazione che ot-
timizza uno o piu indici di prestazione. In questo caso la funzione obiettivo del
problema di ottimizzazione non e disponibile in forma analitica, o comunque non
calcolabile con procedure deterministiche, ma il suo valore in corrispondenza di
diversi valori delle variabili e ottenibile solamente attraverso il risultato di una
simulazione. Potrebbero presentarsi casi in cui oltre la funzione obiettivo anche
i vincoli non siano disponibili e la loro verifica puo essere effettuata solamente
attraverso una simulazione. Infine, anche il valore di alcune variabili (uncontrolla-
ble variables) del problema di ottimizzazione potrebbe essere ottenuto solamente
attraverso una simulazione, ovvero dipendere dai dati di input della simulazione
(controllable variables).
Ovviamente lo scopo del processo di ottimizzazione e ben diverso da quello di
analisi di scenario proprio della simulazione: infatti, come gia evidenziato nell’in-
troduzione allo studio della Simulazione, nel primo caso si deve decidere quale
configurazione del sistema e la migliore tra tutte le ammissibili (approccio “what–
best”); nel secondo caso le configurazioni alternative del sistema sono date e ov-
viamente in numero finito (approccio “what–if”). Si osservi che molto spesso
questo numero non puo essere elevato a causa dell’alto costo computazionale dei
run di simulazione.
Per evidenziare l’importanza e l’attualita della simulazione ottimizzazione, ripor-
tiamo di seguito alcune righe della prefazione del testo “Handbook on simulation
optimization” di Michael C. Fu di recente pubblicazione [Fu, 2014]:
Arguably, the two most powerful operations research/management science
(OR/MS) techniques are simulation and optimization. [ ... ] Optimization
dates back many centuries and is generally considered the older of the two
siblings. Both approaches were propelled forward by the advent of the digital
computer over half a century ago, leading up to the present golden age
when both routinely address complex large-scale realworld problems and
both are implemented in a large variety of computer software packages.
However, combining the two techniques is a more recent development, and
software effectively integrating the two is relatively limited; thus, simulation
optimization remains an exciting and fertile area of research.
Risulta molto evidente come attualmente l’approccio simulazione ottimizzazio-
ne sia considerato uno strumento molto efficace per lo studio di problemi reali
complessi e di grandi dimensioni che sempre piu spesso si presentano in numerosi
settori applicativi.
Nella sua forma piu semplice un problema di Ottimizzazione e Simulazione puo
essere formulato come segue: supponiamo di avere una misura di prestazione ot-
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
L’APPROCCIO SIMULAZIONE OTTIMIZZAZIONE (“SIMULATION–BASED OPTIMIZATION”) 167
tenuta in output da una simulazione e indichiamo tale misura con f . Il valore
che f assume dipende dal valore dei dati di input che denotiamo con x1, . . . , xn.
Per il problema di ottimizzazione questi ultimi saranno le variabili di decisione.
Poiche f e ottenuta dall’output di una simulazione, essa e, in generale, una varia-
bile aleatoria. Lo scopo sara pertanto quello di minimizzare/massimizzare il suo
valore atteso E (f(x1, . . . , xn)) su tutte le combinazioni possibili di x1, . . . , xn.
Naturalmente ci possono essere vincoli di vario tipo sulle variabili di decisione (i
dati di input) che limitano le combinazioni possibili (ammissibili) da considerare.
Sottolineamo ancora il fatto che la funzione obiettivo non e calcolata analitica-
mente o comunque in maniera deterministica, ma ottenuta come output di una
simulazione. Sara quindi necessario effettuare un certo numero di repliche indi-
pendenti della simulazione e utilizzare per la f il valore medio dei diversi valori
ottenuti in queste repliche (ovviamente in corrispondenza dei medesimi dati di
input).
In una forma piu generale un problema di simulazione ottimizzazione, puo essere
formalmente scritto nel seguente modo:
min E [f(x, z, ω)]
s.t. E [g(x, z, ω)] ≥ 0
xl ≤ x ≤ xuzl ≤ z ≤ zu
x ∈ IRn, z ∈ Zn,
dove f e funzione a valori reali, g e h sono funzioni a valori vettoriali. La funzione
f viene valutata attraverso la simulazione di particolari istanze delle variabili
continue x e delle variabili discrete z (input del modello) e della realizzazione
di varibili aleatorie, ovvero il vettore ω, che puo essere anch’esso funzione degli
input x e z. In pratica, l’introduzione del vettore ω serve per evidenziare l’effetto
stocastico della simulazione. Analogamente, i vincoli rappresentati dalla funzione
g sono valutati attraverso la simulazione. Inoltre ci possono essere vincoli che non
dipendono da variabili aleatorie (rappresentati dalla funzione h) e vincoli di tipo
box sulle variabili.
Uno schema complessivo dell’implementazione di una procedure di simulazione
ottimizzazione puo essere cosı sintetizzato: viene simulata una sequenza di confi-
gurazioni del sistema, ciascuna corrispondente ad una particolare scelta dei dati
di input. Ad ogni step della sequenza un algoritmo di ottimizzazione effettua
una ricerca nello spazio dei possibili dati di input (regione ammissibile) fino alla
determinazione della configurazione ottima. Quindi l’algoritmo di ottimizzazione
sceglie una configurazione iniziale e il simulatore effettua le repliche della simu-
lazione. Il risultato della simulazione e inviato all’algoritmo di ottimizzazione
e viene verificato un criterio di arresto. Se questo non e verificato, l’algoritmo
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
168 SIMULAZIONE
di ottimizzazione determina una nuova configurazione da simulare e il processo
continua iterativamente fino al soddisfacimento del criterio di arresto.
Per una trattazione completa della “simulation–based optimization” si rimanda al
gia citato [Fu, 2014]. Riportiamo di seguito in estrema sintesi alcune osservazioni
importanti:
• il tempo totale di un’esecuzione dipende dal numero di configurazioni del
sistema che sono state esaminate (numero di iterazioni dell’algoritmo di
ottimizzazione), ma anche da tempo necessario per simulare ciascuna di
queste configurazioni;
• in caso di varianza elevata, sara necessario un numero di repliche adeguato
per ottenere l’accuratezza desiderata;
• esiste la possibilita di effettuare implementazioni parallele.
Molti simulatori includono “optimization packages”’ (spesso basati su algorim-
ti euristici) gia integrati con il simulatore e che prevedono una visulizzazione
di vario tipo durante le iterazioni dell’algoritmo. Un esempio e costituito da
tool OptQuestr for ARENAr disponibile all’interno del simulatore ARENAr (cfr.
paragrafo 2.7).
Disponendo di un proprio codice che implementa l’algoritmo di ottimizzazione
che si vuole utilizzare, sara ovviamente necessario creare un’interfaccia tra il
simulatore e l’algoritmo di ottimizzazione e questo spesso richiede tecnicalita che
possono essere non da banali da realizzare, soprattutto per la sincronizzazione
dei due tool. A titolo di esempio, il simulatore ARENAr permette di utilizzare
Microsoft Visual Basic (VBA) per realizzare tale interfaccia.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
SOFTWARE DI SIMULAZIONE (SIMULATORI) 169
2.7 SOFTWARE DI SIMULAZIONE (SIMULATORI)
Nel paragrafo 2.1.4 abbiamo fatto cenno ad alcuni pacchetti software di simula-
zionei. L’uso di tali simulatori e oggi diventato assai diffuso all’interno di diverse
realta della vita reale in ragione della loro generalita e flessibilita di uso. Faremo
riferimento a simulatori ad eventi discreti “general purpose”, ovvero non destinati
alla creazione di modelli di simulazione specifici bensı capaci di creare modelli di
simulazione di sistemi complessi anche molto diversi fra di loro. Tali simulatori
hanno ormai raggiunto un ottimo livello di adattabilita e facilita di uso. Sono di
solito basati su un’interfaccia grafica che permette di costruire un modello con
il semplice “drag and drop” di moduli predefiniti all’interno di appositi riquadri
di lavoro, ignorando l’esistenza del linguaggio nel quale il simulatore realizzera
il modello stesso. Inoltre, sono dotati di animazioni grafice anche tridimensio-
nali che permettono di osservare il funzionamento del sistema implementato nel
modello di simulazione. Tali animazioni sono molto utili se si vuole mostrare
tale funzionamento a persone che sono poco propense ad osservare tabelle spes-
so molto lunghe e piene di valori numerici poco significativi ad un occhio poco
esperto.
Abbiamo scelto i due seguenti simulatori:
• ARENAr http://www.arenasimulation.com
Prodotto dalla Rockwell Automation, e un simulatore oggi largamente uti-
lizzato, che attraverso un’interfaccia grafica consente sia di realizzare un
modello di simulazione, sia di effettuare i diversi run di una simulazione e
di analizzare i risultati ottenuti.
• SIMIOTM
http://www.simio.com
Prodotto dalla SIMIO LCC, anch’esso con una efficace interfaccia grafica
e basato su “oggetti intelligenti” che una volta costruiti, possono essere
riutilizzati in progetti diversi.
SIMIO e un simulatore di piu recente generazione e rappresenta un passaggio
dal paradigma della simulazione “process oriented” a quello della simulazione
“objects oriented”. Entrambi i simulatori sono dotati di strumenti di animazione e
visualizzazione 3D che permettono di osservare bene il funzionamento del sistema
reale.
Gli studenti posso scaricare dai rispettivi siti le versioni trial/studente delle
release piu recenti di entrambi.
Per ogni dettaglio sull’uso di ARENA si fa riferimento a [Kelton et al., 2004],
[Altiok, Melamed, 2007], [Rossetti, 2010] e alla User’s Guide.
Per ogni dettaglio sull’uso di SIMIO si fa riferimento a [Kelton et al., 2011],
[Joines, Roberts, 2012] e alla User’s Guide.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
170 SIMULAZIONE
2.7.1 ARENA
Si riporta una descrizione molto sintetica dei moduli che costituiscono i Basic
process di ARENA, ovvero i moduli: CREATE, PROCESS, DISPOSE, DECIDE,
ASSIGN, BATCH e SEPARATE.
• CREATE : genera le entita e le immette nel sistema. La generazione puo
avvenire secondo uno schema fissato o specificando la distribuzione di pro-
babilita dei tempi di interarrivo.
• PROCESS : rappresenta qualsiasi tipo di processamento di entita. Permet-
te anche di considerare l’utilizzazione di risorse (e dei relativi vincoli di
disponibilita) secondo tre schemi: size–delay, seize–delay–release e delay–
release.
• DISPOSE : espelle le entita dal sistema. Permette anche di specificare se
registrare le statistiche dell’ingresso delle entita.
• DECIDE : permette di controllare l’istradamento delle entita sulla base di
condizioni oppure sulla base di probabilita in corrispondenza per ciascuna
uscita.
• ASSIGN : esegue l’assegnazione di variabili e attributi. L’operazione di
assegnazione e effettuata quando un’entita attraversa il modulo.
• BATCH : permette di raggruppare entita. La costruzione del “batch” e
controllata dal numero delle entita o da attributi. Le entita che entrano nel
modulo vengono raggruppate e si avra una sola entita (gruppo) in uscita
dal modulo.
• SEPARATE : ha due funzioni: produrre una o piu copie di una entita,
oppure separare un gruppo di entita precedentemente raggruppate.
2.7.2 Modelli in ARENA
Si riportano di seguito i testi di alcuni esercizi svolti durante le esercitazioni. I
relativi modelli implementati con ARENA sono disponibili sul sito web del corso.
Si tratta di esempi molto introduttivi destinati solamente ad illustrare l’uso dei
moduli presenti nei Basic process di ARENA.
MODELLO 1. Un modello che rappresenti un sistema a coda M/M/1 con
media dei tempi di interarrivo pari a 60 minuti e media dei tempi di servizio pari
a 50 minuti.
MODELLO 2. Ad un istituto di credito arrivano richieste di erogazione di
mutui. La distribuzione dei tempi di interarrivo di tali tempi e esponenziale con
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
SOFTWARE DI SIMULAZIONE (SIMULATORI) 171
media 2 ore. La richiesta viene esaminata da un impiegato il cui tempo per
completare la pratica e distrbuito secondo una ditribuzione triangolare con valori
(1, 1.75, 3) ore. In genere, nell’88% dei casi la pratica va a buon fine, ovvero il
mutuo e concesso. Nel rimanente 12% dei casi, la risposta e negtiva.
MODELLO 3. Un’impresa per la trasformazione di semilavorati deve proces-
sare l’arrivo di 100 semilavorati. Questi dovranno passare attraverso 2 ulteriori
processi di lavorazione: detti Machining Center ed Inspection Station. Dopo tali
processi verra effettuato il controllo di qualita sul prodotto finito valutandone il
grado di finitezza. I dati relativi al problema sono:
• al primo processo produttivo i semilavorati arrivano uno alla volta con
intertempo di arrivo esponenziale (media pari ad 1 minuto);
• i tempi di processamento (espressi in minuti) al Machining Center sono
uniformemente distribuiti nell’intervallo [0.65, 0.70];
• i tempi di processamento (espressi in minuti) al Inspection Station sono
uniformemente distribuiti nell’intervallo [0.75, 0.80];
• con probabilita del 50% i pezzi controllati sono “perfect” ed escono dal
sistema attraverso un proprio canale;
• con probabilita del 20% i pezzi controllati sono “very good” ed escono dal
sistema attraverso un proprio canale;
• con probabilita del 5% i pezzi controllati sono “good” ed escono dal sistema
attraverso un proprio canale;
• con probabilita del 25% i pezzi controllati rientrano nel sistema e vengono
rilavorati sui due processi.
MODELLO 4. In un laboratorio di analisi di un ospedale il centro prelievi e
aperto dal lunedı al sabato dalle 7.30 alle 11.00. Gli utenti arrivano al centro con
intertempi di arrivo distribuiti esponenzialmente con media 7 minuti. Appena
arrivato al centro prelievi un utente deve ritirare un numero da una macchinetta
distributrice e poi deve attendere di essere chiamato allo sportello per l’accetta-
zione dove presentera la richiesta di analisi del proprio medico curante; gli utenti
sono chiamati in ordine di numero crescente e le operazioni per l’accettazione
sono eseguite in tempi distribuiti uniformemente tra 1 e 1.5 minuti. Dopo l’ac-
cettazione l’utente si reca all’ufficio cassa se deve pagare il ticket oppure, se e
esente da tale pagamento, si reca direttamente all’ambulatorio per il prelievo. Il
pagamento del ticket richiede tempi distribuiti uniformemente tra 0.5 e 1 minu-
to, mentre i tempi per effettuare i prelievi presso l’ambulatorio sono distribuiti
esponenzialmente con media 5 minuti.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
172 SIMULAZIONE
MODELLO 5. Una catena di montaggio costruisce apparecchi assemblando
due parti. Le due parti arrivanno alla catena dall’esterno da due ingressi diversi
(I1 e I2 ). Le singole parti arrivano alla catena di montaggio secondo due diverse
distribuzione di probabilita dei tempi di interarrivo: dall’ingresso I1 secondo una
distribuzione uniforme tra 0.4 e 0.7 ore, dall’ingresso I2 secondo una distribuzione
normale a media 0.5 ore e deviazione standard 0.2. Le singole parti, prima dell’as-
semblaggio subiscono un pretrattamento; i tempi necessari per questa operazione
sono distributi secondo la distribuzione normale a media 0.1 ore e deviazione
standard 0.03 per le parti che sono arrivate dall’ingresso I1 e a media 0.15 ore e
deviazione standard 0.04 per le parti arrivate dall’ingresso I2. L’assemblaggio e
effettuato in un tempo costante pari a 0.3 ore. Dopo l’assemblaggio c’e una lavo-
razione finale che richiede un tempo distribuito secondo la distribuzione normale
a media 0.6 ore e deviazione standard 0.13.
MODELLO 6. In una videoteca entrano clienti con tempi di interarrivo di-
stribuiti esponenzialmente con media 5 minuti. Una volta entrati i clienti scelgono
casualmente tra i tre diversi tipi di film (drammatico, commedia, azione). Il
tempo impiegato nella scelta ha distribuzione triangolare con valore minimo 1
minuto, massimo 10 minuti e valore piu probabile 5 minuti. Una volta effettuata
la scelta il cliente puo decidere se acquistare anche una pizza oppure dirigersi
direttamente alla cassa. Nel 25% dei casi la pizza viene acquistata e questo
richiede un tempo distribuito secondo una distribuzione triangolare con valore
minimo 0.5 minuti, massimo 2 minuti e valore piu probabile 1 minuto. Il tempo
necessario al pagamento alla cassa (con o senza la pizza) e anch’esso distribuito
secondo una distribuzione triangolare con valore minimo 1 minuto, massimo 4
minuti e valore piu probabile 2 minuti.
MODELLO 7. Un Computer consta di 1 CPU, 3 terminali, 1 unita Hard Disk
(HD) e 1 unita a nastro (Tape). In ciascun terminale siede un utente il quale
invia jobs alla CPU con intertempi distribuiti esponenzialmente (valore atteso
pari a 100 secondi). I jobs si accumulano in una coda (FIFO) prima di essere
processati dalla CPU. Ciascun job processato dalla CPU la tiene occupata per
un tempo distribuito esponenzialmente (valore atteso pari a 1 secondo). Un job
che lascia la CPU puo:
• con probabilita 20% ritornare al terminale che l’ha inviato;
• con probabilita 72% andare nella coda (FIFO) dell’ Hard Disk;
• con probabilita 8% andare nella coda (FIFO) dell’unita nastro.
Ciascun job processato dall’Hard Disk necessita di un tempo distribuito espo-
nenzialmente (valore atteso pari a 1.39 secondi), dopodiche torna alla coda del-
la CPU. Analogamente ciascun job processato dall’unita a nastro necessita di
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
SOFTWARE DI SIMULAZIONE (SIMULATORI) 173
un tempo distribuito esponenzialmente (valore atteso pari a 12.50 secondi), poi
torna alla coda della CPU. Tutti i job prodotti dagli operatori sono tra loro
statisticamente indipendenti.
MODELLO 8. Un’azienda di vendite per corrispondenza riceve richieste di or-
dini con distribuzione esponenziale a media 9 minuti. Due operatori telefonici
accettano gli ordini. Per entrambi il tempo di servizio e distribuito secondo la
distribuzione triangolare con valore minimo 3 minuti, massimo 11 minuti e valore
piu probabile 7 minuti il primo operatore e valore minimo 2 minuti, massimo
11 minuti e valore piu probabile 8 minuti il secondo. Una volta terminata la
chiamata l’ordine viene passato all’addetto del reparto consegne che provvede a
preparare l’ordine e a caricarlo su un furgone. Tale operazione richiede un tem-
po distribuito secondo la distribuzione triangolare con valore minimo 5 minuti,
massimo 20 minuti e valore piu probabile 10 minuti. Prima di partire per le
consegne, il furgone attende che vi siano caricate le merci corrispondenti a 10
ordini. I tempi per le consegne sono distribuiti secondo la distribuzione triango-
lare con valore minimo 20 minuti, massimo 80 minuti e valore piu probabile 40
minuti. Per ciascun ordine consegnato il fattorino fornisce la notifica di avvenuta
consegna all’azienda che provvede ad inserire il cliente nella lista delle consegne
effettuate. Tale operazione avviene con un tempo distribuito secondo la distri-
buzione triangolare con valore minimo 1 minuto, massimo 4 minuti e valore piu
probabile 2 minuti. L’ultima operazione e l’emissione della fattura da inviare al
cliente che avviene con un tempo distribuito secondo la distribuzione triangolare
con valore minimo 0.8 minuti, massimo 2 minuti e valore piu probabile 1 minuto.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
174 SIMULAZIONE
2.8 RIFERIMENTI DEL CAPITOLO 2
La trattazione di questo capitolo e basata sul testo [Law, Kelton, 2000] che rappre-
senta un testo completo su tutti gli aspetti della simulazione e quindi utile anche
per approfondimenti. Per aspetti piu teorici si puo consultare il capitolo 11 del
testo [Ross, 2003a]. Fra i numerosi ulteriori testi sulla simulazione citiamo, ad
esempio, [Ross, 2002] e [Banks et al., 2001]. Il testo [Banks, 1998] riporta una
collezione di articoli che forniscono un’adeguata panoramica sulle metodologie e
le aree di applicazione della simulazione.
Per quanto riguarda il software di simulazione ARENA trattato nel paragrafo 2.7
si fa riferimento al testo [Kelton et al., 2004].
Segnaliamo inoltre il sito Web dell’Informs College on Simulation
http://www.informs-cs.org
e il sito dei proceedings delle conferenze annuali “Winter Simulation Conference”
http://www.informs-cs.org/wscpapers.html
nel quale si possono trovare molti articoli e survey aggiornate su diversi aspetti
della simulazione e su applicazioni della simulazione.
Banks, J. (1998). Handbook on Simulation. Principles, Methodology, Ad-
vances, Applications, and Practice. Wiley and Sons, New York.
Banks, J., J.S.Carson, Nelson, B., Nicol, D. (2001). Discrete–event
system simulation. Prentice Hall.
Fu, M. C. (2014). Handbook on Simulation Optimziation. Springer.
Law, A., Kelton, W. (2000). Simulation modeling and analysis. McGraw
– Hill, New York, third edition.
Kelton, W., Sadowski, R., Sturrock, D. (2004). Simulation with
Arena. McGraw – Hill, New York, third edition.
Ross, S. (2002). Simulation. Academic Press, San Diego, third edition.
Ross, S. (2003a). Introduction to probability models. Academic Press, San
Diego.
Ross, S. (2003b). Probabilita e statistica. Apogeo, Milano.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
ESERCIZI SULLA SIMULAZIONE 175
2.9 ESERCIZI SULLA SIMULAZIONE
Esercizio 2.9.1 Applicare il metodo della trasformazione inversa per generare
cinque osservazioni casuali dalla distribuzione esponenziale con media pari ad 1,
a partire dai seguenti cinque numeri casuali generati dalla distribuzione uniforme
in [0, 1)
0.495 0.335 0.791 0.469 0.279
Esercizio 2.9.2 Si vuole effettuare una simulazione di un sistema di code M/M/3
in cui la frequenza media di arrivo e pari a 20 e la velocita di servizio e pari a
10, essendo interessati a stimare il numero medio di utenti presenti nel sistema.
Supponiamo che siano state effettuate 5 repliche indipendenti dalle quali si sono
ottenuti i seguenti risultati:
1a replica: 3.28
2a replica: 2.57
3a replica: 3.01
4a replica: 2.93
5a replica: 2.55
Calcolare la stima della media e l’intervallo di confidenza al 95%. Facendo rife-
rimento a tale intervallo di confidenza, trascurando il problema dello start up e
supponendo di voler ottenere un errore assoluto sulla stima della media inferio-
re a 0.2, determinare se le cinque repliche effettuate sono sufficienti ad ottenere
tale precisione. In caso negativo, applicando la procedura iterativa, determinare
il numero delle repliche necessarie per ottenere la precisione voluta sapendo che
nelle eventuali successive cinque repliche si ottengono i seguenti valori:
6a replica: 2.69
7a replica: 2.61
8a replica: 2.85
9a replica: 2.91
10a replica: 2.88
Esercizio 2.9.3 Applicare la procedura di Welch, con un valore di time window
k = 2, al seguente output di una simulazione
1a replica: 3.72 2.90 3.32 2.50 2.12 3.90 3.10 2.60
2a replica: 2.37 3.80 4.00 3.75 2.90 3.12 2.30 2.70
3a replica: 3.12 2.05 2.90 2.55 2.95 3.30 3.00 2.80
4a replica: 2.98 3.90 3.54 3.10 2.78 2.10 1.90 2.35
5a replica: 2.50 1.90 3.10 2.60 2.20 3.00 2.70 2.85
Dei valori Yj(2) ottenuti dalla procedura, (trascurando il problema dello startup)
determinare media e t intervallo di confidenza al 95%.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
Appendice C:Elementi di statistica
In questo capitolo, sono riportati in breve alcuni elementi di Statistica di uso
comune all’interno della Simulazione che lo studente dovrebbe aver gia acquisito
nel suo curriculum di studi attraverso corsi specifici. Tali importanti aspetti non
sono stati trattati esplicitamente, ma considerati come prerequisiti gia in possesso
degli studenti e sono qui riportati per completezza di informazione. Pertanto
questo capitolo non fa esplicitamente parte del programma di esame.
C.1 ELEMENTI DI STATISTICA INFERENZIALE
In questo paragrafo verranno illustrati alcuni elementi di Statistica che sono essen-
ziali per procedere alla costruzione di un modello di simulazione e per effettuare
correttamente una simulazione. Si tratta essenzialmente della stima dei parame-
tri e dei test statistici, elementi chiave per l’analisi dell’input e dell’output di una
simulazione.
Infatti, per effettuare una simulazione di un sistema che presenta elementi stoca-
stici e necessario specificare le distribuzioni di probabilita che regolano i processi
che caratterizzano il sistema stesso. Se e possibile raccogliere dati reali (osser-
vazioni) sulle variabili aleatorie di interesse, essi possono essere utilizzati per
determinare queste distribuzioni facendo uso di tecniche di inferenza statistica
(analisi dell’input). Una volta stabilite tali distribuzioni, la simulazione procede
generando valori casuali da queste distribuzioni, ovvero, durante ogni esecuzione,
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
178 ELEMENTI DI STATISTICA
la simulazione genera osservazioni casuali di variabili aleatorie distribuite secondo
particolari distribuzioni di probabilita. Oltre che per progettare una simulazione,
e necessario l’uso di tecniche statistiche anche per interpretare i risultati ottenuti
da una simulazione (analisi dell’output).
La trattazione degli argomenti riportati in questo paragrafo sara sintetica e consi-
dera esclusivamente quegli elementi che risulteranno utili all’interno dello studio
della simulazione. Per ogni trattamento piu esteso e approfondito si rimanda ai
testi specifici come, ad esempio [Ross, 2003b].
C.2 STATISTICHE CAMPIONARIE E LORO DISTRIBUZIONE
In generale, nello studio di un fenomeno riguardante un insieme di elementi (po-
polazione) che presenta caratteristiche aleatorie, molto spesso si dispone solo di
informazioni su una parte di essi (campione) e si vogliono dedurre proprieta ge-
nerali riguardanti l’intera popolazione. L’inferenza statistica si occupa di questa
problematica e riveste un importante strumento di analisi.
Solitamente viene fatta l’assunzione che esiste una distribuzione di probabilita
della popolazione nel senso che se da essa vengono estratti casualmente alcu-
ni elementi, ad essi sono associate variabili aleatorie indipendenti identicamente
distribuite secondo tale distribuzione. In questo senso, un insieme di variabili
aleatorie X1, . . . , Xn di variabili aleatorie indipendenti tutte con la stessa distri-
buzione si dice campione di questa distribuzione. L’interesse principale risiede
nella possibilita di dedurre caratteristiche della distribuzione non nota sulla base
dei dati a disposizione. Naturalmente ci sono casi in cui della distribuzione della
popolazione non si conosce nulla (se non il fatto che essa e discreta o continua),
mentre in altri casi la distribuzione e nota ma non sono noti alcuni suoi parame-
tri. Esamineremo entrambi i casi, introducendo, innanzitutto alcuni elementi che
riguardano le statistiche campionarie e la loro distribuzione.
Dato un campione X1, . . . , Xn estratto da una popolazione, ogni funzione delle
osservazioni campionarie e chiamata statistica campionaria e i valori ottenuti da
una qualsiasi funzione dei soli valori osservati sono chiamati statistiche. I due
principali esempi di statistiche sono la media campionaria e la varianza campio-
naria.
C.2.1 Media campionaria e varianza campionaria
Sia dato un campione X1, . . . , Xn estratto da una popolazione, ovvero le Xi sono
variabili aleatorie indipendenti identicamente distribuite, e sia µ e σ2 rispetti-
vamente la loro media e la loro varianza (ovvero la media e la varianza della
popolazione).
La media campionaria e data daMedia cam-
pionaria
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
STATISTICHE CAMPIONARIE E LORO DISTRIBUZIONE 179
Xn =1
n
n∑i=1
Xi.
Xn e una variabile aleatoria funzione delle Xi e si verifica facilmente che risulta
E(Xn) = µ e V ar(Xn) =σ2
n.
La varianza campionaria e data da Varianza
campiona-
rias2n =
1
n− 1
n∑i=1
(Xi − Xn
)2e si verifica facilmente che risulta E(s2
n) = σ2.
C.2.2 Distribuzione (approssimata) della media campionaria
La distribuzione della media campionaria puo essere determinata grazie al Teo-
rema del limite centrale. Siano date le osservazioni X1, . . . , Xn estratte da una
distribuzione di probabilita. Le Xi sono variabili aleatorie indipendenti identica-
mente distribuite e sia µ e σ2 rispettivamente la loro media e la loro varianza.
Definendo la variabile aleatoria
Zn =
X1 + · · ·+Xn
n− µ
σ√n
,
e la sua funzione di distribuzione Fn(z) = P (Zn ≤ z), il Teorema del Limite
Centrale afferma che
limn→∞
Fn(z) =1√2π
∫ z
−∞e−
y2
2 dy, (C.2.1)
ovvero Fn(z) converge alla funzione di distribuzione della distribuzione Normale
standard. Questo risultato permette di ottenere l’importante proprieta riportata
nel seguente teorema riguardante la media campionaria Xn.
Proposizione C.2.1 Sia X1, . . . , Xn un campione estratto da una distribu-
zione di probabilita a media µ e varianza σ2. Allora, per n sufficientemente
grande,Xn − µσ√n
(C.2.2)
e una variabile aleatoria distribuita approssimativamente secondo la distribu-
zione Normale standard.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
180 ELEMENTI DI STATISTICA
Questo risultato si ricava immediatamente dal Teorema del Limite Centrale, os-
servando che la (C.2.1) puo essere interpretata nel seguente modo: per n sufficien-
temente grande, la variabile aleatoria Zn e distribuita approssimativamente come
una variabile Normale standard, indipendentemente dalla distribuzione delle Xi.
Ricordando che il valore atteso di Xn e µ e che la sua deviazione standard e
pari a σ/√n, la Proposizione C.2.1 afferma che se si normalizza Xn sottraendo
la sua media e dividendo per la sua deviazione standard, si ottiene una variabile
aleatoria che e approssimativamente distribuita secondo una Normale standard
purche n sia sufficientemente grande.
La problematica che nasce nell’utilizzare questo risultato sta nel fatto che non e
noto quanto deve essere grande n affinche l’approssimazione sia buona; natural-
mente questo dipende dalla distribuzione in questione. Osservazioni sperimentali
hanno portato a formulare la regola empirica comunemente adottata secondo la
quale quando n > 30, si ha in genere una buona approssimazione, qualsiasi sia la
distribuzione della popolazione considerata.
C.2.3 Distribuzioni delle statistiche di popolazioni normali
Assumiamo ora che la distribuzione della popolazione sia Normale a media µ e
varianza σ2 e sia X1, . . . , Xn un campione estratto da tale popolazione. In questo
caso ovviamente si ha che la variabile aleatoria (C.2.2) e una variabile Normale
standard per qualsiasi valore di n, ovvero vale il seguente risultato.
Proposizione C.2.2 Sia X1, . . . , Xn un campione estratto da una distribu-
zione di probabilita Normale a media µ e varianza σ2. Allora
Xn − µσ√n
(C.2.3)
e una variabile aleatoria distribuita secondo la distribuzione Normale standard.
Si osservi che avendo assunto che le Xi sono variabili aleatorie normali, il ri-
sultato ora enunciato non vale piu in senso approssimato, come nel caso della
Proposizionea C.2.1, ma e un risultato esatto che vale per qualsiasi valore di n
(non necessariamente grande).
Si puo inoltre dimostrare che vale il seguente importante risultato:
Teorema C.2.1 Sia X1, . . . , Xn un campione estratto da una distribuzione
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
STIMA DI PARAMETRI 181
Normale a media µ e varianza σ2. Allora
i) Xn e s2n sono variabili aleatorie indipendenti;
ii) Xn e una variabile aleatoria normale con media µ e varianza σ2/n;
iii) (n − 1)s2n
σ2e una variabile aleatoria distribuita secondo la distribuzione
Chi-quadro a n− 1 gradi di liberta (χ2n−1).
Questo teorema ha un’importante conseguenza che enunciamo nella seguente
proposizione.
Proposizione C.2.3 Sia X1, . . . , Xn un campione estratto da una distribu-
zione di probabilita Normale a media µ. Allora
Xn − µsn√n
(C.2.4)
e una variabile aleatoria distribuita secondo la distribuzione t di Student con
n− 1 gradi di liberta (tn−1).
Confrontando la (C.2.3) e la (C.2.4) si nota che nella Proposizione C.2.3 la varia-
bile aleatoria considerata differisce da quella considerata nella Proposizione C.2.2
per la presenza di sn al posto di σ. Quindi possiamo riepilogare i risultati fino ad
ora ottenuti dicendo che, data una popolazione Normale, se si normalizza la me-
dia campionaria Xn sottraendo la sua media µ e dividendo per la sua deviazione
standard σ/√n, si ottiene una variabile aleatoria Normale standard; se invece si
divide per sn/√n, si ottiene una variabile aleatoria con distribuzione t di Student
con n− 1 gradi di liberta.
C.3 STIMA DI PARAMETRI
Supponiamo ora che la popolazione sia distribuita secondo una distribuzione di
probabilita nota, ma caratterizzata da uno o piu parametri incogniti. Siamo
in questo caso interessati a determinare tali parametri incogniti sulla base di un
campione X1, . . . , Xn. Si tratta di un problema di stima di parametri che consiste
nel determinare, sulla base del campione X1, . . . , Xn, un valore per ciascuno dei
parametri in modo che essi costituiscano la migliore approssimazione dei para-
metri incogniti. Esistono diversi metodi di stima che non sono altro che tecniche
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
182 ELEMENTI DI STATISTICA
per ricavare statistiche e che sono detti stimatori. Uno stimatore e quindi unaStimatori e
stime funzione h(X1, . . . , Xn) delle osservazioni campionarie e il valore che tale funzione
assume in corrispondenza di una particolare realizzazione del campione e detto
stima. Se θ e un parametro incognito, si indichera con θ la stima di θ.
In alcuni casi si determina un unico valore θ come migliore approssimazione pos-
sibile del parametro θ e tale valore viene detto stima puntuale. In altri casi, puoStima
puntuale essere preferibile calcolare due valori dello stimatore ovvero θ1 = h1(X1, . . . , Xn)
e θ2 = h2(X1, . . . , Xn) che definiscono un intervallo [θ1, θ2] tale che, in un cam-
pionamento ripetuto, il valore incognito θ apparterra all’intervallo in una deter-
minata percentuale di casi che e detta confidenza dell’intervallo. In questo casoStima per
intervalli si parla di stima per intevalli.
C.3.1 Proprieta degli stimatori
Continuando ad indicare con θ un parametro incognito, con θ una sua stima
e con h = h(X1, . . . , Xn) uno stimatore, riportiamo in questo paragrafo alcune
definizioni standard.
Definizione C.3.1 Si definisce valore dell’errore di campionamento la diffe-
renza θ − θ. Si chiama distorsione di uno stimatore h la differenza E(h)− θ.
L’errore quadratico medio dello stimatore h e dato da EQM(h) = E(h− θ)2.
Definizione C.3.2 Uno stimatore h = h(X1, . . . , Xn) si dice stimatore cor-
retto del parametro θ se risulta E(h) = θ. Se invece si ha E(h) 6= θ si dice che
h e uno stimatore distorto per θ. Uno stimatore h = h(X1, . . . , Xn) si dice
stimatore efficiente del parametro θ se
i) E(h) = θ
ii) V ar(h) ≤ V ar(h1) per ogni h1 stimatore corretto di θ.
C.3.2 Stima di media e varianza
Supponiamo di avere un campione casuale X1, . . . , Xn e di voler stimare stimare
la media µ e la varianza σ2 della popolazione mediante questo campione, ovvero
supponiamo che X1, . . . , Xn siano variabili aleatorie indipendenti identicamente
distribuite con E(Xi) = µ e V ar(Xi) = σ2, i = 1, . . . , n con µ e σ2 non note.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
STIMA DI PARAMETRI 183
Uno stimatore corretto per la media µ e dato dalla media campionaria
µ = Xn =1
n
n∑i=1
Xi,
in quanto risulta E(Xn) = µ. Uno stimatore corretto per la varianza σ2 e dato
dalla varianza campionaria
σ2 = s2n =
1
n− 1
n∑i=1
(Xi − Xn
)2,
in quanto risulta E(s2n) = σ2.
C.3.3 Metodi di stima
Esistono diversi metodi per stimare i parametri incogniti di una distribuzione.
Riporteremo brevemente nel seguito lo stimatore di massima verosimiglianza Stimatore
di massima
verosimi-
glianza
(Maximum Likelihood Estimator – MLE) che e molto utilizzato e per il quale
rimandiamo alla letteratura specifica per una trattazione completa. Riportiamo
di seguito una breve descrizione.
Date n osservazioni X1, . . . , Xn, assumiamo che esse siano ottenute da una distri-
buzione di probabilita continua avente densita fθ(x), dove θ e un parametro che
caratterizza la distribuzione. Nell’ipotesi che le osservazioni Xi sono indipendenti,
una misura della probabilita di aver ottenuto quelle osservazioni proprio da quella
distribuzione (se θ e il valore del parametro incognito) e data dalla funzione di
verosimiglianza
L(θ) = fθ(X1)fθ(X2) · · · fθ(Xn).
Nel caso di distribuzioni di probabilita discrete con funzione di probabilita data
da pθ(x), la funzione di verosimiglianza e definita da
L(θ) = pθ(X1)pθ(X2) · · · pθ(Xn).
Il metodo della massima verosimiglianza consiste nello scegliere come stimatore
del parametro incognito θ il valore θ che massimizza L(θ).
In generale, puo non essere facile massimizzare la funzione di verosimiglianza e
naturalmente la difficolta aumenta se la distribuzione ha piu di un parametro
e quindi la funzione L sara una funzione di piu variabili. Nell’effettuare que-
sta operazione di massimizzazione invece di considerare la funzione L(θ) si puo
considerare la funzione
`(θ) = lnL(θ)
ovvero il logaritmo della funzione di verosimiglianza (log-likelihood function) che Log-
likelihood
function
potrebbe essere piu facile da massimizzare. Infatti, poiche la funzione logaritmo
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
184 ELEMENTI DI STATISTICA
e una funzione strettamente crescente, un valore θ massimizza L(θ) se e solo se θ
massimizza `(θ). Si osservi, comunque, che l’operazione di massimizzazione non
e sempre facile e che potrebbe anche essere necessario ricorrere a metodi numerici
per risolvere l’equazione che si ottiene dall’annullamento della derivata.
Esempio C.3.3 Date n osservazioni X1, . . . , Xn dalla distribuzione esponenziale, deter-
minare con il metodo della massima verosimiglianza il parametro λ della distribuzione.
La funzione di verosimiglianza e data da
L(λ) = (λe−λX1)(λe−λX1) · · · (λe−λXn) = λne−λ∑ni=1Xi = λne−λnXn .
Uguagliando a zero la derivata (rispetto a λ) si ha
dL(λ)
dλ= nλn−1e−λnXn
(1− λXn
)= 0,
da cui si ha λ =1
Xn. Poiche si ha
d2L(λ)
dλ2< 0, λ e un punto di massimo per la funzione
L(λ). Il valore ottenuto per λ non ci sorprende perche la media campionaria e uno
stimatore corretto della media della distribuzione che e 1/λ.
Allo stesso risultato si puo arrivare considerando la funzione log-likelihood, ovvero
`(λ) = lnL(λ) = n lnλ− λn∑i=1
Xi.
Infatti si had`
dλ=n
λ−
n∑i=1
Xi = 0
per λ =1
Xned inoltre
d2`(λ)
dλ2= − n
λ2< 0.
Esempio C.3.4 Date n osservazioni X1, . . . , Xn dalla distribuzione geometrica, deter-
minare con il metodo della massima verosimiglianza il parametro p ∈ (0, 1) della di-
stribuzione. La distribuzione geometrica e una distribuzione di probabilita discreta
caratterizzata da
pp(x) =
{p(1− p)x se x = 1, 2, . . .
0 altrimenti.
La funzione di verosimiglianza e data da
L(p) = pp(X1)pp(X2) · · · pp(Xn) = p(1−p)X1p(1−p)X2 · · · p(1−p)Xn = pn(1−p)∑ni=1Xi .
Si ha
`(p) = lnL(p) = n ln p+
n∑i=1
Xi ln(1− p).
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
STIMA DI PARAMETRI 185
Annulliamo quindi la derivata della `(p), ottenendo
d`(p)
dp=n
p−
n∑i=1
Xi
1− p= 0.
Risolvendo l’equazione si ottiene
p =1
1 + Xn
e poiche risultad2`(p)
dp2< 0,
allora p e un massimo per la `(p) e quindi anche per la funzione di massima verosimi-
glianza L(p).
Esercizio C.3.5 Siano date n osservazioni indipendenti X1, . . . , Xn estratte dalla distri-
buzione di Poisson di parametro λ. Determinare lo stimatore di massima verosimiglianza
del parametro λ.
C.3.4 Stime per intervalli
Nei metodi di stima puntuale e sempre presente un errore θ − θ dovuto al fatto
che la stima θ in genere non coincide con il parametro θ. Sorge quindi l’esi-
genza di determinare una misura dell’errore commesso. Inoltre, dato il campio-
ne X1, . . . , Xn estratto da una distribuzione di probabilita caratterizzata da un
parametro incognito θ, qualunque sia lo stimatore h(X1, . . . , Xn) scelto per sti-
mare θ, esso dipende dal campione, ovvero lo stimatore fornira stime diverse in
corrispondenza di campioni diversi.
Queste due osservazioni fanno nascere l’esigenza di considerare stime per inter-
valli. Infatti, sulla base dei valori di θ ottenuti considerando un campione casuale
X1, . . . , Xn, si puo definire un intervallo in cui sono compresi i valori piu probabili
per il parametro θ, secondo un “livello di confidenza” fissato. Per fare cio si puo
procedere indirettamente utilizzando una statistica campionaria g(X1, . . . , Xn) la
cui distribuzione sia nota e non dipendente da θ. Naturalmente, visto che la g e
nota, fissato un livello di confidenza (1 − α), e possibile determinare due valori
g1 e g2, indipendenti da θ tali che, comunque scelto α ∈ (0, 1),
P (g1 ≤ g ≤ g2) = 1− α.
Lo scopo e quello di tradurre una probabilita su un intervallo per g in una pro-
babilita su intervallo per θ in modo da poter avere
P (h1 ≤ θ ≤ h2) = 1− α,
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
186 ELEMENTI DI STATISTICA
ovvero in modo tale che h1 e h2 rappresentino gli estremi dell’intervallo per θ.
Le distribuzioni note alle quali si fa di solito riferimento sono la distribuzione
Normale, la distribuzione t di Student e la distribuzione Chi–quadro.
Lo scopo sara quello di avere un valore ben superiore a 0.5 in modo che la pro-
babilita che il parametro θ appartenga all’intervallo [h1, h2] sia tale da assicurare
all’evento h1 ≤ θ ≤ h2 (evento che si verifica nel 100(1−α)% dei casi) una carat-
teristica di “sistematicita”, mentre all’evento complementare (che si verifica nel
100α% dei casi) una caratteristica di “accidentalita”.
Formalmente si puo fornire la seguente definizione.
Definizione C.3.6 Intervallo di confidenza. Dato un campione
X1, . . . , Xn, dato α ∈ (0, 1) e date le statistiche h1 = h1(X1, . . . , Xn) e
h2 = h2(X1, . . . , Xn) con h1 < h2, per le quali
P (h1 ≤ θ ≤ h2) = 1− α,
l’intervallo [h1, h2] si dice intervallo di confidenza per θ con livello di confidenza
pari ad (1− α).
Naturalmente tanto piu e piccolo α, tanto piu e ampio l’intervallo di confidenza
ad esso associato.
C.3.5 Stima per intervalli di una media
Siano date le osservazioni X1, . . . , Xn estratte da una distribuzione di probabilita
a media µ e varianza σ2. Assumiamo inizialmente che la media µ sia incogni-
ta mentre la varianza sia nota. Dalla Proposizione C.2.1 sappiamo che per nCaso
varianza
nota
sufficientemente grande, la variabile aleatoria
Zn =Xn − µσ√n
(C.3.1)
e distribuita approssimativamente secondo la distribuzione Normale standard,
indipendentemente dalla distribuzione delle Xi. Quindi per n sufficientemente
grande risulta
P(−z1−α
2≤ Zn ≤ z1−α
2
)= P
(−z1−α
2≤ Xn − µ
σ√n
≤ z1−α2
)' 1− α,
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
STIMA DI PARAMETRI 187
dove z1−α2
e il punto critico (1 − α2 ) per una distribuzione Normale standard e
puo essere immediatamente determinato dai valori tabulati1
Segue che
P
(Xn − z1−α
2
σ√n≤ µ ≤ Xn + z1−α
2
σ√n
)' 1− α.
Quindi, per n sufficientemente grande, il valore incognito di µ si trova nell’inter-
vallo [Xn − z1−α
2
σ√n
, Xn + z1−α2
σ√n
]con un livello di confidenza pari a 1− α. Si parla di intervallo di confidenza del
100(1− α)%.
Ad esempio, per determinare un intervallo di confidenza al 95% per una media,
dato un campione X1, . . . , Xn, si trova il punto critico z1−α2
= z0.975 = 1.96 dalle
tabelle della distribuzione Normale standard e facilmente si ricava tale intervallo[Xn − 1.96
σ√n
, Xn + 1.96σ√n
].
Si noti che l’intervallo di confidenza puo contenere o meno il valore della media
µ. Questo implica che un singolo intervallo di confidenza, da un punto di vista
probabilistico, fornisce informazioni non complete. La corretta interpretazione
dell’intervallo di confidenza e la seguente: se si costruisce un numero molto elevato
di intervalli di confidenza al 100(1 − α)%, indipendenti e ciascuno basato su n
osservazioni, con n sufficientemente grande, una frazione di questi intervalli pari
ad 1 − α conterra µ. Questa frazione si chiama copertura per l’intervallo di
confidenza. Quindi solamente una proporzione α di casi dara luogo ad intervalli
che non contengono µ.
Osservazione C.3.7 Affermare che l’intervallo di confidenza, ad esempio per
una media µ e al 100(1− α)% non significa che la probabilita che µ appartenga
a questo intervallo e pari a (1 − α)%; questo perche prima di osservare i dati si
puo parlare di probabilita che l’intervallo che si otterra contenga µ in quanto gli
estremi dell’intervallo sono variabili aleatorie. Dopo l’osservazione dei dati, ovvero
dopo aver determinato gli estremi dell’intervallo, si puo solamente affermare che
l’intervallo trovato contiene µ con il 100(1 − α)% di confidenza in quanto non
si ha a che fare con nessuna variabile aletoria (µ e incognita, ma costante e gli
estremi dell’intervallo, una volta determinati sono valori numerici).
Quanto esposto fino ad ora assume che la varianza σ2 sia nota. Se invece σ2 non
si conosce sono necessarie altre considerazioni. Infatti, anche se la varianza σ2 Caso
varianza
non nota1Sono ampiamente disponibili tabelle che, per una variabile Normale standard Z, riportano il valore di
P (Z ≤ z) = 1√2π
∫ z−∞ e−
y2
2 dy.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
188 ELEMENTI DI STATISTICA
non e nota, considerando la varianza campionaria s2n, poiche si ha
limn→∞
s2n = σ2,
il Teorema del Limite Centrale continua a valere se nell’espressione (C.3.1) della
Zn sostituiamo σ2 con la sua stima s2n. Questo significa che, per n sufficientemente
grande la variabile
Zn =Xn − µsn√n
,
e approssimativamente distribuita secondo la distribuzione Normale standard.
Quindi di ha
P
(Xn − z1−α
2
sn√n≤ µ ≤ Xn + z1−α
2
sn√n
)' 1− α,
dove z1−α2
e il punto critico (1− α2 ) per una distribuzione normale standard.
Quindi, per n sufficientemente grande, con livello di confidenza 1 − α il valore
incognito di µ si trova nell’intervallo[Xn − z1−α
2
sn√n
, Xn + z1−α2
sn√n
].
Esempio C.3.8 Siano dati i seguenti valori di dieci osservazioni 1.20, 1.50, 1.68, 1.89,
0.95, 1.49, 1.58, 1.55, 0.50, 1.09, da una distribuzione Normale a media e varianza non
note. Si ricava facilmente che X10 = 1.34 e s210 = 0.17. Volendo costruire un intervallo di
confidenza al 90% per la media si ha 1−α = 0.90, ovvero 1− α2
= 0.95. Dalla tavola dei
valori di P (Z ≤ z) per Z variabile Normale standard si ricava che risulta P (Z ≤ c) = 0.95
per c = 1.65, ovvero z1−α/2 = z0.95 = 1.65. Quindi gli estremi dell’intervallo richiesto
sono dati da
Xn ± z1−α2
√s2n
n= X10 ± z0.95
√s2
10
10= 1.34± 1.65
√0.017.
Pertanto l’intervallo di confidenza richiesto e dato da [1.12 , 1.55].
La difficolta nell’utilizzare questo tipo di intervallo di confidenza per µ sta nel
fatto che esso ha valore asintotico, ovvero per n sufficientemente grande e quindi
risulta approssimato. Inoltre il valore di n per cui la Fn(z) approssima bene
la funzione di distribuzione di una variabile Normale standard dipende dalla
distribuzione delle osservazioni Xi. Se si scelgono valori di n troppo piccoli si
ottiene una copertura di un intervallo di confidenza al 100(1 − α)% inferiore
a 1 − α. In questo caso, ovvero per valori piccoli di n si puo utilizzare una
definizione alternativa dell’intervallo di confidenza che fa riferimento non piu alla
distribuzione Normale standard, ma alla distribuzione t di Student a n− 1 gradi
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
STIMA DI PARAMETRI 189
di liberta. Dalla Proposizione C.2.3 si ha che se le Xi sono variabili Normali, la
variabile
Zn =Xn − µsn√n
ha distribuzione t di Student con n − 1 gradi di liberta per ogni n > 1. Quindi,
per ogni n > 1 un intervallo di confidenza esatto al 100(1− α)% e dato da[Xn − tn−1,1−α
2
sn√n
, Xn + tn−1,1−α2
sn√n
],
dove tn−1,1−α2
e il punto critico 1− α2
per una distribuzione t di Student ad n− 1 t intervallo
di confiden-
zagradi di liberta2. Questo intervallo si chiama t intervallo di confidenza.
Poiche risulta tn−1,1−α2> z1−α
2l’intervallo di confidenza definito in riferimento
alla distribuzione t di Student e piu ampio di quello definito in riferimento alla
distribuzione Normale standard.
Esempio C.3.9 Volendo determinare il t intervallo di confidenza nel caso delle osserva-
zioni dell’Esempio C.3.8 e sufficiente ricavare dai valori tabulati il valore di t9,0.95 che e
pari a 1.83. Quindi l’intervallo di confidenza richiesto e dato da [1.10 , 1.58].
Tuttavia anche questo secondo tipo di intervallo di confidenza presenta aspetti
problematici; infatti, nella sua definizione si assume che le osservazioni Xi sono
estratte da una distribuzione Normale e poiche questa assuzione, in generale, non
e verificata, anche questo tipo di intervallo di confidenza e, di fatto, approssimato.
In conclusione, possiamo dire che l’intervallo di confidenza definito in riferimento
alla distribuzione Normale standard e basato sul Teorema del Limite Centrale
e la copertura dipende dalla scelta di n. L’intervallo di confidenza definito in
riferimento alla distribuzione t di Student e approssimato perche influenzato dalla
distribuzione delle Xi che in generale non sono Normali; tuttavia questo secondo
tipo di intervallo di confidenza ha maggiore copertura dell’altro.
2Anche in questo caso sono disponibili tavole che riportano valori tabulati di tali punti critici
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
190 ELEMENTI DI STATISTICA
C.4 TEST DELLE IPOTESI
Nel cercare di costruire un legame tra dati osservati e ipotesi teoriche sulle carat-
teristiche dell’intera popolazione si deve, in genere, prendere una decisione per il
raggiungimento di tale conclusione generale e nasce il problema di esprimere un
giudizio di plausibilita di un’ipotesi che si e specificata per la popolazione. Per
verificare la coerenza tra osservazioni e ipotesi fatta si fa uso di test statistici che
prendono nome di test delle ipotesi. In sintesi, possiamo dire che tali test devono
confrontare i valori osservati e i corrispondenti valori teorici attesi condiziona-
tamente all’ipotesi fatta. Le differenze che vengono riscontrate possono essere
ovviamente ricondotte a due possibilita:
• l’ipotesi specificata e corretta e la differenza riscontrata e puramente casua-
le;
• l’ipotesi specificata e errata e quindi non ci si puo aspettare che i due valori
siano “vicini”.
Il procedimento consiste nel confrontare due ipotesi: l’ipotesi da sottoporre a ve-
rifica e il suo complemento. Si indica con H0 il sottoinsieme dei valori individuati
dall’ipotesi da sottoporre a verifica che viene detta ipotesi nulla, mentre il suo
complemento si indica con H1 e viene detto ipotesi alternativa. Si osservi che se
un test di ipotesi non scarta l’ipotesi H0, questo non vuol dire che H0 e accettata
come vera, ma solamente che essa non deve essere scartata, ovvero che puo essere
considerata possibile.
Si parla di errore di I specie se il test porta a rifiutare un’ipotesi H0 quando questa
e corretta e di errore di II specie se il test porta ad accettare H0 quando questa e
falsa. Si noti che l’obiettivo non e quello di dire se l’ipotesi fatta e vera o falsa, ma
piuttosto di verificare se l’ipotesi fatta sia compatibile con i dati. In genere, c’e
un ampio margine di tolleranza nell’accettare H0, mentre per rifiutarla occorre
che i dati siano veramente poco probabili quando H0 dovesse essere corretta. Per
ottenere questo, si specifica un valore α, detto livello di significativita e si impone
che il test sia tale che, quando l’ipotesi H0 e corretta, la probabilita che essa venga
scartata e non superiore ad α. Quindi un test con livello di significativita pari
ad α deve essere tale che una probabilita di commettere un errore di I specie e
minore o uguale ad α.
Rimandiamo alla letteratura specifica per una trattazione generale dei test d’ipo-
tesi; nel seguito faremo esclusivamente riferimento al problema di decidere se le
osservazioni sono un campione indipendente di una particolare distribuzione di
probabilita con funzione di distribuzione F . Ovvero, si vuole utilizzare un test
delle ipotesi per avvalorare o smentire un’ipotesi fatta sulla distribuzione di pro-
babilita che meglio rappresenta tali dati. Quindi, date le osservazioni X1, . . . , Xn,
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
TEST DELLE IPOTESI 191
l’ipotesi da sottoporre a verifica e la seguente:
H0 ={X1, . . . , Xn sono variabili aleatorie indipendenti,
identicamente distribuite con funzione di distribuzione F}
Esistono diversi test di ipotesi; consideriamone brevemente nel seguito due, ri-
mandando alla letteratura specifica per una trattazione completa.
C.4.1 Test Chi–quadro
Siano X1, . . . , Xn variabili aleatorie indipendenti identicamente distribuite. Per
applicare il test si effettua il seguente procedimento: si suddivide il range della
distribuzione in k intervalli adiacenti
[a0, a1) , [a1, a2) , [a2, a3) , . . . [ak−1, ak) ,
dove puo anche essere a0 = −∞ e/o ak =∞. Si definisce ora Nj come il numero
delle Xi contenute nell’intervallo j-esimo, ovvero in [aj−1, aj) per j = 1, 2, . . . k.
Naturalmente risulta∑k
j=1Nj = n. Si calcola la proporzione pj che le Xi cadono
nell’intervallo j-esimo se le osservazioni fossero dalla distribuzione ipotizzata,
ovvero
pj =
∫ aj
aj−1
f(x)dx nel caso continuo∑aj−1≤xi<aj
p(xi) nel caso discreto
dove f(x) e la densita di probabilita della distribuzione ipotizzata nel caso con-
tinuo e p(xi) sono i valori della distribuzione di probabilita nel caso discreto. Si
definisce
χ2 =
k∑j=1
(Nj − npj)2
npj. (C.4.1)
Poiche npj e il numero atteso che n Xi cadano nell’j-esimo intervallo se l’ipotesi
H0 fosse vera, ci si puo aspettare che χ2 sia piccolo se la distribuzione ipotizzata
realizza un buon fitting delle X1, . . . , Xn e quindi l’ipotesi H0 e scartata se χ2 e
troppo grande.
Per effettuare praticamente il test dobbiamo considerare il caso in cui tutti i
parametri della distribuzione sono noti e il caso in cui questo non avviene. Il test
Chi-quadro puo essere effettuato in entrambi i casi.
Nel primo caso, se α e il livello di significativita del test, si vuole trovare la Caso para-
metri noticosiddetta “regione critica”, ovvero si vuole calcolare un valore z tale che
PH0(χ2 ≥ z) = α,
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
192 ELEMENTI DI STATISTICA
dove la notazione PH0 indica il condizionamento della probabilita al fatto che H0
sia vera. Il test dovra rifiutare l’ipotesi nulla quando il valore osservato per χ2 e
superiore a z. Il valore critico z puo essere ottenuto utilizzando il fatto che, se
l’ipotesi H0 e vera, allora χ2 per n → ∞ converge alla distribuzione Chi–quadro
con k− 1 gradi di liberta. Percio, per n grande, un test con livello di accuratezza
α e ottenuto scartando H0 se
χ2 > χ2k−1,1−α
dove χ2k−1,1−α e il punto critico 1− α per la distribuzione Chi–quadro con k − 1
gradi di liberta3.
Se invece, alcuni parametri fossero non noti, c’e da considerare il fatto che le pjCaso con
parametri
non noti
derivano dalla distribuzione ipotizzata, ma con parametri stimati. In particolare,
supponendo che siano presenti m parametri non noti e che siano stimati con il
metodo della massima verosimiglianza, si usano queste stime per calcolare
χ2 =
k∑j=1
(Nj − npj)2
npj,
dove rispetto alla (C.4.1), le pj sono state sostituite dalle pj che sono ottenute
utilizzando la distribuzione con i parametri stimati. In questo caso si puo dimo-
strare che, se l’ipotesi H0 e vera, allora χ2 per n→∞ converge alla distribuzione
Chi–quadro con k−1−m gradi di liberta (si perdono tanti gradi di liberta quanti
sono i parametri). Il test con livello di significativita α diventa il seguente: si
scarta H0 se
χ2 > χ2k−1−m,1−α,
dove χ2k−1−m,1−α e il punto critico 1 − α per la distribuzione Chi–quadro con
k − 1−m gradi di liberta
L’aspetto piu problematico nell’effettuare il test Chi–quadro e la scelta del nume-
ro e dell’ampiezza degli intervalli. E un problema difficile e non c’e una risposta
definitiva, ma solo delle linee guida da seguire per le quali si rimanda ai testi spe-
cifici. Diciamo solo che e opportuno scegliere gli intervalli in modo tale che risulti
p1 = p2 = · · · = pk (equiprobabilita) e, in questo caso, basandosi su osservazioni
empiriche si e visto che il test Chi-quadro e valido se k ≥ 3 e npj ≥ 5 nel caso con-
tinuo; nel caso discreto e opportuno avere tutti i valori npj approssimativamente
uguali e pari ad almeno a 5.
Nonostante queste difficolta il test Chi–quadro e molto utilizzato perche puo
essere applicato ad ogni distribuzione ipotizzata, mentre altri test non hanno
questa flessibilita.
3Sono disponibili tavole che riportano i valori di questi punti critici
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
TEST DELLE IPOTESI 193
Esercizio C.4.1 Supponiamo di monitorare il numero di automobili che passano ad
un casello autostradale tra le 7.00 e le 7.05 di giorni feriali per 100 giorni, ottenendo il
seguente risultato (numero di auto/frequenza)
numero automobili: 0 1 2 3 4 5 6 7 8 9 10 11
frequenze: 12 10 19 17 10 8 7 5 5 3 3 1
Ovvero, in 12 casi e stato osservato il passaggio di 0 auto, in 10 casi e stato osservato
il passaggio di 1 auto, etc. Supponiamo di essere nella fase di analisi dell’input di una
simulazione e di voler quindi scegliere una distribuzione di input che meglio rappresenti
questi dati. A tale scopo, supponiamo di aver scelto come distribuzione teorica ipotizzata
la distribuzione di Poisson di parametro 3.64. Applicare il test di ipotesi Chi-quadro con
un livello di significativita α = 0.05 per valutare l’adattamento di questa distribuzione
ai dati, ovvero si assuma p(x) =e−λλx
x!, x = 0, 1, 2, . . ., con λ = 3.64. Nell’applicazione
del test considerare i seguenti 7 intervalli: [0, 2), [2, 3), [3, 4), [4, 5), [5, 6), [6, 7),
[7, 11]. Considerare i due casi in cui il parametro λ e stimato oppure e noto.
C.4.2 Test di Kolmogorov–Smirnov
Un altro test di ipotesi che viene spesso preso in considerazione e basato sul con-
fronto tra una funzione di distribuzione empirica e la funzione di distribuzione F
della distribuzione ipotizzata: si tratta del test di Kolmogorov–Smirnov che come
pregio maggiore ha quello di essere valido per ogni dimensione n del numero delle
osservazioni, mentre il test chi–quadro e valido in senso asintotico. L’inconve-
niente maggiore e invece rappresentato dalla limitata applicabilita del test che,
almeno nella sua forma originaria, e valido se tutti i parametri della distribuzione
ipotizzata sono noti, ovvero non stimati dai dati e se la distribuzione e continua.
Esistono estensioni del test al caso in cui i parametri sono stimati solo in riferi-
mento a particolari distribuzioni come la distribuzione normale e la distribuzione
esponenziale.
Siano X1, . . . , Xn variabili aleatorie indipendenti con funzione di distribuzione
F e siano dati i valori x1, . . . , xn i valori osservati di X1, . . . , Xn, ovvero siano
essi una realizzazione delle variabili aleatorie. Si puo definire un’approssimazione
empirica Fe della F nel seguente modo:
Fe(x) =|{i : xi ≤ x}|
n=
numero delle xi ≤ xn
.
Fe e la percentuale dei dati del campione che sono minori o uguali di x, ovvero
Fe e la funzione di ripartizione della variabile aleatoria discreta che assume con
la stessa probabilita gli n valori osservati.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
194 ELEMENTI DI STATISTICA
In pratica, ordinando le xi per valori crescenti in modo che si abbia
x(1) ≤ x(2) ≤ · · · ≤ x(n), risulta
Fe(x) =
0 se x < x(1)
i
nse x(i) ≤ x ≤ x(i+1)
1 se x(n) ≤ x.
Se F (x) e la funzione di distribuzione ipotizzata, una valutazione naturale della
“bonta” dell’approssimazione e data da una qualche misura di “vicinanza” tra le
funzioni Fe e F . Una misura dello scostamento puo essere definita da
D = max{∣∣∣Fe(x)− F (x)
∣∣∣} .D e detta statistica del test di Kolmogorov–Smirnov.
Considerando separatamente il caso in cui Fe(x) − F (x) ≥ 0 e il caso in cui
Fe(x) − F (x) ≤ 0, poiche F e una funzione non decrescente, si possono definire
le quantita D+ e D− e risulta
D+ = maxx
{Fe(x)− F (x)
}= max
1≤i≤n
{i
n− F (x(i))
}D− = max
x
{F (x)− Fe(x)
}= max
1≤i≤n
{F (x(i))−
i− 1
n
}.
Quindi si ha
D = max{D+ , D−
}.
Chiaramente, un valore grande di Dn indica che l’approssimazione non e molto
buona e quindi si scarta l’ipotesi H0 se
D > dn,1−α, (C.4.2)
dove dn,1−α sono costanti che, in generale, possono dipendere dalla distribuzione
ipotizzata e dove α e il livello di accuratezza del test.
E importante notare che se tutti i parametri della F sono noti e non stimati,
allora si puo dimostrare che la distribuzione di D non dipende dalla particolare
distribuzione ipotizzata e questo permette di utilizzare una sola tabella di valori
per dn,1−α per tutte le distribuzioni continue.
Sempre nel caso di tutti i parametri noti, e stato verificato che una buona
approssimazione si ottiene utilizzando al posto di (C.4.2) il seguente test(√n+ 0.12 +
0.11√n
)D > c1−α
dove c1−α sono costanti non dipendenti da n delle quali si dispone di valori
tabellati riportati nella Tabella C.4.1.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
TEST DELLE IPOTESI 195
1− α 0.85 0.90 0.95 0.975 0.99
c1−α 1.138 1.224 1.358 1.480 1.628
Tabella C.4.1 Valori delle costanti c1−α nel caso in cui tutti i parametri sono noti
Esistono alcuni casi particolari in cui il test di Kolmogorov–Smirnov puo essere
applicato anche se i parametri sono non noti, ovvero sono stimati. Due casi impor-
tanti in cui questo accade sono il caso in cui la distribuzione ipotizzata e Normale
con media µ e varianza σ2 non note, ma entrambe stimate rispettivamente con
Xn e s2n e il caso in cui la distribuzione ipotizzata e esponenziale di parametro
incognito λ essendo λ stimato con il metodo della massima verosimilianza ovvero
λ = 1/Xn.
1) Nel caso in cui la distribuzione ipotizzata e Normale con media e varian-
za stimate rispettivamente con Xn e s2n, si definisce la funzione F come
la funzione di distribuzione di una Normale a media Xn e varianza s2n.
Utilizzando questa funzione F si procede analogamente al caso precedente
calcolando D; cambia pero il test in quanto e stato visto che un buon livello
di accuratezza si ottiene scartando H0 se risulta(√n− 0.01 +
0.85√n
)D > c′1−α
dove c′1−α sono i valori riportati nella Tabella C.4.2.
1− α 0.85 0.90 0.95 0.975 0.99
c′1−α 0.775 0.819 0.895 0.955 1.035
Tabella C.4.2 Valori delle costanti c′1−α nel caso in cui la distribuzione ipotizzata e Normale
2) Nel caso in cui la distribuzione ipotizzata e esponenziale in cui il parametro
incognito λ e stimato con il metodo della massima verosimilianza ovvero
λ = 1/Xn, si definisce la funzione F come la funzione di distribuzione della
distribuzione esponenziale di parametro 1/Xn, ovvero F (x) = 1 − e−x/Xn ,
x ≥ 0, e si calcola D. In questo caso il test prevede di scartare H0 se(D − 0.2
n
)(√n+ 0.26 +
0.5√n
)> c′′1−α
dove c′′1−α sono i valori riportati nella Tabella C.4.3.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
196 ELEMENTI DI STATISTICA
1− α 0.85 0.90 0.95 0.975 0.99
c′′1−α 0.926 0.990 1.094 1.190 1.308
Tabella C.4.3 Valori delle costanti c′′1−α nel caso in cui la distribuzione ipotizzata e esponenziale
Esercizio C.4.2 Nel definire una distribuzione di probabilita durante l’analisi dell’input
di una simulazione, si hanno a disposizione i seguenti dati:
81, 72, 94, 66.
Utilizzare il test di Kolmogorov-Smirnov per verificare l’ipotesi che la distribuzione espo-
nenziale con media 100 approssima bene questi dati con un livello di significativita α pari
al 2.5%.
C.4.3 Il p-value
Il risultato di un test delle ipotesi porta a scartare o meno l’ipotesi H0: essa si
scarta quando la statistica del test (ovvero χ2 nel test Chi–quadro o D nel test
Kolmogorv–Smirnov) superano il valore del punto critico. Tuttavia non viene
specificato quanto la statistica del test dista dal valore critico. Si puo allora
definire una probabilita che fornisce una misura della distanza della statistica del
test dal valore critico: si definisce p-value o p-dei-dati come l’estremo inferiore
dei livelli di significativita che porterebbe a rifiutare H0. Quindi se il livello α e
maggiore di tale probabilita (il p-value), l’ipotesi H0 viene scartata, altrimenti
non viene scartata. La probabilita data dal p-value fornisce quindi il livello di
significativita critico scendendo al di sotto del quale la decisione cambia da rifiuto
dell’ipotesi a non rifiuto.
Poiche il p-value e il livello di significativita al di sopra del quale si scarta H0,
valori grandi del p-value indicano un buon adattamento della distribuzione ipotiz-
zata ai dati, mentre valori piccoli del p-value indicano che il fitting non e buono.
Sulla base di cio, molti software operano calcolando il p-value in corrisponden-
za di diverse distribuzioni e decidono che la piu adatta a rappresentare i dati e
quella alla quale corrisponde il p-value piu grande. Questa opzione disponibile in
molti software deve essere tuttavia eventualemte utilizzata con estrema cautela,
in quanto e possibile che sulla sola base del valore del p-value non si abbiano
risposte adeguate.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
4
4.1 ESERCIZI DI RIEPILOGO
Esercizio 4.1.1 Dato un sistema M/M/2 con frequenza media di arrivo pari a λ
e velocita di servizio pari a µ, con λ < 2µ
• Scrivere le equazioni di Kolmogorov relative al processo di nascita e morte
al quale si puo ricondurre questo sistema;
• ricavare un’espressione di pn in funzione di λ e µ.
Esercizio 4.1.2 In un fast–food e previsto il servizio “drive”, ovvero si raggiunge
con l’auto una piattaforma e si viene serviti senza scendere dall’auto. Viene
servita un’auto alla volta e attualmente c’e spazio solamente per 3 auto in attesa.
Gli arrivi dei clienti sono poissoniani con in media 40 auto l’ora e i tempi di
servizio sono distribuiti esponenzialmente con media di 50 auto l’ora. Il profitto
medio che il gestore ricava dall’effettuazione di un servizio e di 0.50 Euro. Inoltre
il gestore puo affittare altri spazi auto a 9 Euro al giorno. Il fast–food e aperto
24 ore su 24. Determinare se e conveniente affittare posti auto aggiuntivi e, in
caso affermativo, calcolare quanti.
Esercizio 4.1.3 In un ciclo produttivo di un’industria c’e una stazione di ispezione
il cui servizio consiste nell’effettuare un controllo su un pezzo semilavorato e che
ha un ciclo operativo costante di 6 minuti. Gli arrivi alla stazione di ispezione
sono poissoniani alla media di 4 pezzi l’ora. Determinare
• il numero medio di pezzi presenti in coda in attesa di essere lavorati dalla
stazione di ispezione;
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
198
• il numero medio di pezzi presenti nel sistema;
• il tempo medio di attesa di un pezzo in coda prima di essere lavorato;
• il tempo medio di permanenza nel sistema.
Senza effettuare il calcolo diretto, ma utilizzando i valori gia calcolati, dedurre i
valori del tempo medio di attesa in coda e il numero medio di pezzi in coda nel
caso in cui i tempi di servizio fossero esponenziali.
Esercizio 4.1.4 Presso un ambulatorio medico arrivano i pazienti ad intervalli
medi di 20 minuti e sono necessari, in media 15 minuti per la visita medica. Si
assume che gli arrivi sono poissoniani e che i tempi di durata delle visite sono
esponenziali. Il medico desidera avere sufficienti posti a sedere tra la sala d’attesa
dell’ambulatorio e quello su cui siede il paziente durante la visita in modo che non
piu dell’1% dei pazienti deve attendere in piedi. Quante sedie dovra acquistare ?
Esercizio 4.1.5 Supponiamo che in 10 repliche indipendenti di una simulazione,
una variabile aleatoria assuma i seguenti valori Xj : 9.98, 10.07, 9.94, 10.22, 9.98,
10.01, 10.11, 10.01, 9.99, 9.92. Determinare la stima della media e l’intervallo t
di confidenza al 95%.
Esercizio 4.1.6 In un ufficio postale ci sono 6 sportelli cassa che operano dalle 8
alle 14. Assumiamo che i clienti arrivino secondo Poisson con frequenza di 2 al
minuto, e che i tempi di servizio di ciascun servente sono uniformemente distri-
buiti tra 2 e 3 minuti. Effettuando 10 repliche indipendenti di una simulazione,
si ottengono i seguenti valori per il numero medio degli utenti in coda
numero replica valore
1 1.49
2 1.64
3 1.95
4 1.45
5 1.51
6 1.55
7 1.88
8 1.36
9 1.78
10 1.46
Determinare se le 10 repliche effettuate sono sufficienti per stimare il numero
medio di utenti in coda con un errore relativo pari a 0.07 e un intervallo di
confidenza al 90%. Dire, inoltre, motivando la risposta, se si possono utilizzare
le formule per le misure di prestazione di un sistema di code M/G/s per studiare
questo problema.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
BIBLIOGRAFIA 199
Bibliografia
Altiok, T., Melamed, B. (2007). Simulation Modeling and Analysis with
ARENA. Academic Press, Amsterdam.
Banks, J. (1998). Handbook on Simulation. Principles, Methodology, Advances,
Applications, and Practice. Wiley and Sons, New York.
Banks, J., J.S.Carson, Nelson, B., Nicol, D. (2001). Discrete–event system
simulation. Prentice Hall.
Bertsekas, D., Gallager, R. (1991). Data Networks. Prentice Hall, Upper
Saddle River, NJ, second edition.
Bhat, U. (2008). An Introduction to Queueing Theory. Birkhauser, Boston.
Billingsley, P. (1979). Probability and mesaure. Wiley and Sons, New York.
Bratley, P., Fox, B., Schrage, L. (1987). A guide to simulation. Springer–
Verlag, New York, second edition.
Cooper, R. (1981). Introduction to queueing theory. North Holland, New York.
Fu, M., editor (2014). Handbook on Simulation Optimziation. Springer, New
York.
Gross, D., Harris, C. (1998). Fundamental of queueing theory. Wiley and
Sons, New York, third edition.
Hillier, F., Lieberman, G. (2001). Introduction to Operations Research.
McGraw – Hill, New York.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
200 BIBLIOGRAFIA
Joines, J., Roberts, S. (2012). Simulation Modeling with SIMIO: a workbook.
SIMIO LCC, Sewickley, second edition.
Kelton, W., Sadowski, R., Sturrock, D. (2004). Simulation with Arena.
McGraw – Hill, New York, third edition.
Kelton, W., Smith, J., Sturrock, D. (2011). Simio & Simulation. Modeling,
Analysis, Applications. McGraw – Hill, Boston, second edition.
Kleinrock, L. (1975). Queueing systems, vol. I: Theory. Wiley, New York.
Disponibile nella traduzione italiana: Sistemi a coda, Introduzione alla teoria
delle code, Hoepli, 1992.
Kleinrock, L. (1976). Queueing systems, vol. II: Computer Applications. Wiley,
New York.
Law, A., Kelton, W. (2000). Simulation modeling and analysis. McGraw –
Hill, New York, third edition.
Ross, S. (2002). Simulation. Academic Press, San Diego, third edition.
Ross, S. (2003a). Introduction to probability models. Academic Press, San Diego.
Ross, S. (2003b). Probabilita e statistica. Apogeo, Milano.
Rossetti, M. (2010). Simulation Modeling and Arena. Wiley, Hoboken.
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
Indice
Prefazione iii
Introduzione v
1 Teoria delle code 1
1.1 Generalita 3
1.1.1 Esempi reali di sistemi a coda 4
1.1.2 Componenti di un sistema di servizio 5
1.1.3 Notazione di Kendall 8
1.2 Problematiche di interesse e relazioni fondamentali 11
1.2.1 Definizioni e notazioni standard 11
1.2.2 Misure di prestazione 13
1.2.3 Relazioni fondamentali 17
1.3 Modelli stocastici dei processi di arrivo e di servizio 22
1.3.1 La distribuzione esponenziale e le sue proprieta 22
1.3.2 Il processo di Poisson 27
1.3.3 Altri modelli per processi di arrivo 32
1.4 Processi di nascita e morte 35
1.4.1 Preliminari 35
1.4.2 Proprieta di Markov. Catene di Markov a tempocontinuo 35
1.4.3 Processi di nascita e morte 37
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
202 INDICE
1.4.4 Le equazioni di Kolmogorov e la soluzionestazionaria 40
1.4.5 Equazioni di bilancio 45
1.4.6 Processi di sole nascite 47
1.4.7 Sistemi a coda e processi di nascita e morte 49
1.5 Sistemi a coda basati su processi di nascita e morte 50
1.5.1 Sistemi con arrivi poissoniani: proprieta PASTA 51
1.5.2 Sistemi M/M/s 54
1.5.3 Sistemi M/M/s/K 67
1.5.4 Sistemi M/M/s/ · /U 73
1.5.5 Sistemi con velocita di servizio e frequenza di arrivodipendenti dallo stato 77
1.5.6 Sistemi M/M/∞ 80
1.6 Modelli di code con disciplina della coda basata su criteridi priorita 82
1.7 Sistemi a coda con distribuzioni non esponenziali 84
1.7.1 Sistemi M/G/1 84
1.7.2 Sistemi con arrivi non poissoniani 90
1.8 Reti di code 91
1.8.1 Generalita 91
1.8.2 Il processo delle partenze per i sistemi M/M/s 93
1.8.3 Serie di code 94
1.8.4 Determinazione della frequenza media effettivadegli arrivi in ciascun nodo di reti aperte 96
1.8.5 Reti di Jackson aperte 97
1.8.6 Reti di Jackson chiuse 104
1.9 Riferimenti del Capitolo 1 108
1.10 Esercizi di riepilogo 110
1.10.1 Esercizi sui sistemi a coda 110
1.10.2 Esercizi sulle reti di code 114
Appendice A: Formulario – Sistemi a coda 117
A.1 Modelli basati su processi di nascita e morte 117
A.1.1 Sistemi M/M/1 117
A.1.2 Sistemi M/M/s 118
A.1.3 Sistemi M/M/1/K 118
A.1.4 Sistemi M/M/s/K 119
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
INDICE 203
A.1.5 Sistemi M/M/1/ · /U 119
A.1.6 Sistemi M/M/s/ · /U 120
A.1.7 Sistemi con λ e µ dipendenti dallo stato – singoloservente 120
A.1.8 Sistemi con λ e µ dipendenti dallo stato –multiservente 121
A.1.9 Sistemi M/M/∞ 121
A.2 Modelli con distribuzioni non esponenziali 122
A.2.1 Sistemi M/G/1 122
A.2.2 Sistemi M/D/s 122
A.2.3 Sistemi M/Ek/1 122
2 Simulazione 123
2.1 Generalita sui modelli di simulazione 124
2.1.1 Elementi di un modello di simulazione 124
2.1.2 Classificazione dei modelli si simulazione 126
2.1.3 Simulazione ad eventi discreti 126
2.1.4 Schema dello studio di un problema basato sullasimulazione 128
2.1.5 Applicazioni tipiche della simulazione 132
2.2 Analisi dell’input 133
2.2.1 Generalita 133
2.2.2 Scelta delle distribuzioni di input 134
2.2.3 Distribuzioni empiriche 135
2.2.4 Distribuzioni teoriche 136
2.2.5 Scelta di una distribuzione teorica 137
2.2.6 Scelta delle distribuzioni di input in assenza di dati 142
2.3 Generazione di osservazioni casuali 143
2.3.1 Generazione di numeri pseudocasuali condistribuzione uniforme 143
2.3.2 Generazione di osservazioni casuali da unadistribuzione di probabilita 145
2.4 Progettazione di una simulazione e analisi dell’output 151
2.4.1 Dati di output 151
2.4.2 Transitorio e stato stazionario 152
2.4.3 Tipi di simulazioni 153
2.4.4 Simulazioni con terminazione: analisi del transitorio 154
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019
204 INDICE
2.4.5 Simulazioni senza terminazione: analisi dello statostazionario 158
2.5 Tecniche per la riduzione della varianza 162
2.5.1 Variabili antitetiche 162
2.5.2 Condizionamento 164
2.6 L’approccio simulazione ottimizzazione (“simulation–based optimization”) 166
2.7 Software di simulazione (SIMULATORI) 169
2.7.1 ARENA 170
2.7.2 Modelli in ARENA 170
2.8 Riferimenti del Capitolo 2 174
2.9 Esercizi sulla simulazione 175
Appendice C: Elementi di statistica 177
C.1 Elementi di Statistica inferenziale 177
C.2 Statistiche campionarie e loro distribuzione 178
C.2.1 Media campionaria e varianza campionaria 178
C.2.2 Distribuzione (approssimata) della mediacampionaria 179
C.2.3 Distribuzioni delle statistiche di popolazioni normali 180
C.3 Stima di parametri 181
C.3.1 Proprieta degli stimatori 182
C.3.2 Stima di media e varianza 182
C.3.3 Metodi di stima 183
C.3.4 Stime per intervalli 185
C.3.5 Stima per intervalli di una media 186
C.4 Test delle ipotesi 190
C.4.1 Test Chi–quadro 191
C.4.2 Test di Kolmogorov–Smirnov 193
C.4.3 Il p-value 196
4 197
4.1 Esercizi di riepilogo 197
M. Roma – Sistemi di Servizio e Simulazione – SAPIENZA Universita di Roma – a.a. 2018-2019