Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che...

82
2 Simulazione 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 realt` a 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 probabilit` a per generare casualmente eventi del sistema e dal sistema simulato si ottengono osservazioni statistiche sulle prestazioni dello stesso. Naturalmente affinch´ e ci` o 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 Universit`a di Roma a.a. 2018-2019

Transcript of Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che...

Page 1: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 2: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 3: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 4: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 5: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 6: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 7: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 8: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 9: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 10: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 11: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 12: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 13: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 14: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 15: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 16: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 17: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 18: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 19: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 20: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 21: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 22: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 23: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 24: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 25: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 26: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 27: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 28: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 29: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 30: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 31: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 32: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 33: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 34: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 35: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 36: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 37: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 38: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 39: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 40: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 41: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 42: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 43: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 44: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 45: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 46: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 47: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 48: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 49: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 50: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 51: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 52: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 53: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 54: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.
Page 55: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 56: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 57: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 58: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 59: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 60: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 61: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 62: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 63: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 64: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 65: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 66: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 67: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 68: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 69: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 70: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 71: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 72: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 73: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 74: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 75: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 76: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 77: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 78: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 79: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 80: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 81: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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

Page 82: Simulazione - uniroma1.itroma/didattica/SSS18-19/parteE.pdfo calcolata in base ad altri eventi che accadono nel sistema. Un esempio e dato dal tempo di servizio in un sistema a coda.

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