PREZENTAREA CURSULUI Cursul este structurat pe 4 capitole care acopera in intregime programa prezentata. CAPITOLUL I Generalităţi despre modelare si simulare I.1. Introducere I.2. Construcţia unui model de simulare I.3. Concepte de bază în modelarea sistemelor I.3.1. Realizari iterative ale sistemelor I.3.2. Modelarea sistemelori cu evenimente discrete I.3.3. Metodologia de realizare a experimentelor de simulare I.4. Generalităţi despre limbajul GPSS I.5. Exemple de programe GPSS. I.5.1. Model de simulare pentru un sistem de aşteptare cu o staţie I.5.2. Model de simulare pentru un sistem de aşteptare cu staţii paralele I.5.3. Model cu staţii paralele şi preferinţe CAPITOLUL II Numere aleatoare II.1. Noţiuni introductive II.2. Necesitatea simulării numerelor aleatoare II.3. Metode congruenţiale liniare CAPITOLUL III Simularea variabilelor aleatoare neuniforme III.1. Metoda inversă III.2. Metoda compunerii sau amestecării III.3. Metoda respingerii III.4. Alte metode III.5. Simularea repartiţiilor înrudite cu repartiţia normală III.6. Simularea unor variabile particulare III.6.1. Repartiţia exponenţială III.6.2. Repartiţia Gamma III.6.3. Metode pentru simularea variabilei ( ) υ , 1 , 0 Gamma , 1 > υ III.6.4. Repartiţia Beta III.6.5. Repartiţia normală III.7. Simularea unor variabile discrete III.8. Validarea generatorilor
Transcript
PREZENTAREA CURSULUI
Cursul este structurat pe 4 capitole care acopera in intregime programa prezentata.
CAPITOLUL I Generalităţi despre modelare si simulare
I.1. Introducere I.2. Construcţia unui model de simulare I.3. Concepte de bază în modelarea sistemelor I.3.1. Realizari iterative ale sistemelor
I.3.2. Modelarea sistemelori cu evenimente discrete I.3.3. Metodologia de realizare a experimentelor de simulare
I.4. Generalităţi despre limbajul GPSS I.5. Exemple de programe GPSS.
I.5.1. Model de simulare pentru un sistem de aşteptare cu o staţie I.5.2. Model de simulare pentru un sistem de aşteptare cu staţii paralele I.5.3. Model cu staţii paralele şi preferinţe
Simularea variabilelor aleatoare neuniforme III.1. Metoda inversă III.2. Metoda compunerii sau amestecării III.3. Metoda respingerii III.4. Alte metode III.5. Simularea repartiţiilor înrudite cu repartiţia normală III.6. Simularea unor variabile particulare
III.6.1. Repartiţia exponenţială III.6.2. Repartiţia Gamma III.6.3. Metode pentru simularea variabilei ( )υ,1,0Gamma , 1>υ
CAPITOLUL I GENERALITĂŢI DESPRE MODELARE SI SIMULARE
1.1. Introducere
Cuvântul simulare este de origine latină (simulatio) şi înseamnă capacitatea de a
reproduce ceva. În informatică, termenul de simulare a fost introdus de John von Neumann la începutul anilor `40, o dată cu apariţia primelor calculatoare electronice. J.von Neumann împreună cu grupul de savanţi de la Şcoala Los Alamos (Ulam, Metropolis etc) au dezvoltat primele aplicaţii ale calculatoarelor. Tot ei au introdus denumirile Cercetări operaţionale (pentru a desemna aplicaţiile legate de dirijarea operaţiilor militare pe arii geografice mari ale globului pământesc!) precum şi metoda Monte - Carlo (pentru a desemna aplicaţii ale calculatoarelor bazate pe utilizarea numerelor aleatoare). În accepţiunea actuală a informaticii, simularea cuprinde o serie de aplicaţii care realizează imitarea comportamentului unor părţi ale lumii reale (simularea stochastică), luând în considerare şi comportamentul aleator al acesteia. Din domeniul simulării face parte şi metoda Monte Carlo. Definiţia 1. Simularea este o tehnică de realizare a experimentelor cu calculatorul, care implică utilizarea unor modele matematice şi logice care descriu comportarea unui sistem real (sau a unor componente ale sale) de-a lungul unor perioade mari de timp.
Deci simularea se realizează pe baza unui model special, numit model de simulare, cu ajutorul căruia se realizează experimentele prin intermediul calculatorului. Modelul de simulare se construieşte pe scheletul unui model matematic şi se finalizează într-un algoritm. De aceea în cele ce urmează vom prezenta schema generală a unui model de simulare, pornind de la descrierea principalelor elemente ale unui model matematic.
Model matematic. Prin definiţie, un model este un analog ce reprezintă o parte a lumii înconjurătoare într-un mod uşor de perceput de către noi. Modelul ne reprezintă uneori realitatea prin scheme, figuri geometrice sau alte obiecte ce ne sunt familiare şi pe care le înţelegem uşor (i.e. la fel de bine cum le ''vedem'' sau le ''pipăim''). Modelul matematic ne reprezintă realitatea folosind elemente sau abstracţiuni matematice.
Elementele constitutive ale unui model matematic sunt următoarele: a) Variabile (V ) şi Parametri ( P ) care pot fi de intrare ( PIVI , ), dacă pot fi percepute
(măsurate sau înţelese uşor), sau de ieşire ( PEVE, ), dacă dimpotrivă, măsurarea sau perceperea lor este dificilă. Variabilele şi parametrii pot lua valori numerice sau logice. Deosebirea dintre variabile şi parametri constă în aceea că parametrii nu îşi schimbă valorile pe perioade mari de timp, în timp ce variabilele îşi schimbă valorile chiar pe intervale mici de timp. Scopul modelului este de a exprima variabilele şi parametrii de ieşire în funcţie de variabilele şi parametrii de intrare, cu eventuala satisfacere a unor condiţii de performanţă de către sistem (de ex. condiţii de optim). Unele VI pot fi aleatoare, caz în care şi unele variabile sau parametri de ieşire vor fi de asemenea aleatoare.
b) Relaţiile funcţionale constituie o altă categorie de elemente ale unui model matematic. Ele sunt de forma
( ) 0,,, =PEVEPIVIFi (adică implicite) şi pot fi la rândul lor de două tipuri: - ecuaţii, care sunt satisfăcute numai de anumite valori ale variabilelor sau parametrilor; - identităţi, care sunt satisfăcute de orice valori ale variabilelor şi parametrilor; ele exprima
condiţii de echilibru sau legi de conservare. Ecuaţiile si identităţile pot fi relaţii algebrice sau transcendente, diferenţiale sau integrale,
detrministe sau stochastice, etc. c) Caracteristicile operative constituie o altă categorie de elemente ce compun un model
matematic şi ele pot fi:
- ipoteze de lucru (referitoare la relaţiile funcţionale); - ipoteze statistice (referitoare la VI -aleatoare).
d) Tehnica de rezolvare este un alt element constitutiv al unui model matematic. Ea este o tehnică matematică ce realizează separarea elementelor de ieşire în funcţie de elementele de intrare, adică:
( ) ( )IIiE PVfPV ,, = . Cu alte cuvinte, tehnica de rezolvare a modelului exprimă sub formă explicită variabilele
şi parametrii de ieşire în funcţie de variabilele şi parametrii de intrare. Tehnicile matematice de rezolvare a modelelor sunt însă sărace atât ca varietate, cât şi ca performanţă. De exemplu, ecuaţiile modelului se pot reazolva numai dacă sunt liniare sau uneori pătratice, iar dacă sunt de grad superior ele se pot rezolva numai dacă au forme particulare. La fel, ecuaţiile diferenţiale sau cu derivate parţiale se pot rezolva cu metode matematice deductive numai în anumite cazuri particulare. De aceea în construcţia modelelor matematice se fac de multe ori ipoteze simplificatoare care permit aplicarea tehnicilor de care dispune matematica. (Acesta este scopul utilizării de către model a caracteristicilor operative!). Din aceste motive, modelarea matematică este aproximativă şi ea nu permite rezolvarea realistă a problemelor practice. Utilizarea calculatorului permite îmbunătăţirea performanţelor modelelor matematice prin utilizarea metodelor numerice. Dar şi în aceste condiţii modelele matematice nu pot descrie corect realitatea în toată complexitatea ei deoarece nu toate relaţiile dintre obiectele lumii reale se pot exprima prin formule matematice. Într-o atare situaţie modelul matematic trebuie completat cu descrieri care să imite anumite comportări ale lumii reale. Acestea se realizează prin descrieri algoritmice de tipul if-then-else- sau if-then- combinate cu alte structuri algoritmce (cicluri, secvenţe etc.). În acest fel, modelul matematic se completează, se extinde sub formă de algoritm şi devine model de simulare. Simularea măreşte deci mult posibilitatea de tratare realistă a problemelor aplicative. Construcţia unui model de simulare, care în fapt este un algoritm complex, dezvoltat pe scheletul unui model matematic, este o sarcină nu prea uşoară; o vom trata mai jos. Clasificări ale modelelor matematice. Mai întâi să vedem însă cum pot fi clasificate modelele matematice: (i) Clasificarea modelelor matematice după natura variabilelor utilizate de model:
continue/discrete; statice/dinamice (dacă timpul nu intervine sau dacă apare explicit ca variabilă a modelului); deterministe/stochastice (după cum nu conţin sau conţin măcar o variabilă de intrare ca variabilă aleatoare).
(ii) Clasificare topologică, după structura determinată de părţile în care se descopune modelul (când este cazul): cu o componentă/cu mai multe componente în serie , în paralel, în reţea.
Un model, fie el matematic sau de simulare, constituie de fapt o clasă de modele.
I.2. Construcţia unui model de simulare
Modelul de simulare (MS) presupune utilizarea calculatorului şi el se construieşte pe baza/scheletul unui model matematic; mai precis, el completează modelul matematic descriind unele relaţii prin algoritmi, deci în final MS este un algoritm. Acest algoritm trebuie să descrie corect evoluţia sistemului şi să permită efectuarea de experienţe (prin rulări pe calculator), care să înlocuiască experienţele ce ar trebui realizate asupra sistemului real.
Structura unui model de simulare
Construcţia unui MS utilizează două concepte de bază: ceasul simulării şi agenda simulării.
Ceasul simulării. Ceasul asigură eşalonarea corectă în timp a evenimentelor create de
model şi uneori ajută la implementarea condiţiei de terminare a simulării. El este de două feluri: a) ceas cu creştere constantă; b) ceas cu creştere variabilă.
Notă: Evenimentele corespund unor modificări în sistem adică modificari ale valorilor unor variabile care se calculează sau se generează prin instrucţiuni ale modelului (chiar şi dacă sunt aleatoare!).
Ceasul porneşte cu valoarea zero la începutul simulării. Dacă modelul se bazează pe ceasul cu creştere variabilă, atunci ceasul este crescut cu valoarea ce corespunde apariţiei primului eveniment următor; apoi programul prelucrează evenimentul, după care ceasul creşte din nou reluându-se ciclul simulării, până când ceasul atinge o valoare dată iniţial, care corespunde perioadei pe care se realizează simularea.
maxT
Ceasul cu creştere constantă presupune că de fiecare dată creşterea se realizează cu o cuantă de timp constantă; apoi se prelucrează toate evenimentele apărute pe intervalul de timp de lungime , după care se reia ciclul simulării. Simularea se termină de asemenea când ceasul atinge valoarea .
c
Tc
max
Terminarea simulării se poate realiza şi impunând condiţia ca modelul să prelucreze un anumit număr dat de evenimente de un tip precizat. (De exemplu, în cazul unui model de simulare a unui sistem de aşteptare, se poate impune condiţia ca simularea să se termine când s-a simulat servirea unui număr dat de clienţi).
A rămas oarecum neprecizat ce se înţelege prin prelucrarea unui eveniment. Acest fapt se înţelege uşor dacă precizăm şi conceptul de agendă a simulării.
Agenda simulării. Agenda se compune din elementele memorate de modelul de simulare. Variabilele modelului iau diverse valori pe parcursul simulării; de aici rezultă că pe parcursul simulării apar multe evenimente. Memorarea în totalitate a acestor evenimente, împreună cu caracteristicile lor, nu este nici recomandată dar nici necesară dacă simularea se realizează pe perioade mari de timp. De aceea, se memorează (în agendă) numai ceea ce este strict necesar. Evenimentele sunt de diverse tipuri; unele variabile descriu şi stări ale sistemului (sau ale unor componente ale acestuia). Evenimentele sunt create sau generate la momente de timp ulterioare valorii ceasului. De aceea agenda se compune din două părţi: agenda evenimentelor curente şi agenda evenimentelor viitoare . AEC AEV
Deci agenda , unde: AEVAECA ⊕=AEC = agenda evenimentelor curente (care au timpul de apariţie egal cu valoarea ceasului); iar AEV = agenda evenimentelor viitoare (care au timpul de apariţie ulterior ceasului).
Algoritmul simulării prelucrează deci numai evenimentele din ; prelucrarea unui eveniment înseamnă fie determinarea apariţiei unui nou eveniment (memorat în ) sau modificarea unei stări, fie distrugerea unui eveniment ("ştergerea'') lui din agendă. Prelucrările evenimentelor ţin seama şi de stările sistemului la acel moment.
AECAEV
Algoritmul simulării gestionează/actualizează agenda prin interacţiunea acesteia cu ceasul; într-un ciclu al simulării ceasul este acualizat, după care se selectează din agenda evenimentele care fac parte din şi se prelucrează aceste evenimente până când devine vidă. Atunci, ceasul este crescut din nou şi se reia ciclul simulării. Deci agenda simulării se modifică pe parcursul simulării conform următoarei relaţii de dinamică
AAECAEC
elimgen AEAEAA O⊕=
unde sunt evenimentele generate pe parcursul unui ciclu, iar sunt evenimentele eliminate cu ocazia prelucrării . (Semnele
genA elimAAEC ⊕ şi O se autexplică).
I.3. Concepte de bază în modelarea sistemelor
Prin sistem se înţelege, sub forma cea mai vagă, o mulţime de obiecte O interconectate prin intermediul unei relaţii R , adica Sistem = ( )RO , , unde
OOR ×⊂ Definiţia formală generală a unui sistem este: Definiţia 2. Un sistem este următoarea structură de mulţimi: ( )λδΣΩ= ,,,,,, YXTS unde: T = timpul de bază al sistemului (ceasul sistemului); X = mulţimea intrărilor; = mulţimea segmentelor de intrare (forma intrărilor!); segment =
,= =grafic=Ω
( )t1 Xt0 ,:ω ( 1,0 ttω ) ( )( ){ }10,, ttttt ≤≤ω . Se foloseşte de regulă notaţia
= mulţimea stărilor interne ale sistemului = memoria sistemului; Conceptul de stare este esenţial în modelarea sistemelor; el descrie structura internă (intimă!) a sistemului;
Σ
= funcţia de tranziţie a stărilor definită ca δΣΩ×Σδ : .
care este axioma semigrupului sau proprietatea de separabilitate a stărilor. Mulţimea (graficul) se numeşte traiectorie a stărilor; ( σω, )
Y = mulţimea ieşirilor (Sistemul este presupus deschis, intrările şi ieşirile fiind exterioare sistemului). Mulţimea Y conţine răspunsul sistemului la intrarea de forma ω , când la momentul iniţial sistemul se află în starea 0t σ . = funcţia de răspuns, de forma λ YTX ××Σλ : ; ( )tx,,σλ conduce la un segment de ieşire ce reprezintă forma răspunsului sistemului la intrarea x , la momentul când starea la momentul , , este
t0t tt <0 σ ;
Notă: din definiţie rezultă că pentru o stare iniţială σ , (la momentul ) când are loc intrarea de forma ω se realizează o traiectorie unică a stărilor. Cu alte cuvinte, traiectoria stărilor satisface relaţiile:
0t
( ) Σωσ 10, ,: ttSTRAJ , astfel încât ( ) σ=ωσ 0, tSTRAJ
( ) (( )ttSTRAJ ωσδ=ωσ ,, , ( )10 , ttt ∈∀ Proiecţia obţinută prin funcţia de tranziţie a stărilor compusă cu funcţia de
răspuns, este traiectoria de ieşire STRAJ
( ) YttOTRAJ 10, ,:ωσ
Dacă (ieşirea depinde numai de ( )σλ=λ σ ) atunci rezultă relaţia simplă: ( ) ( )( )tSTRAJtOTRAJ ωσωσ λ= ,, .
Nivele de reprezentare a sistemelor. Există mai multe nivele de reprezentare a sistemelor şi anume:
- Reprezentarea la nivel de comportare. Sistemul este o cutie neagră (black box), exprimat formal prin relaţia:
( ){ }Σ∈σ=ρΩ∈ωρω= ωσ ,,, ,OTRAJRs
Acest nivel de reprezentare este cel mai vag. El descrie numai relaţiile de intrare/ieşire ce
se pot observa din afara sistemului (care deci se comportă ca o cutie neagră). - Reprezentarea la nivel de structură de stare. La acest nivel se intră în structra internă
a sistemului, adică se stabilesc elementele acestuia : ( )λδΣΩ ,,,,,, YXT , definite ca mai sus.
- Reprezentarea modulară (ca structură compusă). Dacă sistemul este complex, atunci se identifică subsisteme ale acestuia precum şi interconexiunile dintre ele în sensul că ieşirile unor subsisteme sunt intări în alte subsisteme, intrările unor subsisteme fiind intrări ale sistemului şi ieşirile unor subsisteme fiind ieşiri ale sistemului.
Nivelele de reprezentare a sistemelor, pornind de la forma cea mai vagă şi continuând până la forma detaliată, legitimează Metodologia ''TOP DOWN'' de proiectare a sistemelor de orice fel, în particular, metodologia de proiectare a sistemelor informatice, şi în general metodologia de proiectare ierarhică, descendentă a programelor. 1.3.1 Realizări iterative ale sistemelor
Considerăm T o mulţime numită timp de bază, timp care poate fi discret când T⊆ Z sau de puterea continuului când T⊆ R. Numim interval de observaţie orice submulţime < t0, t1 > ⊂ T, unde < > ∈ { [ ], [ ), ( ], ( )}. Numim segment peste o mulţime arbitrară A, relativ la intervalul de observaţie < t0, t1 > orice funcţie ω:< t0, t1 > → A. Notăm dom ω = < t0, t1 > şi uneori ω< t0, t1 > în loc de ω. De asemenea, notăm cu (T,A) mulţimea tuturor segmentelor peste A relativ la intervale de observaţie din T. Spunem că ⊂ (T,A) este închisă prin segmentare (IPS) ⇔ oricare ar fi ω ∈ asfel încât dom ω = < t0, t1 > şi ∀
drept) şi ω<τ = ω⎮< t0,τ > (numit segment stâng). Numim obiect abstract [ 40,54] o mulţime A de perechi ordonate (ω, ρ) cu condiţia închiderii prin segmentare, în care ω şi ρ sunt legate prin relaţii şi/sau funcţii intrare/ieşire. Perechea (ω, ρ) se numeşte pereche intrare / ieşire; ω se numeşte intrarea lui A, iar ρ se numeşte ieşirea lui A.
Un obiect abstract este dat prin una sau mai multe relaţii sau funcţii intrare / ieşire, adică prin sisteme de ecuaţii diferenţiale sau cu derivate parţiale sau sisteme de ecuaţii cu diferenţe sau este dat prin algoritmi capabili să genereze toate perechile intrare / ieşire ale sale.
Numim sistem o mulţime de obiecte abstracte interconectate { A1, A2, … , An } în care o parte din intrările sau ieşirile unor obiecte pot coincide cu ieşirile respectiv intrările altor obiecte. Astfel un sistem se reduce tot la o mulţime de segmente legate prin relaţii şi/sau funcţii intrare/ieşire.
Se numeşte maşină secvenţială structura M = (X, Σ, Y, δ, λ) în care: X=alfabet de intrare Y=alfabet de ieşire Σ=o mulţime de stări λ : Σ → Y funcţia de ieşire a maşinii δ : Σ × X → Σ funcţia de tranziţie a stărilor.
Se numeşte realizare iterativă a sistemului A = {(ω, ρ)} structura G(A) = (T, X, G, Σ, δG, λ) în care:
Ω
T=timp de bază
X=mulţimea valorilor de intrare
Y=mulţimea valorilor de ieşire Σ=mulţimea stărilor sistemului ΩG=mulţimea generatoare a segmentelor de intrare λ : Σ → Y funcţia de ieşire δG : Σ × ΩG → Σ funcţia de tranziţie a stărilor ΩG ⊂ ( T, X )0 = {ω : <0,t> → X | t ∈ T }
δG este astfel încât extensia δ G : Σ × → Σ să verifice proprietatea de separabilitate a
stărilor (proprietatea compunerii), unde este închiderea la compunere a lui G care se
demonstrează uşor că este : = ∪ , unde
+ΩG+ΩG
iG
Ω+ΩG
0>i
Ω
= { ω · ω · … ·ω | ω ∈ iGΩ
1k 2k ik jk Ω G , j ∈ i,1 }
iar
δG (q,ω) , ω∈Ω G
δ G (q,ω)=
δ G (δG (q,x),ω′ ), ω = x ·ω′ ∈ , x∈+ΩG Ω G, ω′ ∈ ,q∈Σ +ΩG
)'),,(()',( ωωδδωωδ qq GGG = (proprietatea compunerii) Se numeşte model iterativ [1,40] al sistemului A = {(ω, ρ)} cu realizarea iterativă G(A) = (T, X,
ΩG , Σ , Y, δG, λ), structura SG(A) = (T, X, , Σ, Y, +ΩG δ G, λ), în care este închiderea la
compunere a lui ΩG, iar
+ΩG
δ G este extensia lui δG definită mai sus şi verifică proprietatea compunerii. Exemplu Fie maşina secvenţială M = (XM, ΣM, YM, δM, λM). Luăm realizarea iterativă GM = (T, X, ΩG , Σ , Y, δG, λ) în care:
T = Z, X = XM, Y = YM, Σ = ΣM, λ = λM, < > = [ ) este interval de observaţie şi deci < t0, t0 + 1> = [t0, t0 + 1) = { t0 } δG : Σ × ΩG → Σ este dată prin relaţia: δG (q,ω) = δM(q,ω(0))
ΩG = {ω| ω : < 0,1> → X} = {ω| ω : {0} → X} = {ωx| x ∈ X, ωx(0) = x}, mulţimea segmentelor de lungime 1.
Observaţie
Evident, ΩG ∼ X şi dacă X+ = ∪1≥k
kX
unde X = { x1 x2 … xk | xi ∈ X,i =k
k,1 }
avem ∼+ΩG
+X . Propoziţie
Fie M o maşină secvenţială şi δ M extensia funcţiei de tranziţie a stărilor la secvenţe de
simboluri, dată prin δ M : ΣM × → ΣM, +MX
),( αδ qM , α∈XM
δ M (q,α) =
δ M (δM (q,x),α′ ), α = xα′∈ , x∈ XM, α′ ∈ +MX +
MX
Atunci δ M (q, αα′ ) = δ M (δ M (q, α), α′ ) pentru ∀ q ∈ ΣM, ∀ α, α′ ∈ . +MX
Demonstraţie: Se face prin inducţie după lungimea secvenţei α, notată |α|. Propoziţie
Fie M o maşină secvenţială definită ca în exemplul de mai sus şi GM o realizare iterativă a
sa. Atunci, dacă ω ∈ şi dom ω = < 0, n > avem
+ΩG
))1( −nω)...1()0(,(),( = qq MG ωωδωδ Demonstraţie Se face prin inducţie după lungimea lui ω, notată l(ω). Propoziţie
Fie M o maşină secvenţială cu XM, ΣM, YM spaţii vectoriale de dimensiune finită peste un corp K, λM : ΣM → YM dată prin λM(q) = Cq iar δM : ΣM × XM → ΣM , dată prin relaţia δM(q,x) = Aq + Bx, unde A, B, C sunt transformări liniare peste spaţii vectoriale de dimensiune finită, A:ΣM → ΣM, B:XM → ΣM, C:ΣM → Y. Atunci modelul iterativ extins
)(MGS este modelul unui sistem liniar numit sistem discret invariant în timp specificat de A, B, C, K, în care G(M) este realizarea iterativă a lui M construită ca în exemplul de mai sus.
Demonstraţie
Se construieşte din aproape în aproape funcţia de tranziţie a stărilor modelului iterative. Se obţine:
∑−
=
−−+=1)(
0
1)()( ))((),(ω
ωω ωωδl
i
illG iBAqAq
relaţie care se demonstrează prin inducţie după lungimea lui ω. 1.3.2 Modelarea sistemelor cu evenimente discrete
În acest paragraf vom prezenta un exemplu de modelare a sistemelor cu evenimente discrete cu scopul de a evidenţia o metodă de lucru şi de a reliefa dificultatea unei astfel de modelări. Am văzut în introducere că modelul este o reprezentare izomorfă a realităţii, reprezentare care oferă o imagine intuitivă şi riguroasă a fenomenului studiat şi în plus facilitează descoperirea unor conexiuni şi legităţi care nu pot fi observate la o analiză sumară. Reprezentarea unui model nu trebuie neapărat să fie perfectă, ci să contribuie la clarificarea gândirii, la culegerea şi înregistrarea de informaţii necesare elucidării consecinţelor presupunerilor făcute asupra sistemului modelat. Modelul unui sistem este o reprezentare simplificată a acestuia orientată spre studiul anumitor tipuri de acţiuni ale sistemului real. Definitie Se numeşte sistem cu evenimente discrete (SED) structura
M = (XM, SM, YM,δ M , λM, t ) în care: XM = mulţimea evenimentelor externe SM = mulţimea stărilor secvenţiale ale sistemului YM = mulţimea valorilor de ieşire
δ M : ΣM × (XM ) → SM este funcţia de cvasitranziţie a stărilor, dată prin: { }∪ Φ
δ M((s,e),φ) = δφ(s), δφ : SM → SM fiind funcţie de tranziţie autonomă.
ΣM = {(s,e) | s ∈ SM , 0 ≤ e ≤ t (s)} λM : ΣM → YM este funcţia de ieşire a sistemului.
t : SM → [0, ∞) este funcţia de avans a timpului. Pentru s ∈ S, t (s) reprezintă timpul maxim cât sistemul rămâne în starea s.
Neapariţia evenimentelor externe a fost notată cu φ. În cele ce urmează modelăm iterativ un sistem cu evenimente discrete. Vom căuta mai întâi o realizare iterativă G(M) = (T, X, ΩG, Σ, Y, δG,λ), în care vom lua T = Z , intervalul de observaţie de forma [ ) şi X = XM ∪ {φ} . Pentru aceasta vom lua ca mulţime generatoare a segmentelor de intrare mulţimea ΩG = Ωφ ∪ ΩX unde Ωφ = {φτ | τ > 0, φτ : [ 0,τ) →{φ}} este mulţimea segmentelor de lungime finită fără evenimente externe iar ΩX = {xτ | τ > 0, xτ : [ 0,τ) → XM ∪ {φ}, xτ (0) = x, xτ (t) = φ, ∀ t ∈ (0,τ)} este mulţimea segmentelor de lungime finită pentru care evenimente externe pot apare doar la momentul t = 0. Vom defini δG : Σ × ΩG →Σ dată prin
⎩⎨⎧
>+≤++
=−+ )(),),0),(((
)(),,()),,((
)( stesstees
essteG
G τφδδττ
φδτφ
τ
)),0),),,(((()),,(( ττ φδδδ xesxes MGG =
Observaţie:
În cele ce urmează vom nota pentru uşurinţă δG (s,e, rφ ) în loc de δG ((s,e), rφ ) şi analog δG (s,e, xr) în loc de δG ((s,e), xr). Acest model general de sisteme cu evenimente discrete înglobează sistemele de aşteptare, sistemele informatice, reţelele de calculatoare etc.
Propoziţie
Fie δφ : SM → SM funcţie de tranziţie autonomă într-un SED şi δ φ : SM× Ζ+ → SM o extensie a sa dată prin
δ φ (s,0) = s
δ φ (s,n+1) = δφ(δ φ (s,n)), n ≥ 0
δ ∅ (s,n) reprezintă starea la care ajunge sistemul pornind din s după ce face n tranziţii. Atunci:
a) δ φ (s,n) = (s), ∀ n ∈N* , ∀ s ∈ SM nφδ
b) δ φ (s,n+p) = δ φ (δ φ (s,n),p), ∀ p,n ∈ Ζ+, ∀ s ∈ SM Demonstraţie
a) Se face prin inducţie după n ≥ 1. b) Se foloseşte a) sau se demonstrează prin inducţie după p ≥ 1, folosind următoarea
observaţie. Observaţie
Luând n = 0 în definiţia lui δ φ se obţine δ φ (s,1) = δφ(δ φ (s,0)) = δφ(s). Definiţie Fie M un SED şi σ : SM × Z+ → [0,∞) dată prin: σ(s,0) = 0,
∑ ≥∀=−
=
1
01)),,((),(
n
inistns φδσ
Funcţia σ(s,n) dă timpul total folosit de sistem pentru a face n tranziţii după ce a plecat din starea s, adică timpul în care sistemul stă în cele n stări atinse prin n tranziţii succesive, timp măsurat din momentul plecării din starea s (acest timp se poate măsura şi din momentul intrării în starea s
luând ).1)),,((),(1∑=
≥∀=n
i
nistns φδσ
Propoziţie Fie M un SED şi σ definită ca mai sus. Atunci: a) Funcţia us : Z+ → [0,∞) dată prin us(n) = σ(s,n), n ≥ 0 este nedescrescătoare pentru ∀ s ∈ SM.
b) σ(s,n+m)= σ(s,n)+ σ( φδ (s,n),m), ∀ s ∈ SM, ∀ m,n ∈ Ζ+. Demonstraţie Pentru un SED M notăm ms,e,τ = max {n∈ N | σ(s,n) < e + τ} care reprezintă numărul maxim de tranziţii făcute de sistem în intervalul de timp <0,e + τ> plecând din starea s, atunci când acest maxim există.
Observaţie
a)Dacă pentru orice n∈ N, σ(s,n) < e + τ atunci {n∈ N | σ(s,n) < e + τ}= N şi deci nu există ms,e,τ. b)Dacă ∃ j∈ N astfel încât j∉ {n∈ N | σ(s,n) < e + τ} atunci pentru orice m ≥ j, m ∉ {n∈ N | σ(s,n) < e + τ} şi astfel {n∈ N | σ(s,n) < e + τ} este finită, de unde rezultă că în acest caz există ms,e,τ.
Propoziţie Fie M un SED. Atunci:
a)Dacă ∃ ms,e,τ avem σ(s, ms,e,τ) < e+τ ≤ σ (s, ms,e,τ+1) b)Dacă ∃ s ∈ SM pentru care Lns
n=
∞→),(limσ , există şi este finită, atunci ∀ e,τ astfel încât
e + τ > L, nu există ms,e,τ. c)Pentru ∀ s ∈ SM, ∃ ms,e,τ, ∀ e,τ astfel încât e + τ > 0 ⇔ ∞=
∞→),(lim ns
nσ
Demonstraţie Propoziţia 2.2.6. [40] Fie M un SED pentru care există ms,e,τ ∀ s ∈ SM, ∀ e,τ ∈ Ζ+.
Atunci σστδ ττφ ),,(),,( ,,,, eses msemsm −+ =
=max {n ∈ N | σ(s,n + ms,e,τ)< e + τ + σ}. Demonstraţie Propozitie Fie M un SED pentru care există ms,e,τ, ∀ s ∈ SM, ∀ e,τ ∈ Ζ+. Atunci
σστδτστττφ ),,(),,(,,,,
,,,, eses msemseses mmm−++ +=
Demonstraţie Propoziţie Într-un SED avem următoarele proprietăţi:
a) ms,e,τ = 0 ⇔ e + τ ≤ t (s).
b) τ +e > t (s) ⇒ 1)(,0),(,, +=−+ steses mm τδτ
φ
Demonstraţie a)Implicaţia de la stânga la dreapta se obţine punând ms,e,τ = 0 în 2.2.5.a) şi folosind relaţia
σ (s,1) = t (s) care rezultă luând ∑−
=
=1
0)).,((),(
n
iistns φδσ
Aceasta reprezintă timpul total necesar sistemului pentru a face n tranziţii din momentul intrării în starea s. Implicaţia de la dreapta la stânga se demonstrează prin absurd folosind forma de mai sus a lui σ. Propoziţie Fie M un SED pentru care există ms,e,τ, ∀ s ∈ SM, ∀ e,τ ∈ Ζ+.
Propoziţie Fie M un SED pentru care există ms,e,τ, ∀ s ∈ SM, ∀ e,τ ∈ Ζ+. Atunci )),,((),( στστ φφδδφδ qq GGG =+ , ∀ q = (s,e) ∈ Σ, ∀ τ,σ ∈ Ζ+. Demonstraţie Se evaluează cei doi membri ai relaţiei folosind 2.2.9. şi se obţine egalitatea cerută. Propoziţie Fie M un SED pentru care există ms,e,τ, ∀ s ∈ SM, ∀ e,τ ∈ Ζ+. Atunci )),,((),( στστ φδδδ xqxq GGG =+ ,∀ q ∈ Σ, ∀ τ,σ ∈ Ζ+. Demonstraţie Definiţie Un SED M se numeşte viabil ⇔ ∀ s ∈ SM, ∞=
∞→),(lim ns
nσ
Teorema
Fie M un SED viabil. Atunci structura G(M) = (T, X, ΩG , Σ , Y, δG, λ) este o realizare iterativă a lui M, iar
),,,,,,()( λδGGMG YXTS ΣΩ= +
este un model iterativ pentru sistemul cu evenimente discrete M în care: T = Z+, cu intervale de observaţie de forma : < > = [ ) . X = XM ∪ {φ}, φ∉XM
∑ = ∑M = {(s,e) | s ∈ SM, 0 ≤ e ≤ t (s)}
Y = YM λ = λM ΩG = Ωφ ∪ ΩX
δG : ∑ × ΩG → ∑ dată prin :
)),0),),,(((()),,((
)),(),,(()),,(( ,,,,
ττ
ττφτ
φδδδ
στδφδ
xesxes
msemses
MGG
esesG
=
−+=
Demonstraţie Din definiţia lui ΩG se observă că aceasta poate fi luată ca mulţime generatoare a segmentelor de intrare. Folosind propoziţiile de mai sus se demonstrează că δG verifică proprietatea de separabilitate a stărilor.
I.3.3. Metodologia de realizare a experimentelor de simulare (Metodologia simulării)
Etapele realizării unui experiment de simulare sunt: 1. Formularea problemei; constă în a preciza întrebările la care trebuie să răspundă modelul şi
a preciza domeniul lumii reale ce trebuie analizat. Aici se precizează şi forma răspunsului la întrebări (de ex. dacă se vor produce grafice, tabele sau chiar texte scrise).
2. Realizarea unor experimente preliminare; în această etapă se stabilesc, pe baza observaţiilor sau datelor culese din lumea reală, sau existente, referitoare la aceasta (istorice), variabilele şi parametrii şi care din acestea sunt de intrare sau de ieşire.
3. Prelucrarea (interpretarea) primară a datelor preliminare; acum se disting variabilele aleatoare, se estimează parametrii şi se testează ipoteze statistice. (Statistica matematică joacă un rol important în simulare).
4. Formularea unui model matematic preliminar, incomplet; în această etapă se precizează relaţii funcţionale, ipoteze de lucru (care să concorde cu datele existente, colectate) şi se indentifică ce relaţii nu pot fi exprimate matematic şi care sunt dificultăţile ce ar trebui înlăturate pentru a răspunde la întrebările formulate.
5. Evaluarea modelului etapă de decizie; se urmăreşte să se evalueze complexitatea modelului, dacă el poate răspunde în timp real şi complet la întrebări; în această etapă se pot eventual revizui unele din etapele precedente făcându-se simplificări sau completări ale modelului matematic propus.
6. Construcţia modelului de simulare (sub formă de algoritm detaliat); modelul se construieşte conform celor precizate anterior urmărindu-se ca el sa fie general. La construcţia algoritmului simulării se va avea în vedere şi cum se va realiza programarea: într-un limbaj evoluat sau înt-un limbaj specializat de simulare. Există un număr mare de limbaje de simulare toate având la bază, mai mult sau mai puţin, filozofia DEVS. Un astfel de limbaj (ce va fi prezentat în secţiunea următoare) este GPSS= General Purpose Simulation System. O versiune românească a acestui limbaj este SIMUB (limbajul de SIMulare al Universităţii din Bucureşti).
7. Validarea modelului este de asemenea o etapă importantă. În afară de validarea formală (testare sintactică şi formală a programului), trebuie să se decidă dacă modelul rezolvă corect problema formulată. Acest lucru se realizează prin compararea rezultatelor simulării fie cu rezultate practice cunoscute, fie prin compararea cu soluţia obţinută cu ajutorul unui model matematic. (De ex. în teoria cozilor se cunoaşte soluţia matematică a modelului în anumite ipoteze restrictive; se va obţine soluţia prin simulare pe baza unei selecţii de volum mare a variabilelor de ieşire. Dacă soluţia simulată se apropie de soluţia matematică, atunci, în virtutea legii numerelor mari, rezultă corectitudinea modelului de simulare, acesta fiind deci valid şi în ipoteze diferite de cele cosiderate la validare).
8. Planificarea experienţelor de simulare este următoarea etapă necesară. Aşa cum am precizat, simularea este un experiment realizat cu calculatorul pe baza unui model ce descrie sistemul real. Orice experiment (cu caracter statistic ca în cazul simularii) trebuie să se desfăşoare pe baza unui plan experimental. Pentru planificarea experimentelor de simulare se folosesc metode ale statisticii matematice (aşa-zisele experimental designs).
9. Prelucrarea şi interpretarea experienţelor de simulare este o etapă la fel de importantă ca şi cele precedente. Programul de simulare rulat pe calculator produce de regulă valori de selecţie asupra variabilelor de ieşire, dar care nu definesc selecţii statistice în sensul obişnuit (nu sunt selecţii bernouliene deoarece valorile de selecţie sunt de regula dependente). De aceea prelucrarea experimentelor de simulare se face cu mijloace statistice din cele mai sofisticate (prelucrarea seriilor dinamice). Toate limbajele de simulare posedă un minim de facilităţi pentru prelucrări statistice simple (calcule de medii, dispersii, coeficienţi de corelaţie, histograme etc,).
I.4. Generalităţi despre limbajul GPSS
Limbajele de simulare sunt în acelaşi timp limbaje de programare şi limbaje de modelare; ele implementează elementele esenţiale ale simulării: manipularea ceasului şi gestionarea memoriei. Utilizatorul are grija descrierii evenimentelor şi prelucrarea lor.
Vom face o foarte scurtă prezentare a limbajului GPSS (General Purpose System Simulator), urmând ca cititorul interesat să recurgă la un manual detaliat sau să folosească
facilitatea help existentă în implementarea GPSS/PC. Acest limbaj a fost dezvoltat pentru prima data de firma IBM la inceputul anilor '60.
Entităţile limbajului GPSS. Limbajul GPSS se compune din 16 tipuri de entităţi (elemente de abstractizare). La
fiecare entitate se asociază un numar de proprietăţi sau atribute în majoritatea lor adresabile intern (de către GPSS), dar unele sunt adresabile şi de către utilizator; atributele sunt: standard numerice (numere) sau standard logice (valori logice).
Entităţile de bază sunt: 1) Blocuri (entităţi care descriu {activităţi); 2) Tranzacţii (elemente circulante); ele sunt create (printr-un bloc special GENERATE) şi
circulă prin model (sistem!), ca urmare a acţiunii altor blocuri. Blocurile au asociate caracteristici numerice sau logice şi posedă argumente necesare
descrierii activităţilor. Tranzacţiile posedă un număr de parametri standard dar şi parametri introduşi de utilizator. Prin programul de simulare utilizatorul poate accesa parametrii tranzacţiilor.
Entităţi de echipament: 3) Staţii de serviciu sau facilităţi (ele corespund subsistemelor cu o componentă care tratează
de regulă o singura tranzacţie!); 4) Multistaţii de serviciu sau depozite; acestea tratează de regulă mai multe tranzacţii (de ex.
liftul, autobuzul, pot transporta mai mulţi clienţi consideraţi ca tranzacţii etc.); 5) Comutatorii logici sunt variabile logice care permit utilizatorului să realizeze ramificarea
după dorinţă a fluxului de execuţie în programul de simulare. Entităţi de calcul ce permit efectuarea (limitat!) a unor calcule:
6) Variabile aritmetice; permit evaluarea unor expresii aritmetice, iar rezultatul este memorat de o variabilă aritmetică;
7) Variabile booleene; permit evaluarea unor expresii booleene, iar rezultatul este memorat de o variabilă booleană;
8) Funcii descrise prin tabele sau prin segmente de dreaptă (liniare pe porţiuni); 9) Funcţii analitice descrise prin expresii mai complicate (ele există numai în SIMUB şi au
rolul de a facilita simularea diverselor variabile aleatoare prin metoda inversă); Entităţile statistice, permit colectarea unor statistici sau estimaţii privind VE sau PE ;
printre acestea se remarcă: 10) Cozile care sunt entităţi statistice ce memorează tranzacţiile întârziate în sistem; 11) Tabele de frecvenţe care descriu histograme ale VE ; 12) Tabele de frecvenţă bidimensionale sau tabele de contingenţă (acestea sunt disponibile
numai în SIMUB); Entităţi de referinţă care memorează anumite informaţii pe care le doreşte utilizatorul şi anume: 13) Cuvinte de salvare care memorează valori ce corespund câte unui cuvânt de memorie al
calculatorului; 14) Matrice de cuvinte de salvare care memorează matrici;
Entităţi de tip lanţ, care sunt de două tipuri: 15) Lanţuri ale utilizatorilor în care se pot depune sau scoate tranzaţii după dorinţa
utilizatorului (Există şi lanţuri ale sistemului, manipulate intern de către GPSS, care nu pot fi accesate de utilizator).
16) Grupuri care separă tranzacţii (cu anume scopuri de prelucrare dorite de utilizator). NOTĂ. Toate entităţile trebuie definite/identificate şi declarate de către utilizator la
începutul programului. Pentru fiecare entitate se alocă memorie, la instalarea limbajului se alocă şi un număr maxim de entităţi care se pot utiliza într-un program.
Entităţile se invocă prin instrucţiuni cu structura fixă. (Limbajul este un interpreter).
Limbajul permite descrierea modulară a unui model de simulare general; cel mai important este modulul principal (ciclul de bază). Ieşirile sunt standard; se afişează/scriu valori ale atributelor, statistici, tabele de frecvenţă etc. Limbajul GPSS este un interpreter, ceea ce înseamnă că fiecare instrucţiune scrisă corect (şi recunoscută ca atare de către GPSS) se execută imediat.
Structura instrucţiunii GPSS.
Instrucţiunea GPSS are un format fix.
Nume simbolic (etichetă)
Numele blocului (cuvânt cheie)
Argumente/parametrii (separate prin virgulă)
Opţional pentru referiri
Verb imperativ care desemnează acţiunea
blocului A, B, C, …
1) Blocurile descriu acţiuni (sunt entităţi active spre deosebire de tranzacţii care sunt entităţi pasive, în sensul că ele sunt mişcate prin model ca urmare a acţiunii blocurilor). Tipurile de blocuri sunt: a) Blocuri de acţiune: SEIZE/RELEASE (tranzacţia curentă ocupă sau eliberează o
facilitate); PREEMPT (tranzacţia intrată în acest bloc poate ocupa o staţie de serviciu indicată de parametrul-argument sau o poate prelua dacă este ocupată de altă tranzacţie); ENTER/LEAVE (ocupă sau eliberează o poziţie dintr-o multistaţie); QUEUE/DEPART (intră în coadă sau pleacă din coadă); LINK/UNLINK (introduce sau scoate tranzacţia dintr-un lanţ al utilizatorului); SPLIT (descompune tranzacţia în mai multe tranzacţii care formează o familie); ASSEMBLE/GATHER (unifică tranzacţiile dintr-o familie fără, sau cu păstrarea atributelor de bază din familie); MATCH (sunt 2 blocuri conjugate; ele sincronizează deplasarea tranzacţiilor care întâlnesc aceste blocuri); ADVANCE (întârzie tranzacţiile; simulează de ex. durata de serviciu a tranzacţiei); BUFFER (înseamnă prelucrarea tranzacţiilor în ordinea priorităţilor); JOIN/REMOVE (introduc sau scot tranzacţia într-un grup); SCAN (verifică dacă există o tranzacţie în grup cu o anumită proprietate);
b) Blocuri de creare şi distrugere de tranzacţii: GENERATE/TERMINATE (generarea se face cu generatori specializaţi de numere aleatoare, tranzacţia putând fi prevăzută cu priorităţi);
c) Blocuri de control logic al tranzacţiilor: TEST (controlează dacă 2 atribute ale tranzacţiei satisfac o condiţie); TRANSFER (asigură transfer condiţionat sau nu al fluxului de execuţie la blocul indicat prin argument); GATE (modifică drumul tranzacţiilor în funcţie de condiţii referitoare la atribute ale facilităţilor, multistaţiilor, comutatorilor logici sau tranzacţiilor); EXAMINE (modifică fluxul în funcţie de apartenenţa la un grup); LOOP (repetă execuţia de la blocul menţionat ca argument până la blocul LOOP);
d) Blocuri de modificare a caracteristicilor tranzacţiilor şi a valorilor unor entităţi de referinţă: ASSIGN (atribuie valori numerice parametrilor tranzacţiilor şi/sau le modifică; când se generează o tranzacţie poate fi prevăzută cu ''locuri'' pentru un număr de maximum 12 parametri referibili de către utilizator; numărul parametrilor poate fi standard sau precizat de utilizator la generare); INITIAL (iniţializează cuvintele sau matricile de cuvinte păstrate); PRIORITY (modifică priorităţile; lucrează în legătură cu blocul BUFFER); LOGIC (poziţionează pe true sau false un comutator logic ce poate juca rol de semafor); MARK (utilizat pentru a marca într-un cuvânt special asociat fiecărei tranzacţii, valoarea ceasului); SAVEVALUE/MSAVEVALUE (memorează într-un cuvânt/matrice valoarea unui atribut standard numeric precizat); COUNT (determină numărul de entităţi ce satisfac o anumită condiţie şi-l memoreaza într-un parametru al tranzacţiei specificat ca argument al lui COUNT); SELECT (selectează prima entitate
dintr-o gamă de entităţi ce satisfac o anumită condiţie); ALTER (modifică condiţionat sau nu, parametri sau prioritatea uneia sau mai multor tranzacţii dintr-un grup putând modifica opţional şi drumul tranzacţiei); HELP (permite includerea unor proceduri de calcul ale utilizatorului; este un bloc pretenţios căci necesită interfaţa între GPSS şi limbajul în care este scrisă procedura; în SIMUB se foloseau subrutine FORTRAN);
e) Blocuri pentru obţinerea de statistici: TABULATE (pentru construirea histogramei); (în SIMUB se foloseşte şi BTABULATE pentru construirea tabelelor de contingenţă).
f) Blocuri pentru listări: PRINT (permite scrierea parţială a unor statistici; cele mai multe statistici se scriu automat de GPSS la sfârşitul simulării).
2) Tranzacţiile se generează cu GENERATE şi se distrug cu TERMINATE. La generare
tranzaţiile primesc automat nişte parametri interni şi un număr de parametri (limitat de exemplu la 100) declaraţi şi care pot fi accesaţi de utilizator; în funcţie de implementare, un program GPSS poate utiliza un număr maxim de tranzacţii (acesta este un dezavantaj al limbajelor de simulare datorat faptului că limbajul gestionează agenda, facilitând astfel detaliile de programare pe care într-un limbaj evoluat le rezolvă utilizatorul).
3) Staţiile de serviciu pot fi ocupate de o singură tranzacţie la un moment dat. 4) Depozitele sau multistaţiile au o capacitate declarată în prealabil şi în ele pot intra mai multe
tranzacţii (cât permite capacitatea). 5) Comutatorii logici sunt de fapt variabile booleene, ce trebuie iniţializate (precizează condiţii
satisfăcute de ''echipamente''). Toate entităţile de echipament sunt identificate printr-un număr (întreg). Ele posedă atribute standard numerice sau logice. 6) Variabilele au cuvântul cheie VARIABLE, au un nume şi o expresie aritmetică (de
dimensiune limitată) care conţine atribute standard numerice, cuvinte de salvare, parametri etc.
7) Variabilele booleene BVARIABLE, sunt construite analog (de către utilizator) folosind atribute standard logice, inclusiv comutatori logici).
8-9) Funcţia (FUNCTION) desemnează o funcţie de o variabilă dată printr-o listă sau o funcţie liniară precizându-se argumentul prin care este referită.
10) Coada este o entitate statistică pentru care se reţine automat lungimea medie şi lungimea maximă care se listează la sfârşitul simulării.
11-12) Tabela de frecvenţă se defineşte la începutul programului cu TABLE şi reprezintă o histogramă a unei variabile de ieşire din model, de regulă atribut numeric. La definire se precizează argumentul tabelat şi numărul de clase de frecvenţe.
13-14) Cuvintele sau matricile de cuvinte păstrate sunt invocate prin INITIAL, SAVEVALUE sau MSAVEVALUE şi ele se folosesc pentru a reţine anumite valori pe parcursul simulării. Sunt referite printr-un număr de identificare.
15) Lanţurile sistemului sunt: lanţul evenimentelor curente (LEC), lanţul evenimentelor viitoare (LEV), lanţuri ale tranzacţiilor întrerupte, lanţuri ale tranzacţiilor în aşteptare pentru sincronizări, lanţuri de întârziere (asociate echipamentelor). Orice tranzacţie activă la un moment dat se găseşte într-un lanţ. Tranzacţiile care se deplasează în sistem se găsesc în lanţul evenimentelor curente. Mecanismul de actualizare a lanţului evenimentelor curente şi al celor viitoare este cel cunoscut (descris prin relaţia dintre AECA, şi AEV ). Tranzacţiile terminate (prin blocul TERMINATE) sunt distruse; blocul START precizează printr-un argument al său câte tranzacţii se vor prelucra (termina). Prelucrarea tranzacţiilor din LEC înseamnă transferul lor în alt lanţ sau distrugerea. Tranzacţiile întârziate ( de ex. prin ADVANCE) sunt introduse în LEV.
Lanţurile sistemului sunt controlate automat de către sistemul GPSS. Lanţurile utilizatorului (LU) sunt definite de acesta şi ele (sau atribute ale lor) pot fi
referite în programul GPSS. Tranzacţiile pot fi introduse sau scoase din LU prin LINK şi
UNLINK şi ele se folosesc la implementarea unor discipline de serviciu (în cazul prelucrării tranzacţiilor din cozi).
16) Grupurile oferă utilizatorului un instrument de clasificare a tranzacţiilor care au anume
proprietăţi. De ex. toate staţiile care la un moment dat depăşesc un anumit grad de utilizare a unei staţii pot fi făcute membre ale unui grup.
Un program GPSS este o listă ordonată de BLOCURI ale căror argumente se referă la diverse entităţi.
I.5. Exemple de programe GPSS.
I.5.1. Model de simulare pentru un sistem de aşteptare cu o staţie
Programul GPSS pentru acest model este: 001 GENERATE 10,2 ;Sosiri aleatoare la 210 ± minute 003 SEIZE FRIZER ;Ocupă staţia de serviciu "FRIZER" 005 ADVANCE 12,3 ;Durata servirii 312 ± minute 007 RELEASE FRIZER ;Eliberează staţia 009 TERMINATE 1 ;Ieşire din sistem un client
Comentariu. Cifrele din faţa tabelului reprezintă etichetele instrucţiunilor/blocurilor, iar textele de la terminarea liniilor reprezintă comentarii. După denumirile blocurilor apar parametri acestora.
Modelul nu este complet deoarece în faţa blocului SEIZE pot să apară aşteptări de aceea înaintea acestui bloc trebuie introdus un bloc QUEUE şi după el un bloc DEPART ca mai jos (etichetele indică locul unde se introduc aceste blocuri respectând ordinea).
002 QUEUE COADA1 ;Tranzacţia ocupă COADA1 003 SEIZE ……… 004 DEPART COADA1 ;Tranzacţia pleacă, se actualizează
aşteptările Pentru a prelucra 200 tranzacţii se va da la început comanda:
000 START 200 La sfârşitul simulării apelând procedura GPSSREPT.EXE se va afişa un raport final
privind simularea care conţine: - durata totală a prelucrării tranzacţiilor (durata simulării); - gradul de ocupare a staţiei de serviciu (media şi maxima); - durata medie a aşteptărilor etc.
I.5.2. Model de simulare pentru un sistem de aşteptare cu
staţii paralele
Modelul de mai jos simulează un sistem de aşteptare cu 3 staţii paralele (ghişee ale unei bănci) în care serviciul se realizează cu disciplina FIFO (First In First Out) adică serviciul se execută în ordinea sosirii clienţilor şi nu se consideră priorităţi ale unor clienţi. Clientul va trece pe la staţia multiplă de servire (depozit), dar el va ocupa o singură staţie care-l serveşte (STORAGE 3). În final el va trece şi la staţia CASIER unde primeşte un alt serviciu (primeşte sau plăteşte banii pe baza formelor elaborate de ghişeul parcurs anterior!). Desigur, pentru executarea modelului care urmează ar trebui adăugată comanda START corespunzătoare şi apelat în final modulul de ieşire GPSSREPT.EXE. Modelul este dat de programul:
003 ENTER GHISEE ;Ocupă un loc în GHISEE (dacă există!) 004 DEPART COADA ;Clentu pleacă din COADA (spre servire!) 005 ADVANCE 15,3 ;Clientul este servit (durata de servire!) 006 LEAVE GHISEE ;Se eliberează un loc (un ghişeu) 007 SEIZE CASIER ;Ocupă loc la CASIER 008 ADVANCE 2 ;Servirea clientului este 2 minute 009 RELEASE CASIER ;Eliberează casierul (clientul pleacă) 010 TERMINATE 1 ;Plecare client
I.5.3. Model cu staţii paralele şi preferinţe
Se modelează serviciul la o frizerie unde lucrează 3 frizeri (trei staţii de serviciu!):
Figaro, Gică şi Tică. Clienţii care sosesc sunt primiţi în proporţie de 60 % de Figaro (care este preferat!) iar restul de 40 % sunt serviţi de ceilalţi doi, dar din aceştia 50 % merg la Gică şi 50 % merg la Tică. Modelul nu conţine comenzile de START şi cele de ieşire.
001 GENERATE 6,1 ;Sosesc clienţii 002 TRANSFER 4,,ALTII ;60 % clienţi merg la Figaro 003 QUEUE COADAFIG ;Intră în coada lui Figaro 004 SEIZE COADAFIG ;Clientul se duce la Figaro; 005 DEPART COADAFIG ;Se culeg date despre această coadă 006 ADVANCE 8,1 ;Clientul este servit 007 RELEASE FIGARO ;Figaro devine liber 008 TRANSFER ,CASA ;Clientul merge să plătească ALTII TRANSFER ,LAGICA ;Tică şi Gică sunt la fel 010 SEIZE TICA ;50 % din clienţi merg la Gică 011 ADVANCE 12,2 ;Tică serveşte mai încet ca Figaro 012 RELEASE TICA ;Tică este eliberat 013 TRANSFER ,CASA ;Clientu plăteşte LAGICA SEIZE GICA ;50 % şi 40 % clienţi merg la Gică 015 ADVANCE 12,2 ;Gică serveşte ca şi Tică 016 RELEASE GICA ;Se eliberează Gică CASA SEIZE CASIERA ;Se ocupă staţia casiera 018 ADVANCE 1 ;Plata într-un minut 019 RELEASE CASIERA ;Se eliberează casiera 020 TERMINATE 1 ;Tranzacţia pleacă
Notă. Cele trei staţii sunt folosite explicit, nu ca multistaţie. Folosim coada numai la Figaro căci numai acolo (presupunem!) ne interesează statistici finale.
În programele GPSS etichetele sunt opţionale şi ele pot fi mnemonice (care denumesc obiecte modelate ca de ex. LAGICA sau ALTII), sau pot fi constante întregi care sugerează ordinea instrucţiunilor; sunt obligatorii etichetările blocurilor la care se referă alte blocuri (vezi TRANSFER mai sus).
Parametri unor entităţi pot desemna nume simbolice (de ex. staţia GICA sau coada COADAFIG).
CAPITOLUL II. NUMERE ALEATOARE
II.1. Noţiuni introductive
Vom aminti mai întâi câteva noţiuni de bază. Presupunem că X este o variabilă aleatoare şi fie funcţia sa de repartiţie. Densitatea de repartiţie, când aceasta există, este derivata funcţiei de repartiţie, adică
( ) 0=∞−F ( ) 1=∞+ ( )F ( )bFaF ≤ dacă ba < . Între ( )xF şi ( )xf are loc relaţia:
( ) ( )∫∞−
=x
dxxfxF
De obicei X reprezintă o caracteristică a unei mulţimi de obiecte care variază aleator de la un obiect la altul; acea mulţime se numeşte populaţie statistică. Dacă luăm la întâmplare un număr
de obiecte din populaţie şi considerăm valorile ale caracteristicii n nXXX ,...,, 21 X ce corespund acestor obiecte, spunem că aceste valori determină o selecţie de volum asupra lui nX . În statistica matematică selecţia este considerată ca fiind o mulţime de
variabile aleatoare independente şi identic repartizate ca şi nXXX ,...,, 21
X (aceasta numindu-se selecţie bernouliană). Independenţa a două variabile aleatoare se defineşte astfel: X este o variabilă aleatoare independentă de Y dacă ( ) ( ) ( )yYPxXPyYxXP <<=<< , sau, dacă notăm
funcţia de repartiţie comună a variabilelor ( ) ( yYxXyx <<= ,, )PF X şi Y şi notăm , ( ) )xX <= ( )yF(XPF1 (Y )yP <−2 funcţiile de repartiţie marginale, ale lui X respectiv Y ,
atunci condiţia de independenţă se scrie ( ) ( ) ( )yFxFyxF 21, = (2.1.)
Să observăm că se face distincţie între variabila aleatoare X şi o valoare x a acesteia (care reprezintă un număr real fixat!).
Valorile efectiv obţinute prin procesul de selecţie constituie o mulţime de numere care se spune că reprezintă o realizare a selecţiei bernouliene. Deoarece fiecare din aceste numere pot să reprezinte oricare valoare a lui
n
X , se consideră că valorile de selecţie (numite şi variabile de selecţie) sunt independente şi identic repartizate ca şi X .
Când pentru o variabilă aleatoare X se cunoaşte sau spunem că se cunoaşte repartiţia de probabilitate sau repartiţia statistică a lui
F fX . Există multe tipuri de repartiţii de
probabilitate. Cele mai importante vor fi introduse în acest capitol, dar şi în capitolele următoare.
Repartiţia uniformă
Una dintre repartiţiile de probabilitate importante, dar care este naturală, este repartiţia uniformă pe un interval [ care are densitatea de repartiţie de forma ]ba,
( ) [ ]⎩⎨⎧ ∈
=restîn
, xdacã0
bakxg ,
abk
−=
1 (2.2.)
Variabila având densitatea de repartiţie (2.2) se spune că este repartizată uniform pe . Deci toate valorile variabilei V sunt egal probabile. Funcţia de repartiţie corespunzând
densităţii (2.2) este
V[ ba, ]
( ) ( ) [ ]∫∞−
⎪⎪⎩
⎪⎪⎨
⎧
>
∈
<
−−
==x
abaxduugxG
b xdacã
ba, xdacã
a xdacã
1
0
(2.3.)
Să notăm cu U variabila aleatoare uniformă pe [ ]1,0 , pe care o vom numi pe scurt variabilă uniformă 0 1− . Densitatea de repartiţie şi funcţia de repartiţie a lui U sunt respectiv
( ) [ ]⎩⎨⎧ ∈
=restîn
0,1 xdacã01
xf , ( ) [ ]⎪⎩
⎪⎨
⎧
>∈<
=1 xdacã0,1 xdacã0 xdacã
1
0xxF (2.4.)
Valorile de selecţie asupra variabilelor aleatoare uniforme se numesc numere aleatoare. Nu este posibil să se producă cu calculatorul, printr-un algoritm, secvenţe de numere aleatoare care să fie uniform repartizate pe intervalul [ ]1,0 şi care să fie independente stochastic. De aceea numerele produse cu calculatorul vor fi numite numere pseudoaleatoare şi ele vor putea fi folosite drept numere aleatoare dacă au un comportament cât mai aleator. Un algoritm care produce un număr aleator (pseudoaleator) se numeşte generator de numere aleatoare (pseudoaleatoare). Iterând un generator se poate obţine o selecţie (secvenţă) de numere aleatoare. Când este nevoie de numere aleatoare cu cât mai bune calităţi se pot aplica metode care combină doi sau mai mulţi generatori, rezultând un alt generator care produce secvenţe de numere pseudoaleatoare mai bune. În multe aplicaţii sunt însă suficient de bune numerele produse de generatori simpli, aşa cum se va vedea mai jos.
II.2. Necesitatea simulării numerelor aleatoare
Următoarea teoremă stabileşte legătura între repartiţia uniformă şi repartiţia uniformă pe un interval [ oarecare.
10 −]BA,
Teorema 1. Dacă este o variabilă uniformă U 10 − atunci ( )UabaV −+= este o variabilă uniformă pe [ şi reciproc, dacă V este o variabilă aleatoare uniformă pe ]ba, [ ]ba, atunci variabila
abaVU
−−
=
este uniformă . 10 − Fiind dată o variabilă aleatoare oarecare X , pentru care se cunoaşte funcţia de repartiţie , este adevărată următoarea teoremă care evidenţiază importanţa repartiţiei uniforme F 10 − , adică importanţa numerelor aleatoare. Teorema 2 (Teorema lui Hincin). Variabila aleatoare ( )XF este uniformă iar dacă notăm cu inversa funcţiei atunci
10 −1−F F ( )UF 1− are funcţia de repartiţie (sau cu alte
cuvinte ). F
( ) X=UF −1
Consecinţă practică. Dacă am putea produce valorile de selecţie asupra lui şi am cunoaşte funcţia de repartiţie a lui
nUUU ,...,, 21
U F X , atunci am putea produce valorile de selecţie asupra lui nXXX ,...,, 21 X cu formula ( )ii UX 1−F= , ni ≤≤1 . Dacă funcţia se poate calcula cu un algoritm atunci valorile de selecţie ar putea fi produse (generate) cu calculatorul folosind următorul algoritm:
1−F
iX
Metoda inversă. (Algoritm)
repetă de ori următoarele instrucţiuni ngenerează o valoare de selecţie U uniformă ; 10 −calculează ( )UFX 1−= ;
Spunem că acest algoritm simulează o selecţie de volum asupra lui n X . Dar paşii cei mai importanţi ai algoritmului sunt ultimii doi paşi care produc o valoare de selecţie X , căci dacă putem genera numere U uniforme 10 − şi independente, atunci prin iterarea primului pas am produce valorile selecţiei de volum asupra lui n X .
În metoda inversă, dacă am putea genera valori de selecţie , iU ni ≤≤1 uniforme 10 −
independente, atunci şi ar fi independente. ( )ii UFX 1−=Un algoritm pentru generarea (simularea) unei selecţii , iU ni ≤≤1
10 − de variabile
uniforme va trebui deci să producă în primul rând un număr uniform , iar prin iterări acelaşi algoritm să fie în măsură să producă selecţia de volum ; deci selecţia se va produce printr-o relaţie iterativa de forma
10 −n
( )ii UgU =+1 , iar algoritmul trebuie să permită producerea de numere independente stochastic. Valorile de selecţie se numesc numere aleatoare uniforme sau simplu numere aleatoare.
iU iU
II.3. Metode congruenţiale liniare
Fie X o variabilă aleatoare discretă repartizată uniform având repartiţia
⎟⎟
⎠
⎞
⎜⎜
⎝
⎛
mmm
mX 111
21: …
… (2.5.)
unde . Conform teoremei 1, variabila U uniformă +∈ Nm 10 − se poate obţine astfel
( )1,0∈=mXU (2.5'.)
unde în ultima relaţie împărţirea se execută în real. Valori de selecţie asupra variabilei X se pot obţine prin relaţia congruenţială liniară
( )( )mcaXX ii mod1 +=+ (2.6.) care spunem că defineşte generatorul mixt congruenţial ( )mcaX ,,,0 unde cele patru constante sunt numere naturale date.
Pentru ca generatorul (2.6) să producă numere aleatoare, trebuie alese constantele , , , astfel încât să satisfacă condiţiile:
0Xa c m1. Numerele iX să reprezinte o repartiţie uniformă de forma (2.5); acest lucru se realizează
dacă şirul { }iX are o perioadă λ mare adică dacă cel mai mic număr λ pentru care
λ+= ii X are loc pentru un X λ cât mai mare.Acest lucru se întâmplă dacă cele patru constante se aleg astfel încât relaţia (2.6) să producă toate resturile modulo m , adică numerele 1,..., −m , iar m este mare. 1,0
2. Numerele iX să fie independente stochastic; acest lucru nu este practic posibil deoarece
1+iX depinde de iX conform relaţiei (2.6); este însă posibil ca numerele iX , 1+iX să fie slab dependente stochastic. Această condiţie înseamnă că numerele respective au un coeficient serial de corelaţie ρ apropiat de 0 , ρ definit astfel:
( ) [ ] [ ] [ ][ ] [ ]1
11111 ,
+
+++
−==ρ=ρ
ii
iiiii XVarXVar
XEXEXXEXXCorr (2.7.)
[ ] [ ][ ] [ ] [ ]{ }222iiiii XEXEXEXEXVar −=−=
unde cu se notează valoarea medie a variabilei aleatoare [ ]YE Y .
Din motive impuse de aritmetica sistemelor de calcul, modulul se ia de forma unde
mepm = p este un număr prim. Următoarea teoremă precizează condiţii pentru alegerea
constantei care defineşte generatorul (2.6). aTeorema 3. Perioada maximă a generatorului (2.6) este dacă şi numai dacă ep=λ
( )pa mod1≡ când 2>p( )4mod1≡a când 2=p , pa <<1
(2.8.)
O valoare a lui poate fi , a 1+= kpa ek ≤≤2
10 ,...,, XX. Alegerea lui folosind numai această
teoremă poate însă să producă un şir de întregi care să nu fie aleator; acest şir care conţine sigur numărul , poate să conţină subşiruri de lungime mare care să fie sau crescătoare, sau descrescătoare, sau să prezinte o periodicitate mare.
a
1−λX0
Presupunând că (ceea ce este desigur permis) se observă că şirul { este de forma
00 =X }nX
( ) ( )mb
caXn
n mod1−= , 1−= ab (2.6'.)
Dacă în relaţia precedentă luăm 1+= ba atunci (2.6') devine ( )( )mbCbCncX ss
nnn mod... 12 −+++= (2.6''.) unde este cel mai mic număr natural care satisface relaţia s
( )mbs mod0≡ (2.6'''.) Numărul se numeşte potenţa generatorului mixt congruenţial. Din relaţia (2.6'') rezultă că este de forma unui polinom în modulo şi deci generatorul este cu atât mai bun cu cât potenţa sa este mai mare. S-a constatat că din punct de vedere practic, potenţa trebuie să fie cel puţin 5.
s
s
nXb m
Un caz interesant este acela când unde este apropiat de cuvântul calculatorului pe care se implementează generatorul (2.6). În acest caz se arată că o alegere bună a lui este de forma , . Deci un generator pentru care şi satisfac condiţiile teoremei 3 generează un şir de numere aleatoare întregi de perioada mare şi care nu prezintă regularităţi (nu conţin subşiruri periodice de lungimi mari sau subşiruri cu monotonii periodice). Constanta poate fi deocamdată arbitrar aleasă.
em 2=
...>
ea
a1...22 21 +++= ffa
c
21 > ff mλ
Să vedem cum putem alege constantele generatorului (2.6) astfel încât ρ să fie suficient de mic. Numerele produse de (2.6) au o repartiţie de forma (2.5) în care valorile lui X sunt
. Deci în acest caz avem 1,...,1,0 −m
[ ]m
xXE
m
x∑−
==
1
0 , [ ] [ ]{ }2
1
0
2
XEm
xXVar
m
x −=∑−
=
( )21
0
1
0
2
21
0
1
01
⎟⎠
⎞⎜⎝
⎛−
⎟⎠
⎞⎜⎝
⎛−
=ρ=ρ
∑∑
∑∑−
=
−
=
−
=
−
=
m
x
m
x
m
x
m
x
xxm
xxxsm, ( )( )mcaxxs mod)( += .
Evaluarea lui se face cu dificultate. Se arată că ρ
c±⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠⎞
⎜⎝⎛+⎟
⎠⎞
⎜⎝⎛−≈ρ
2
6611mc
mc
a,
ma
≈c (2.9.)
O valoare conciliantă a lui , care să asigure o valoare mică a lui a ρ este ma ≈ iar din
condiţia se deduce că 0≈ρmc
trebuie să satisfacă ecuaţia
0661 2 =+− xx adică 211324865.0361
21
=−=mc
(2.10.)
Ultima relaţie dă o condiţie pentru alegerea lui . c Notând , se arată că ( kiik XXCorr +=ρ , ) 1ρ≤ρk
,..., , ceea ce asigură o dependenţă
slabă a numerelor depărtate din secvenţa . 1>k
, 21 XX În concluzie, pentru ca generatorul (2.6) să fie bun trebuie să alegem un modul cât mai mare, să selectăm un astfel încât să satisfacă teorema 3 împreună cu o potenţă mai mare decât 5 şi să ia o valoare apropiată de
ma
m , iar să satisfacă (2.10). Numărul de start , care se mai numeşte şi sămânţa generatorului poate fi orice număr natural mai mic decât m .
c 0X
Pentru a produce un număr aleator întreg X , generatorul (2.6) necesită deci efectuarea unei înmulţiri şi a unei adunări şi apoi calcularea restului modulo . Dacă am alege m 0=c atunci am obţine generatorul multiplicativ congruenţial ( )m,aX 0,,0 care are avantajul că necesită numai operaţia de înmulţire el fiind deci de forma
( )( )maXX nn mod1 ≡+ (2.11.) În acest caz trebuie să fie neapărat pozitiv căci altfel şirul produs de (2.11) ar avea toate elementele egale cu zero deci nu ar fi un şir de numere aleatoare uniforme întregi. Condiţia de perioadă maximă se modifică şi ea şi într-un caz interesant din punct de vedere practic se exprimă prin teorema următoare.
0X
Teorema 4. Perioada maximă a generatorului ( )maX ,0,,0 se obţine când: (i) 0X este un întreg pozitiv prim cu m ; (ii) a este rădăcină primitivă modulo m şi dacă ip , ti ≤≤1 sunt numere prime atunci
perioada este ( ) 22 −=λ eee dacă 3≥e
( ) ( )11 −=λ − ppp ee , dacă 2>p( ) ( ) ( )( )te
tete
te ppcmmmcpp λλ=λ ,...,........ 1
11
1
(2.12.)
Alţi generatori de numere uniforme sunt: 1. Generatorul aditiv congruenţial, sau Fibonacci decalat notat ( )mjXXXk k ;;,...,,; 10
şi bazat pe relaţia ( )( )mXXX knjnn mod−− −= (2.13.)
iar o alegere bună a constantelor este 24=j , 55=K , , care asigură o em 2= 31≥e
perioadă . 1255 −≥λ2. Generatorul congruenţial inversiv notat ( ) 1 cu prim0 ,,, −mcaX −m , definit de relaţia
( )( )mcaXX nn mod11 += −− (2.14.)
unde este inversul al lui în raport cu operaţia de înmulţire a claselor de resturi, când acesta există sau este zero altfel. (Dacă
1−iX ( )mmod iX
prim−m şi atunci inversul există!).
0≠iX
3. Generatorul matricial congruenţial care este de forma ( )( )mCA nn mod1 += −XX (2.15.)
unde este un vector nX ldimensionad − iar , sunt matrici . Înmulţirea matriceală va produce vectori cu componente corelate. Acest tip de generatori se utilizează pe calculatoare paralele.
A C dd ×
nX
4. Generatori bazaţi pe registre de deplasare care utilizează reprezentarea binară a numerelor întregi în calculator. Dacă notăm cu ia , pi ≤≤1 cifrele binare ale lui 1−nx şi considerăm cifrele iC nu toate egale cu zero, atunci generarea cifrelor ia ale lui lui nX se realizează prin relaţia
( )( )2mod... 1111 −+−−− +++ ipippi acaca= pi ca (2.16.) În practică se foloseşte forma mai simplă
( )( )2modqpipi aa +−= +ia =⊕
, 0>> qp (2.16'.) sau dacă notăm operaţia binară or-exclusiv ca suma de biţi modulo atunci relaţia precedentă devine
2
qpipii aaa +−− ⊕= (2.16''.) Să observăm că relaţia de recurenţă referitoare la biţii este aceeaşi cu relaţia de recurenţă a numerelor aleatoare interpretate ca
ia
iX tuplurip − de biţi, adică
qpnpnn XXX +−− ⊕= (2.17.) Un generator bun este de exemplu generatorul qnpnn XXX 33 −− ⊕= , , 521=p 32=q ,
care are o perioadă . Formula (2.17) se poate extinde matricial sub forma 15212 −=λqpnpnn +−− ⊕= XXX (2.17'.)
5. Amestecarea de generatori se obţine astfel: Să presupunem că se dau doi generatori 1G ,
2G şi folosim o tabelă [ ]kT ,...,1 (de ex 64−k sau 128=k ). Algoritmul de amestecare a celor doi generatori este
Iniţializăm cu numere produse cu ; [ ] UiT =
2Gi iU 1G
Generăm cu un indice aleator { }kj ,...,2,1∈ ; Luăm ; Generăm U cu şi luĂm [ ]jTY = 1G [ ] UjT =
Generarea indicelui aleator j , care este un întreg uniform cu valori în { se face astfel
}k,...,2,1
Generează U cu 2G( ) 1+=j × kUtruncIa .
Generatorul rezultat din amestecare este mai aleator decât părinţii săi , , iar perioadele satisfac relaţia
G 1G 2G( ) ( ) ( )( )21 ,..... GGcmmmcG λλ=λ .
CAPITOLUL III. SIMULAREA VARIABILELOR ALEATOARE NEUNIFORME
III.1. Metoda inversă Această metodă a fost deja introdusă ca o consecinţă directă a teoremei lui Hincin. Ea se aplică în cazul în care funcţia de repartiţie se poate inversa uşor. În tabela următoare se prezintă câteva repartiţii de probabilitate ce se pot simula cu metoda inversă.
În formulele lui din tabel se precizează numai valorile lui f x pentru care , presupunându-se că în rest.
0)( >xf0(xf ) =
În ultima coloană se scrie direct formula cu care se simulează variabila aleatoare X , deoarece s-a înlocuit cu u , economisindu-se în acest fel o operaţie de scăadere. Înlocuirea lui cu u se va face ori de câte ori se va ivi ocazia.
u−1u−1Folosind metoda inversă se poate construi uşor un algoritm pentru simularea variabilei
discrete
⎟⎟⎠
⎞⎜⎜⎝
⎛
m
m
pppaaa
X……
21
21: , 11
=∑=
m
iip (3.1.)
Pentru simularea variabilei aleatoare X calculăm mai întâi valorile funcţiei de repartiţie a acesteia. Avem
⎪⎪⎪⎪
⎩
⎪⎪⎪⎪
⎨
⎧
≤
<≤
<≤<≤
<
+++
+=
+
x
ax
axax
ppp
ppp
xF
kk
m
1k
32
21
1
2`
21
1
a dacã
a dacã
a dacãa dacã
a xdacã
1
...
0
)(
…
…
…
(3.2.)
Algoritmul pentru simularea lui X este Algoritmul SIMDISCRV
Intrare tabela [ ]mT ..1 , . Intrare [ ] ∑=
=i
jjpiT
1
[ ]ma ..1
Generează U uniform 10 − cu randomU =: ; 1:=i ; while [ ]iFU > do 1: += ii ;
Ia . [ ]iaX =III.2. Metoda compunerii sau amestecării
Această metodă se aplică variabilelor aleatoare X a căror repartiţie de probabilitate satisface următoarea definiţie.
Definiţie 1. Funcţia de repartiţie este o amestecare (sau compunere sau mixtură) discretă a mulţimii de funcţii de repartiţie
)(xF{ } mi≤≤1i xF )( cu repartiţia disretă
⎟⎟⎠
⎞⎜⎜⎝
⎛
mpppm
J……
21
21: , 1
1=∑
=
m
iip (3.3.)
dacă
∑=
=m
iii xFpxF
1)()( (3.4.)
Funcţia de repartiţie este o amestecare continuă a familiei de funcţii de repartiţie , cu funcţia de repartiţie continuă
)(xF( ){ } RYYxG ∈, ( )yH a lui Y dacă ea este de forma
( )∫=R
ydHYxGxF )(,)( (3.5.)
ultima integrală fiind integrală Stieltjes. Dacă notăm cu X variabila aleatoare care are funcţia de repartiţie şi cu variabila aleatoare care are funcţia de repartiţie atunci amestecarea discretă poate fi interpretată astfel:
)(xF iX)(xFi
iXX = cu probabilitata , ipde unde rezultă următorul algoritm pentru simularea lui X : Algoritmul COMPDISCR
Generează un indice j având repartiţia (3.3); Generează având funcţia de repartiţie ; jx )(xFj
Ia . jXX =:O interpretare analogă are şi amestecarea continuă, de unde se deduce următorul algoritm de simulare a variabilei aleatoare )(~ xFX → Algoritmul COMPCONT
Generează variabila Y care are funcţia de repartiţie ; ( )yHGenerează variabila aleatoare care are funcţia de
repartiţie ; YZ
( )YxG ,Ia . YZX =:
Desigur, în algoritmii precedenţi se presupun cunoscute metode de generare a variabilelor
aleatoare , , J iX Y , . Dacă în definiţia compunerii considerăm în loc de funcţii de repartiţie, densităţi de repartiţie atunci formulele (3.4) şi (3.5) se scriu respectiv:
YZ
∑=
=m
iii xfpxf
1)()( (3.4'.)
∫=R
dyyhyxgxf )(),()( (3.5'.)
III.3. Metoda respingerii
Această metodă (denumită în mod pesimistic ca metodă a respingerii), ar putea fi denumită şi metoda acceptării-respingerii. Ea are forma generală următoare dacă ne referim la simularea unei variabile aleatoare X .
Presupunem că se cunosc următoarele elemente: - Se cunoaşte un procedeu de simulare a unei variabile aleatoare N care ia valori naturale
pozitive; - Se cunosc metode pemtru simularea unor variabile aleatoare SSi ∈ , 1≥i , unde S este o
familie de variabile aleatoare dată; - Se cunoaşte un predicat )nS care se poate calcula iS( SS ,...,, 21P ∀ , ni ≤≤1 ; (Acest
predicat sau condiţie trebuie să poată fi evaluat fără calcule de mare complexitate!); - Se cunoaşte funcţia ψ astfel încât { }( )nSSSX ,...,, 21ψ= , ( ) trueSSS n =,...,, 21P .
Când pentru o variabilă aleatoare X este posibil să se construiască elementele precizate anterior, atunci se poate construi un algoritm pentru simularea lui X sub forma generală următoare: ALGRES=Algoritm General de Respingere repeat
Simulează o valoare n a lui ; NSimulează valorile de selecţie din S; nSSS ,...,, 21
until ( ) trueSSS n =,...,, 21P ;
Ia . ( )nSSSX ,...,, 21ψ=Să observăm că dacă ( ) falseSSS n =,...,, 21P atunci mulţimea de variabile aleatoare
se respinge, de unde provine şi numele de metodă de respingere. În aceeaşi ordine de idei, se observă că dacă
{ nSSS ,...,, 21 }( )( )trueSnSSPpa == ,...,, 21P , numită probabilitate de
acceptare, este apropiată de 1, atunci algoritmul este bun; în caz contrar algoritmul este lent. Teorema 1 (înfăşurătoarei). Fie dată o variabilă aleatoare X care are densitatea de
repartiţie , , pe care dorim să o simulăm. Fie )(xf Rx∈ Y o altă variabilă aleatoare a cărei densitate de repartiţie este astfel încât densităţile , au acelaşi suport (adică iau valori diferite de zero pe aceeaşi mulţime ). Presupunem că există o constantă
)(xh f h SRS ⊂ α ,
, astfel încât , ∞<α<0 )(xhf α)(x ≤ Sx∈∀ . În aceste condiţii dacă este o variabilă uniformă independentă de
U10 − Y , atunci densitatea de repartiţie a variabilei Y , condicţionată
de ))(() Yhα(Yf≤0 este . U≤ f Probabilitatea de acceptare este α= 1ap , de unde rezultă că pentu o metodă a infăşurătoarei ne banală trebuie să avem 1>α . Procedura de respingere este formată din următoarele elemente: (variabila aleatoare constantă); 2=N { }YU ,=S ; dacă ( ) trueYU =,P
))(()(0 YhYfU α≤≤ ; ( ) YYU =ψ , , (adică proiecţia). XTeorema 2. Fie o variabilă aleatoare care are funţia de repartiţie de forma:
∫∞−
=x
QcxF )( ( )φ tdRt )()( (3.6.)
unde este funcţia de repartiţie a unei variabile aleatoare )(zQ Z , [ ]MZ ,0∈ , este o funcţie ce ia valori în (unde
)(xφ[ ]M, M0 poate lua şi valoarea ∞ ) iar este funcţia de
repartiţie a unei variabile aleatoare )(yR
Y , RY ∈ , iar variabilele Z , Y sunt independente stochastic. În aceste condiţii funcţia de repartiţie a variabilei Y , condiţionată de este
. )(Yφ≤Z
)(xF Se observă din nou că probabilitatea de acceptare a algoritmului este . Este evident că teorema descrie un procedeu de respingere ale cărui elemente sunt uşor de precizat.
)(BPpa =
Teorema rămâne valabilă dacă condiţia (3.6.) se scrie în termeni de densităţi şi anume: ( ))() rxcQ φ= )(( xxf , )()( xRxr ′= (3.7.)
O formă duală a teoremei se obţine dacă este de forma )(xF
( )( )( )φ tdRtQ )(∫∞−
−=x
cxF 1)(
)(Yφ
(3.8.)
iar predicatul devine iar condiţia (3.7.) se scrie în termen de densităţi. Z ≥( ( )) )()( xrxQ1)( cxf −=
)(~ 00 xGZ →
φ (3.7'.) Teorema 3 (a şirului descendent). Presupunem date variabilele aleatoare
. şi , independente stochastic. Atunci, următoarele afirmaţii sunt adevărate
)(~ xGZi → 1≥i
(1) Dacă x este fixat şi numărul k este fixat, atunci
( ) [ ] [ ]!)(
)!1()(...2 ZZ k≥≥≥
1
11 kxG
nxGZZxP
kk
k −−
=<≥−
− (3.9.)
(2) Dacă x este fixat şi K este indicele aleator la care se ''rupe'' subşirul descendent (ca la pct (1)), atunci
( ) ( ) )(1)2(mod xGeKPKP −==== imparnumãr Z
(3.9'.) (3) Dacă subşirul descendent este kk ZZZ <≥≥1 ...≥0 −1 (adică se rupe la un K aleator şi
începe cu )(G ), atunci ~ 00Z → x
( ) ∫∞−
−==<x
tGa tdGepKxZP )(1)2(mod 0
)(0 (3.10.)
unde este constanta de normare ap1
0)( )(
−∞
∞−
−⎥⎦
⎤⎢⎣
⎡= ∫ xdGep xG
a (3.10'.)
Algoritmul de respingere rezultat din teoremă este urmărorul Algoritmul RESPSIRD bazat pe şiruri descendente. repeat
1+= KK ; 10 : ZZ = ; Generează ; )(~1 xGZ → end; until ; 12mod =K Ia . *: ZX = Să analizăm acum performanţa algoritmului. Observăm că ne dă informaţie asupra vitezei algoritmului: cu cât este mai mare, cu atât mai repede ajungem să acceptăm un (când ). Dar nu este de ajuns pentru a caracteriza performanţa algoritmului; ar trebui să vedem şi câte numere , sunt necesare pentru a obţine un
. Dacă este numărul total de valori de selecţie , necesare până la
acceptarea unui , atunci avem:
ap
iZ
ap 0Zimparnumãr =K
acceptat *N
0Z
ap
iZ 0≥i
0 −Z 0≥i
( ) ⎟⎟⎠
⎞⎜⎜⎝
⎛+= ∫
∞
∞−
)(110
)(* xdGep
NE xG
a
(3.11.)
III.4. Alte metode Metodele de simulare a variabilelor aleatoare neuniforme X pot fi obţinute fie prin transformarea unor variabile aleatoare uniforme 10 − , fie sub forma
( )nYYYTX ,...,, 21= (3.12.) unde variabilele pot fi simulate. De fapt, metodele prezentate în primele subsecţiuni ale acestui capitol sunt toate de acelaşi tip: simularea lui
iYX se reduce la simularea unor variabile mai simple
, , unde şi n poate fi aleator. iY ni ≤≤1 În tabelul următor se dau câteva repartiţii de probabilitate a căror simulare se realizează prin transformări simple de variabile U uniforme 10 − dar nu prin metoda inversă. În ultima coloană se precizează transformarea T .
Nr. crt.
Numele repartiţiei
Densitatea de repartiţie Transformarea T
1. Normală
( )1,0(N2
2
21)(
x
exf−
π= , ∞<<∞− x ∑
=
=12
1iiUX
2. Modul [ ]
⎩⎨⎧ ∈−
=altfel
1,1- xdacã01
)(x
xf 21 UUX −=
3. Repartiţia maximului
[ ]⎩⎨⎧ ∈
=−
altfel0,1 xdacã
0)(
1nnxxf { }UUU ,...,,max n21
4. Repartiţia minimului
( ) [ ]⎩⎨⎧ ∈−
=−
altfel0,1 xdacã
01)(
1nxnxf { }UUU ,...,,min n21
5. Repartiţia
)(kErlang ⎪⎩
⎪⎨
⎧
≥
<
Γ=
0 xdacã
0 xdacã
(k)1
0)(xf , +∈ Nk ⎟⎟
⎠
⎞⎜⎜⎝
⎛−= ∏
=
k
iiUX
1
ln
Definiţia 2. Dacă , iZ ki ≤≤1
)(kErlang
sunt variabile independente, atunci variabila
este o variabilă .
)1(Exp
∑=
=n
iiZX
1
Când o familie de variabile aleatoare are proprietatea că suma a două variabile independente este tot o variabilă din aceeaşi familie spunem că aceasta este o familie stabilă de variabile aleatoare. Pentru stabilitatea unor familii de variabile aleatoare putem da următoarele enunţuri. Teorema 4. Dacă ),1,0(~ υ→ GammaX , ),1,0(~ μ→ GammaY şi X şi Y sunt independente stochastic, atunci ),1,0(~ μ+υ→+= GammaYXZ (familia variabilelor aleatoare Gamma standard este stabilă). Teorema 5. Dacă ( )σ→ ,~ mNX şi , atunci ( λ→ ,~ pNY )
( )22,~ λ+σ+→+= pmNYXZ (familia variabilelor normale este stabilă).
III.5. Simularea repartiţiilor înrudite cu repartiţia normală Repartiţia 2χ
Dacă , sunt variabile normale independente atunci 1Z fi ≤≤1 )1,0(N
∑=
=χj
iiZ
1
22 (3.13.)
se numeşte variabilă hi pătrat centrată cu grade de libertate notată . f 2fχ
Dacă , atunci variabila dată prin relaţia (3.11.) se numeşte variabilă hi
pătrat necentrată cu grade de libertate şi cu parametrul de excentricitate , notată , unde
( 1,~ ii mNZ →
f)
δ 2,δχ f
∑=
δ=δf
iiim
1
222 (3.13'.)
Observaţie. Variabila centrată este o variabilă de tip 2fχ ⎟
⎠⎞
⎜⎝⎛
2,
21,0 fGamma .
Formula (3.24) permite simularea directă, pornind de la definiţie a variabilelor , . 2fχ 2
,δχ f
Repartiţia Student
Dacă este independentă de variabila hi pătrat atunci variabila )1,0(~ NZ → 2
fχ
f
Ztf
f 2χ=
(3.14.)
se numeşte variabila a lui Student cu grade de libertate. Dacă în (3.14) în loc de se
introduce atunci se obţine variabila notată care se numeşte variabila t Student
necentrată cu grade de libertate şi cu parametrul de excentricitate
t f 2fχ
2,δχ f δ,ft
f δ . Variabilele de tip t Student se simulează direct cu formule de tipul (3.14).
Repartiţia Snedecor Dacă variabilele , sunt independente atunci variabila 2
1fχ 2
2fχ
21
12
21f
fff f
fF
χ
χ= (3.15.)
se numeşte variabila a lui Snedecor centrată cu , grade de libertate (în notaţia indicilor contează ordinea!). Dacă în (3.15) se foloseşte câte una din
F 1f 2f
1,1 δχ f
0, ,1fF
, sau ambele,
atunci se obţin variabilele Snedecor simplu necentrate , , cu parametri
corespunzători de excentricitate, sau variabila Snedecor dublu necentrată .
Variabilele pot fi de asemenea simulate direct prin formule de tipul (3.15.).
2,2 δχ f
2,0;2 δfF1;2 δf,1f
F
F2,1;2,1 δδffF
F
Repartiţia lognormală Variabila aleatoare pozitivă Y se numeşte lognormală ( )σμ,LN de parametri μ , σ dacă variabila este normală ( )YX log= ( )συ,N .
Dacă se dau media şi dispersia pentru variabila lognormală m 2s Y atunci media μ şi dispersia pentru 2σ X se calculează cu formulele:
( ) ⎟⎟⎠
⎞⎜⎜⎝
⎛+−=μ 1log
21log 2
2
msm
⎟⎟⎠
⎞⎜⎜⎝
⎛+=σ 1log 2
22
ms
(3.16.)
deci simularea variabilei Y se realizează cu algoritmul simplu Algoritmul LNORM pentru simularea lognormalei Intrare , ; calculează m 2s μ, σ ;
(acesta este un pas pregatitor!); Generează ( )1,0~ NZ → ;
Calculează σ+μ= ZX ; (unde 2σ=σ );
Ia . (XeY = ( )σμ→ ,~ LNY )
Familia de variabile de tip Johnson Din aceasta familie fac parte o serie de variabile ce se obţin prin transformări de variabile normale . Aceste transformări sunt: ( σ→ ,~ mNX )- [ ] ξ++λ= xeY 1 = variabila logit normală; - ξ+λ = variabila lognormală; = XeY
- ( ) ξ+ = variabila 1sinh − normală, λ= XY sinh2
)sinh(xx eex
−−= , unde 0>λ , 0≥ξ ,
aceste constante alegându-se de obicei astfel încât [ ] 0=YE , [ ] 1=YVar .
III..6. Simularea unor variabile particulare
III.6.1. Repartiţia exponenţială Este suficient să ne ocupăm de repartiţia variabilei a cărei densitate de repartiţie este
)1(~ ExpZ →
⎩⎨⎧
>≤
=0 xdacã0 xdacã
e0
)( x-xf , (3.17.)
căci simularea variabilei ( )λ→ ExpX ~ se realizează cu λ
=ZX . Metoda inversă pentru
simularea lui Z , adică , nu este recomandabilă deoarece dacă este apropiat de zero atunci nu se poate calcula. De aceea vom prezenta o metodă de respingere datorată lui John von Newmann care se bazează pe teorema următoare.
( )Ulog−=Z U(log )U
Teorema 6. În teorema 3. (a şirului descendent) se consideră 00 UZ = , , unde , sunt uniforme . Dacă notăm cu numărul aleator de subşiruri descendente respinse până când se acceptă un subşir, atunci
ii UZ = 1≥i
0U iU 10 − N( )1Exp~0ZN →X += , unde este cel
acceptat (din ultimul subşir descendent). 0Z
De aici se deduce următorul algoritm pentru simularea lui Z . Algoritmul EXRJ de simulare a exponenţialei prin metoda respingerii Iniţializează 0:=N ; repeat generează , uniforme 0U 1U 10 − şi independente;
if 02mod =K then 1: += NN ; (se numără subşirurile respinse);
until 12mod =K ;
Ia . *: UNZ +=Pe lângă probabilitatea de acceptare , performanţa algoritmului este
caracterizată şi de numărul al variabilelor uniforme
11 −−= epa*N { } 0≥iiU necesare a fi generate până când
se obţine o valoare de selecţie exponenţială Z . Conform (3.11.) avem
( ) ( ) 8,31
111
111 2
1
1
0
* ≈−
=−+−
=⎟⎟⎠
⎞⎜⎜⎝
⎛+= −∫ e
eee
dzep
NE x
a
(3.11'.)
adică trebuie simulate în medie aproape 4 numere uniforme pentru a obţine o valoare de selecţie exponenţială Z .
III.6.2. Repartiţia Gamma
Repartiţia are densitatea dată de formula ( υλα ,,Gamma )
( ) ( ) ( )⎪⎩
⎪⎨
⎧
≥
<
α−υΓλ= α−λ−−υυ
0 xdacã
0 xdacã0)( 1 xex
xf
iar formulele
λ+α=
XY , ( )λα−= YX
spun că este necesar şi esenţial să construim metode pentru simularea unei repartiţii Gamma standard, adică , . Dacă , atunci repartiţia gamma devine repartiţie . Dacă
( )υ,1,0Gamma)(k 0
+∈υ R1
+∈=υ NkErlang <υ< atunci simularea variabilei se realizează
cu o metodă de respingere bazată pe înfăşurarea cu o densitate ( υ,1,0Gamma )
( )υ,1,0Weibull cu densitatea de repartiţie dată de
( )⎪⎩
⎪⎨⎧
≥<
υ= υ−υ 0 xdacã
0 xdacãx
01- xe
xh
Algoritmul pentru simularea lui X este: Algoritmul GAMRESM1 (algoritm de respingere pentru variabila gamma standard cu parametru
subunitar)
Intrare ; Calculează υυ
=1:c ; υ−
υ
υ=ζ 1: ; ; ( )1: −υζ= ea
repeat Generează U uniform 10 − ;
Calculează ( )UZ ln−= ; cZY = ;
{se generează ( )υ→ ,1,0~ WeibullY } Generează alt U uniform 10 − ;
until ; YZaeU −≤Ia . YX =:
O metodă de compunere-respingere pentru cazul 10 <υ< . Vom scrie densitatea de repartiţie ( )υ,1,0Gamma , 10 <υ< dată de expresia
( )⎪⎩
⎪⎨
⎧
≥
<
υΓ= −−υ 0 xdacã
0 xdacã1
0)( 1 xex
xf
sub forma )()()( 2211 xfpxfpxf +=
unde
[ ]⎪⎩
⎪⎨
⎧∈
=altfel
0,1 xdacã
0
)()( 11 p
xfxf ,
( )⎪⎩
⎪⎨
⎧∞∈
=altfel
1, xdacã
0
)()( 22 p
xfxf
( )( )υΓυΓ
=:1
1p , ( ) ( )
( )υΓυΓ−υΓ
=−=;11 12 pp
( ) ∫∞
−−υ=υΓ0
1 dxex x , ( ) ∫ −−υ=υΓ1
0
1;1 dxex x
Să notăm , )(~ 11 xfX → ( )xfX 22 ~→ . Atunci are loc următoarea teoremă: Teorema 7. Variabila se simulează folosind teorema 3. a subşirului descendent unde 1X
υ= 100 UZ , cu { } uniforme ii UZ = 0≥iiU 10 − , iar se simulează cu ajutorul teoremei 3.2.
duale unde densitatea este de forma 2X
)(2 xf( ) )()(1)(2 xrxQcxf == , , 0>x
( ) ( )( )vec
:11Γ−υΓ
= ,
⎩⎨⎧
≥<
=+ 1 xdacã
1 xdacãe0
)( 1x-xr ,
⎩⎨⎧
≥<
=υ 1 xdacã
1 xdacãx-1
0)( 1-xQ .
Presupunând că s-au calculat în prealabil ( )υΓ=g ; ( )υΓ= ;11g ; gg
p 11 = ,
ggg
p 12
−= ,
υ=
1a , v
b−
−=1
1, ( )1
1gge
c−
= , algoritmul de simulare al variabilei
, este ( )υ→ ,1,0~ GammaX 10 <υ<Algoritmul GAMCOMNL1 pentru simularea variabilei gamma de parametru subunitar prin compunere şi respingere:
Intrare , ; 1p 2pGenerează ; 1-0 uniform~→Uif then begin 1pU ≤
Generează ; Ia )(~ 11 xfX → 1XX = ; end
else begin Generează ; Ia )(~ 22 xfX → 2XX =
end. Variabilele şi se generează cu algoritmii de respingere prezentaţi în continuare. 1X 2X
Algoritmul GAMRESNL1 pentru simularea variabilei . 1Xrepeat
0≥i Z , Y necesare pentru a genera în final un X este ( ) ( ) ( ) ( ) 2*
11*22
*11
* +=+= NEpNEpNEpNE .
III.6.3. Metode pentru simularea variabilei ( )υ,1,0Gamma , 1>υ Vom prezenta doi algoritmi de respingere bazaţi pe metoda înfăşurătoarei.
G1 Să considerăm ( )υ→→ ,1,0~)(~ GammaxfX , 1>υ şi să luăm ca înfăşurătoare
densitatea ⎟⎠⎞
⎜⎝⎛υ
→1~)( Expxh , adică
⎪⎩
⎪⎨
⎧
≥
<
υ
=υ
− 0 xdacã
0 xdacã
1
0)( x
exh
Analizând raportul
( ) υ−
−−υ
υΓ
υ== x
x
e
exxhxfxr
1
)()()( , 1>υ , se constată că acesta are ca punct de
maxim υ=x , deci constanta din teorema 3.1. (a înfăşurătoarei) este α ( ) ( )υΓυ
=υ=αυ−υ 1er .
Probabilitatea de acceptare a algoritmului de simulare (care este simplu de construit) este deci
( )1
22
1υυΓ
υ−υe( )
−υπ
≈=epa
)
. S-a folosit aproximarea lui Stirling pentru adică ∞→υ
( ) ( ) ( 1−211 υπ≈υΓ −υ−−υ e1−υ , deci această probabilitate scade de ordinul υ .
G2 Algoritmul precedent este lent pentru υ foarte mare. De aceea vom prezenta în acest caz o altă metodă de respingere, bazată pe înfăşurarea cu o densitate Cauchy nonstandard trunchiată pe [ de forma )∞,0
( )( )c
xkxh 211
)(−υ−
+= , 0≥x
(3.18)
unde k este o constantă de normare. Putem enunţa teorema: Teorema 8. Dacă înfăşurăm densitatea ( )υ,1,0Gamma , 1>υ , cu densitatea dată
de (3.18) atunci pentru avem )(xh
12 −υ≥c
( ) ( ) ( )1111)()()( −υ−−υ−υ
υΓ=α≤= e
kxhxfxr (3.18'.)
Pentru a descrie algoritmul dedus din teorema precedentă presupunem calculate în prealabil constantele , 1−υ=b bc +υ= , 12 −υ=S . Atunci algoritmul pentru generarea variabilei , ( ) 1>υ,1,0Gamma→~X υ este Algoritmul GAMCAUCHY (simularea variabilei ( )υ,1,0Gamma , 1>υ prin înfăşurarea cu o densitate Cauchy)
repeat repeat
Genereaz random:=U şi ia ( )( )5.0: −π⋅= UtgsT ;
(T este Cauchy standard pe ( )+∞∞− , ; Ia TbY += ; (Y este Cauchy non standard pe ( )+∞∞− , ;
until ; (Se aplica respingerea pentru a obţine ; 0>Y trunchiatã−YGenerează random:1 =U ;
until ⎟⎟
⎠
⎞
⎜⎜
⎝
⎛++−⎟
⎠⎞
⎜⎝⎛
≤c
TTbYb
eU
21loglog
1 ; Ia YX = .
Se observă că constanta nu intervine în construcţia algoritmului dar ea este necesară
pentru a calcula
k
α=
1ap . Un calcul simplu arată că
1
121
2
−
⎥⎦
⎤⎢⎣
⎡⎟⎟⎠
⎞⎜⎜⎝
⎛−υ
−υ−+
π= arctgk . Referitor
la repartiţia , mai observăm şi faptul că dacă descompunem sub forma
, ,
( )υ,1,0Gamma[ ] +∈υ= Nk
1>υ υpk +=υ [ )1,0∈−υ=p k şi considerăm variabilele ( )υ→ ,1,0~X Gamma ,
, )(kErlang~Ek → ( )p,1,0GammaY ~→ , atunci simularea lui X se realizează cu relaţia , , Y kEEX k += Y independente.
III.6.4. Repartiţia Beta
Variabila X are repartiţia , , dacă densitatea sa de repartiţie este ),( baBeta 0>a 0>b
( ) [ ]⎪⎩
⎪⎨
⎧∈−
=−−
altfel
1,0 xdacã
0
1),(
1)(
11 ba xxbaBxf , (3.19.)
unde
∫ −− −=1
0
11 )1(),( dxxxbaB ba , )()()(),(
bababaB
+ΓΓΓ
= . (3.19'.)
Următoarea teoremă permite simularea variabilei . ),( baBeta Teorema 9. Dacă ( )aGammaX ,1,0~1 → , , independent de , atunci variabila
),1,0(~2 bGammaX → 1X
2X
21
1
XXXX+
= (3.20.)
este o variabilă . ),( baBetta Simularea variabilei poate fi deci făcută cu (3.20). Dar deoarece această formulă presupune generarea prealabilă a două variabile gamma, rezultă că această metodă are o complexitate mare. De aceea, în cazuri particulare se pot folosi următoarele teoreme:
),( baBeta
Teorema 10. Fie şi +∈ Nba, 1−+= ban şi să considerăm variabilele
uniforme . Să construim statisticile de ordine nUUU ,...,, 21
10 − )()2()1( ... nUUU <<< prin ordonarea
valorilor { } . Atunci n≤iiU ≤1 ( )ba,BU a ~)( → .
Teorema 11. Dacă 10 << a , 10 << b şi , sunt variabile aleatoare uniforme
şi independente, şi dacă
1U 2U
10 − aUV =1
1 , b1
UT 2= , atunci repartiţia variabilei TV
VX+
=
condiţionată de este . 1<+TV ),b(aBeatTeorema 12. Dacă 10 << a , şi , sunt variabile uniforme 1>b 1U 2U 10 −
independente şi considerăm aUV , 1
1= 11
2 −= bUT , atunci repartiţia variabilei V condiţionată de este . 1<+TV ),( baBeta
Algoritmii ce rezultă din aceste ultime teoreme sunt uşor de construit. Vom prezenta algoritmul ce rezultă din teorema 11. Algoritmul BETA14 pentru simularea variabilei , ),( baBeta 1,0 << ba
repeat Generează , uniforme 1U 2U 10 − şi independente;
Ia aUV1
1= , bUT1
2= ;
until ; 1<+TV
Calculează TV
VX+
= .
Probabilitatea de acceptare a acestui algoritm de respingere este ),( baBba
abpa += .
Algoritmul de respingere construit pe baza teoremei 3.11 are probabilitatea de acceptare . ),( baaBpa =
III.6.5. Repartiţia normală
Ne vom opri la simularea variabilei ( )1,0~ NZ → căci dacă ştim să simulăm această
variabilă atunci se simulează cu formula ( σ→ ,~ mNX ) σ+= ZmX .
O metodă de compunere-respingere. Să considerăm variabila normală standard ( )1,0N , notată cu , care are densitatea 1X
⎪⎩
⎪⎨
⎧
≥<
π
= − 0 xdacã0 xdacã
2
0)(
2
21
x
exf (3.21.)
Considerăm şi variabila 12 XX −= cu densitatea . Densitatea se
scrie
)(2 xf )1,0(~ NX →
)(21)(
21)( 21 xfxfxf += , adică este o compunere discretă a densităţilor şi .
Pentru a construi un algoritm de simulare a lui
)(1 xf (2 xf )
X va trebui să construim mai întâi un algoritm de simulare a lui . Vom înfăşura densitatea cu . Rezultă teorema: 1X )(x1f )1()(xh ~ Exp→
Teorema 13. Dacă înfăşurăm cu avem )(1 xf )(xhπ
=α≤e
xhxf 2)()(
şi deci putem
aplica teorema de respingere a înfăşurătoarei pentru simularea variabilei . 1X Algoritmul pemtru simularea lui ( )1,0~ NX → este: Algoritmul RJNORM de compunere-respingere pentru simularea normalei )1,0(N
Algoritmul corespunzător metodei polare se deduce cu uşurinţă şi el produce simultam două valori de selecţie , şi , independente. )1,0(N 1Z 2Z
III.7. Simularea unor variabile discrete Dacă X este o variabilă discretă care ia ca valori şirul { } ∞≤≤nna 1 şi se cunoaşte funcţia de frecvenţă , iaXPi ∈= ()f ( −f calculabilă, atunci folosind proprietatea 0)(lim =
∞→if
i,
(dedusă din ) se poate construi uşor un algoritm bazat pe metoda respingerii (măcar
aproximativ!).
1)(1
=∑∞
=iif
În această secţiune vom prezenta algoritmi pentru simularea unor repartiţii particulare.
Simularea unor repartiţii bazate pe probe Bernoulli
Fie un eveniment aleator observabil care are probabilitatea constantă . Într-o experienţă întâmplătoare se poate produce cu probabilitatea
A 0)( >= APpA p sau evenimentul contrar
A cu probabilitatea . O astfel de experienţă se numeşte probă Bernoulli. Când se produce spunem că avem de-a face cu un succes, iar când nu se produce spunem că se realizează un eşec. Să asociem unei probe Bernoulli variabila aleatoare
p−1q =A A
Z astfel încât 1=Z dacă se produce şi dacă se produce A 0=Z A , adică Z are repartiţia
⎜⎜⎝
⎛pq
Z10
: ⎟⎟⎠
⎞, , ( ) pZE = ( ) ( )pppqZVar −== 1 . (3.22.)
Funcţia de repartiţie a lui Z este
⎪⎩
⎪⎨
⎧
≥<≤
<=<=
1 xdacã1x0 dacã
0 xdacã
1q0
)( xZP)(xF . (3.22'.)
De aici rezultă că algoritmul de simulare a lui Z prin metoda inversă este Algoritmul BERN de simulare a unei variabile (probe) Bernoulli
Generează ; random:=Uif then 0:=Z else 1:=Z . pU >
Să observăm din nou că dacă 21
=p , atunci suntem în cazul particular al aruncării la
întâmplare cu banul.
Repartiţia binomială Se spune că variabila aleatoare discretă NX ∈ este o variabilă binomială , , dacă
),( pnBinom+∈ Nn 1<< p0 X este numărul de succese în n probe Bernoulli independente, adică
, unde sunt variabile identic repartizate Bernoulli, independente. Simularea
variabilei
∑=
=n
iX
1iZ iZ
X se face deci simplu, prin numărarea de succese în probe Bernoulli independente. Avem ,
n( ) α=α= n pC α−α nqXP pq −= 1 , adică ( )α=XP este termenul general al dezvoltării
binomului , de unde derivă şi denumirea de repartiţie binomială. ( )nqp +Media şi dispersia lui X sunt npXE =)( şi npqXVar =)( , iar din teorema limită
centrală se deduce că pentru suficient de mare (n ∞→n ) variabila )1,0(~ Nnpq
npXWn
~ BinomX →
→−
= .
De aici rezultă următorul algoritm simplu de generare a lui pentru mare. ),( pn n
Algoritmul BINNORM Generează ; )1,0(~ NW →
Calculează { }npqWnpX +=: . (Notaţia { }E înseamnă ''cel mai apropiat întreg de E .
Repartiţia Pascal Variabila X are repartiţia , , ),( pkPascal +∈ Nk 10 << p , dacă X este numărul de eşecuri până la apariţia a succese într-un şir oarecare de probe Bernoulli independente. De aici rezultă că variabila se simulează cu următorul algoritm care numără eşecurile realizate până la realizarea a succese într-un şir de probe Bernoulli independente.
k~→ ),( pk
kPascalX
Algoritmul PASCAL Intrare , +∈ Nk p , 10 << p ; 0=X ; 0:=j ; ( X numără
eşecurile şi j succesele); repeat
Generează random:=U ; if then pU < 1: += jj else 1: += XX ;
until . (kj = X este valoarea de selecţie generată). Să observăm că , α−
−+α=α= qpCXP kkk1
1)( ,...2,1,0=α , care este termenul general al
dezvoltării în serie a expresiei ( ) kk qp −−1 din care cauză repartiţia se mai numeşte şi repartiţia binomială cu exponent negativ. Avem şi următorul rezultat.
),( pkPascal
Teorema 15. Dacă ( )pkPascalX ,~ 11 → şi ( )pkPascalX ,~ 22 → , sunt variabile independente, atunci ( )p,2kXX kPascalX ~ 121 +→+= (repartiţia Pascal este stabilă). Avem şi că:
pkqXE =)( , 2)(
pkqXVar = (3.23.)
formule ce se folosesc la validarea algoritmului.
Repartiţia geometrică )( pGeom Este un caz particular de repartiţie Pascal, când 1=k . Simularea variabilei
se poate realiza cu algoritmul PASCAL sau cu metoda inversă cu formula )(~ pGeomX →
⎥⎦
⎤⎢⎣
⎡=
)log()log(
qUX de unde se deduce şi
pqXE =)( , 2)(
pqX =Var , formule care se folosesc la
testarea algoritmului de simulare.
Repartiţia hipergeometrică Această repartiţie se introduce după cum urmează: Să considerăm experimentul cu urnă în care bile se extrag la întâmplare din urnă, fără întoarcere. În acest caz numărul
nX de bile albe extrase este o variabilă hipergeometrică. Să
notăm cu U evenimentul care reprezintă extragerea unei bile albe şi cu V evenimentul care constă în extragerea unei bile negre; atunci probabilităţile de a extrage în prima extragere o bilă
albă respectiv neagră, sunt respectiv NAUPp == )( ,
NBVP =)( . Probabilităţile de extragere a
unei bile albe sau negre în a doua extragere sunt condiţionate de rezultatele primei extrageri adică
( )11
−−
=NAUUP , ( )
1−=
NAVUP , ( )
1−=
NBUVP , ( )
11−−
=NBVVP . Se observă că la
fiecare extragere compoziţia urnei se schimbă şi probabilitatea de a extrage o bilă albă sau neagră este variabilă în funcţie de extragerile anterioare. Variabila hipergeometrică se notează
, , de unde ),,( npNH 10 << p Nn < { }NpA = ({ }x , Rx∈ este cel mai apropiat întreg de x ), . Probabilitatea ca în extracţii succesive fără întoarcere, să se extragă bile
albe este
A−NB = n a
nN
anB
aA
CCC −
aXP == )( , na ≤≤0 , Nn < . De asemenea, avem npXE =)( ,
( ) ( ) N+ 11 ( − )[ ]pNpnnp−
−1NXE =2 şi
1)1()(
−−
−NNp=XVar nnp .
X este Algoritmul de simulare a variabilei hipergeometrice Algoritmul HIPERGEOM
Intrare , A B , , n BAN += ; (Acesta este un pas pregăritor);
NAp =Calculează ; IniţializeazAă 0:=j , 0:=X ;
repeat 1: += jjrandom:=U ; Ia ; Generează
if U then begin p<:= XX 1 1:=S
0:+ ; ; (S-a estras o bilă albă);
end else ; (S-a extras o bilă neagră); =S
: , NAp =: ; 1−= N AN SA ==:, Calculează
( )npN ,,HXuntil nj = ~→; ( ).
Repartiţia Poisson
are repartiţia , dacă ( )λPoisson Variabila aleatoare X , NX ∈ 0>λ
( ) λ−α
αλ
=α=XP ( ) +λ= 22XEλ=)(XE λ =)(XVare!
, . Avem şi 0> , . λλ ,
Repartiţia este repartiţia evenimentelor rare în sensul următor: evenimentele sunt independente şi se produc la intervale aleatoare de timp astfel încât un eveniment se produce pe intervalul de timp [ cu probalilitatea
( )λ
tt Δ+,
Poisson
]t ( )tOt Δ+Δλ unde şi ( )ΔtO 0=lim0→Δt
( ) 00
=ΔΔttO
( este neglijabilă în raport cu ( tO Δ ) tΔlim→Δt
) iar probabilitatea ca pe acelaşi interval
să se producă mai mult de un eveniment (condiţia de raritate) este ( )tO Δ (adică este neglijabilă). Numărul de evenimente rare ce se produc pe unitatea de timp este o variabilă aleatoare
. Numărul λ este intensitatea cu care se produc evenimentele rare. ( )λPoisson Intervalul de timp la care se produc două evenimente rare consecutive are o repartiţie
fapt care spune că
θ
( )λExp 1−=X j dacă . Ţinând acum seama de faptul că ∑θj
i∑−
≤θj
i
1
<λ== ii 11
λ−=θ i
iU
log , uniforme , atunci rezultă că iU 10 −
1−= jX dacă ∏∏=
λ−−
=
>≥j
ii
j
ii UeU
1
1
1(3.24.)
Pe baza relaţiei (3.24.) se poate construi cu uşurinţă un algoritm pentru simularea lui X .
III.8. Validarea generatorilor Validare înseamnă nu numai verificarea corectitudinii formale a programelor, ci în special dacă algoritmul implementat (programat) produce valori de selecţie asupra variabilei aleatoare în cauză, adică dacă pe baza unei selecţii , de volum suficient de mare se verifică ipoteza statistică , care se mai numeşte şi ipoteză de concordanţă.
nXXX ,...,, 21 n)(xFH →~: X
Mai întâi ar trebui să verificăm intuitiv dacă repartiţia empirică sau repartiţia de selecţie se aseamănă cu cea teoretică. Mai precis trebuie să construim grafic histograma şi să vedem dacă ea se aseamănă cu densitatea de repartiţie. Construcţia histogramei se face astfel: - se determină mai întâi ( )nXXXm ,...,,1min= 2 , ( )nXXXM ,...,,max 21= care reprezintă
intervalul de variaţie [ ]Mm, ; - se alege un întreg pozitiv k , (practica recomandă 4015 ≤≤ k ) care reprezintă numărul de
intervale ale histogramei; - se împarte intervalul [ ]Mm, în k intervale [ ]iii aaI ,1−= , ki ≤≤1 , ma =0 , Mak = ; - se determină if = numărul valorilor de selecţie ce cad în intervalul [ )ii aa ,1− , ki ≤≤1 ;
se numesc frecvenţe absolute; if
- se determină frecvenţele relative nf
r ii = .
Repartiţia empirică este
⎜⎜⎝
⎛rI
1
1⎟⎟⎠
⎞
k
k
rrII
……
2
2 . (3.25.)
Pentru reprezentarea grafică a valorilor din (3.25.) procedăm astfel: luăm pe abscisă intervalele şi construim dreptunghiuri care au ca bază aceste intervale şi ca înălţimi . Graficul obţinut reprezintă histograma selecţiei date. (Ea este histograma frecvenţelor relative; asemănător se obţine şi histograma frecvenţelor absolute).
iI ir
nXXX ,...,, 21
Testul 2χ
Odată construită histograma putem aplica testul de concordanţă pentru verificarea ipotezei . Pentru aceasta calculăm întâi probabilităţile
2χ)(: xFH →~X ( )11 aFp = ,
, , ( ) )1−i F 2 ≤ i( −ia ( )11 −−= kk aFp . =i aFp 2−≤ k
Calculăm apoi ( )∑
=
k
j
if1
−=
i
i
npnp 2
2χ care se ştie că are repartiţia hi pătrat cu grade
de libertate, ( este variabila corespunzătoare). Fiind dat riscul de genul I, , (o probabilitate
mică, apropiată de zero) determinăm (numită
1−k
21−kχ α
2,1 α−χ k −α cuantilă superioară) astfel încât
( ) α=χ>χ α−−2
,12
1 kkP . Ipoteza H se acceptă dacă şi se respinge în caz contrar. 2,1
2α−χ<χ k
Un test simplu
Cel mai simplu mod de testare a unui generator care simulează o variabilă neuniformă X se poate realiza astfel: - Se determină câteva momente teoretice ale lui X ca de exemplu )(XE şi
)( ; m =
2 XVar=σ- Cu ajutorul selecţiei simulate nX se calculează momentele empirice de selecţie: XX ,...,, 21
n
XX
n
ii∑
== 1 , 21
2
2 Xn
Xs
n
ii
−=∑= .
Se consideră că generatorul este bun dacă pentru suficient de mare ( ) valorile lui
n 1000>nX şi sunt tot mai apropiate de m şi , ca o consecinţă a legii numerelor mari.
Se presupune că momentele teoretice m , sunt cunoscute.
2s 2σ2σ
CAPITOLUL IV. SIMULAREA VECTORILOR ALEATORI
IV.1. Generalităţi
( )′= kXXX ,...,,1 X 2Un vector (vector coloană) ale cărui componente sunt variabile aleatoare, se numeşte vector aleator. Cazul în care componentele sunt independente stochastic este desigur banal, de aceea cazul cel mai interesant este când componentele sunt dependente sau corelate.
iX
iX Funcţia de repartiţie a vectorului X este
( ) ( )kk xXxXxXPxxFxF kx <<<== ,...,,,...,,)( 221121 (4.1.) şi ea se mai numeşte funcţia de repartiţie comună a lui sau funcţia de repartiţie multidimensională. Desigur,
kXXX ,...,, 21
( ) 0,..., =−∞∞−F , ( ),..., = 1+∞∞+F , dacă ii yx < atunci (adică este monotonă pe componente). Funcţia
definită de relaţia ( ) (...,,......, i yFxF ≤
( )kx,...,2
),...i Fxxf ,1
( ) ( )kk
k
k xxxxxx
Fxxxfxf ,...,,...
,...,,)( 2121
21 ∂∂∂∂
== (4.2.)
când derivata parţială există, se numeşte densitate de repartiţie multidimensională a lui X . Funcţia de repartiţie ( ) ( ) ( iiiii xFxFxXP ==+∞ )+∞+∞∞+=< ,...,,,,...,
iX se
numeşte funcţie de repartiţie marginală a lui , iar ( ) ( )iiii xFxf ′= ( când derivata există) se numeşte densitate de repartiţie marginală a lui . Se pot defini funcţii de repartiţie marginale pentru
iX1,...,3,2 −k componente etc. De exemplu, funcţia de repartiţie
1 corespunde componentelor , ale vectorului iX X X .$ Au loc relaţiile
( ) ( )
( )
( )∫
∫ ∫ ∫
∞−
∞− ∞− ∞−
=
==
==
x
x x kx
kk
k
duuf
dududuuuuf
xxxFxF1 2
2121
21
...,...,,...
,...,,
(4.3.)
Dacă variabilele , sunt independente atunci iX ki ≤≤1( ) ( ) ( ) ( )kkk xFxFxFxxxF ...,...,, 221121 = (4.4.)
sau analog pentru densităţi ( ) ( ) ( ) ( )kkk xfxfxfxxxf ...,...,, 221121 = (4.4'.)
Când există integralele
( ) ( )∫+∞
∞−
== iiiii dxxfxXEm 1 ,
( ) ( )∫ ∫+∞
∞−
+∞
∞−
== jijiijjijiij dxdxxxfxxXXEm
acestea se numesc momente de ordinul I, respectiv de ordimul II. Expresiile
( )( )[ ] ( )( )∫ ∫+∞
∞−
+∞
∞−
−−=−−=σ jijjiijjiiij dxdxmxmxmXmXE
când integralele există, se numesc momente centrate; în acest caz, se numeşte covarianţa lui cu şi se notează
ijσ
iX jX ( )jiij XXCov ,=σ . Se observă că
. Este uşor de arătat că dacă este independent de ,
( ) = ( ) σ= 2, iiiiii XVarXXCov=σ
jXiX
ji ≠ , atunci ( ) ( ) ( )jijijiij XEXEmmXXEm === ,
( ) 0, ==σ jiij XXCov . Este satisfăcută inegalitatea lui Schwarz, adică
( ) ( ) ( )jiji XVarXVarXXCov ≤, sau jiij σσ≤σ .
Se numeşte coeficient de corelaţie al variabilelor şi mărimea iX jX
( ) ( )( ) ( )ji
jijiij XVarXVar
XXCovXXCorr
,, ==ρ . (4.5.)
Se observă că şi el are următoarea interpretare: [ 1,1−∈ρ ] 0≈ρ ij atunci depinde foarte puţin de ; dacă 0 , atunci şi sunt invers corelate iar dacă , atunci şi sunt direct corelate; dacă 1
iX
jX
jX<ρ ij iX jX 0>ρ ij
iX ≈ρ ij atunci şi sunt puternic corelate.
iX jX
Pentru un vector aleator mediile definesc un vector numit
vectorul valoare medie notat im ( )′= kmmmm ,...,, 21
( )XE iar covarianţele ijσ definesc o matrice ijσ=Σ
(care este pozitiv definită, ) şi se notează 0Σ ( )XCov X ′=Σ , .
Presupunem că X are densitatea de repartiţie şi să considerăm
descopmunerea lui
)(xf
X în doi subvectori ( ) ( )( )′= 21 , XXX unde ( )1X este un vector iar ldimensionar − (2)X este un vector ldimensionaq − , kr < , . Definim
densitatea de repartiţie condiţionată kq =r +
( ) ( )( )2x1xf a lui ( )1X când ( ) x= ( )22X astfel
( ) ( )( ) ( )( ) ( )( ) ( )∫
∞+
∞−
=121
2121
,
,
dyxyf
xxfxxf , unde ( ) ( ) ( )( )21 , xxfxf = . (4.6.)
Cu această densitate se definesc, când există, ( ) ( ) ( )[ ]221 xXXE = ,
( ) ( ) ( ) ( ) ⎟⎠⎞
⎜⎝⎛ =
′ 2211 , xXXXCov adică media condiţionată, respectiv matricea de covarianţă
condiţionată. Uneori este necesar să simulăm ( ) ( ) ( )( )221 xXX = , când se cunoaşte repartiţia condiţionată. Metoda inversă este valabilă şi în cazul simulării vectorilor aleatori. Fie
( ) ( kkk xXxXxXPxxxFxF )<<<== ,...,,,...,,)( 221121 , funcţia de repartiţie a lui X şi fie funcţia de repartiţie marginală a lui , şi ( ) ( )1111 xXPxF <=
( ) ( )iiii xXxXPxxF1X
<<= ,...,,..., 111...1 funcţia de repartiţie marginală a lui ( )′iXX ,...,1 . Fie de asemenea funcţiile de repartiţie condiţionate
Să notăm cu inversa funcţiei 1−iF ( )ii xxxF ,...,, 21 în raport cu . Atunci are loc
următoatea teoremă. ix
Teorema 1. Dacă iU , ki ≤≤1 sunt variabile aleatoare uniforme 10 −
independente atunci vectorul aleator ( )′= YY 1 kYY ,...,, 2 ale cărui componente sunt ( )1
111 UFY −= , , …, ( )21
122 ,UYFY −= ( )iiiii UYYYFY ,,...,, 12
1−
−= , ki ≤≤2
(4.6'.)
are funcţia de repartiţie . ( )xFDe aici se deduce următorul algoritm pentru simularea vectorului aleator X când
se cunosc inversele ale funcţiilor de repartiţie condiţionate 1−iF ( )ii xxxF ,...,, 21 ix,1− în
raport cu . ix
Algoritmul AIMULT pentru metoda inversă multidimensională.
Generează ; Calculează random:U = ( )UFY 111−= ;
for to do begin 2:=i k
Generează random:=U ;
Calculează ( )UYYYFY iii ,,...,,: 1211
−−= ;
end.
IV.2. Simularea vectorilor uniformi
Vectorul aleator ldimensionak − ( )′= kVVVV ,...,, 21 are repartiţie uniformă pe domeniul mărginit kRD ⊂ dacă densitatea sa de repartiţie este
⎩⎨⎧
∉∈
=D vdacãD vdacã
0)(
kvf ,
mesDk 1= , ∫=
D
dvmesD (4.7.)
Rezultă că are repartiţie uniformă pe intervalul V [ ] [ ] [ kk bababaI ,...,, 2211 ×× ]×= , dacă +∞<<∞− i ba <i ki ≤≤1
( ) ( )
⎪⎪⎩
⎪⎪⎨
⎧
∉
∈−= ∏
=
I vdacã
vdacã
0
1
1
Iabvf
k
iii (4.7'.)
Din (4.7') rezultă că dacă este uniform pe V I atunci sunt variabile uniforme pe i independente deoarece
iV[ ii ba , ]
( ) ∏=
=k
iiik vfvvvf
121 ,...,, ,( ) ( )
[ ]
[ ]⎪⎩
⎪⎨
⎧
∉
∈−=
ii
iiiiii
ba
baabvf
, vdacã
, vdacã
0
1
i
i (4.7''.)
De aici se deduce următorul algoritm simplu pentru simularea vectorului IpeuniformV →~ .
Algoritmul VUNIFINT de simulare a vectorului uniform pe interval
Intrare , , ia ib ki ≤≤1 .
for to do begin 1:=i kGenerează random:=U ; Ia ;
end.
Pentru a simula un vector aleator W uniform pe kRD ⊂D ( D fiind un domeniu
oarecare), nu mai este adevărată (4.7''.). Dar să observăm că dacă şi V e uniform pe
I⊆I , atunci restricţia lui V la ( vectorul notat W ), este un vector uniform pe . De
aici rezultă următorul algoritm de respingere pentru simularea lui W . D D
Algoritmul VECTUNIFD de simulare a vectorului uniform pe domeniu
oarecare:
Construieşte un interval minimal [ ] [ ] [ kk bababaI ,...,, 2211 ×× ]×=
astfel încât ; (Acest interval există deoarece este mărginit).
ID ⊆ D
repeat Generează V uniform pe I ;
until ; DV ∈Ia VW =: .
Probabilitatea de acceptare a acestui algoritm este ( )∏
=
−== k
iii
a
ab
mesDmesImesDp
1
de unde se
deduce că pentru a obţine un interval I cât mai bun, acesta trebuie să fie minimal astfel încât . ID ⊆
IV.3. Simularea vectorilor normali
Vectorul aleator ldimensionak − X are repartiţia normală dacă densitatea sa este
( Σμ,N )
( )( ) [ ]
( ) ( )μ−−Σ′μ−−
Σπ=Σμ
xx
k exf1
21
21
2 det2
1,, , . kRx∈ (4.8.)
şi avem ( ) μ=XE , ( ) Σ=′XXCov , . (4.8'.)
Să notăm cu Z vectorul normal ( )IN ,0 unde I este matricea unitate dimensională iar 0 este vectorul nul din −k kR . Atunci densitatea a lui ( )zf Z
satisface proprietatea
( ) ( )∏=
=k
iii zfzf
1
, ( ) 2
2
21 iz
ii ezf−
π= (4.9.)
adică componentele ale lui iZ Z sunt independente şi normale ( )1,0N . Deci, simularea lui se realizează simplu astfel ( INZ ,0~→ ) Algoritmul VECNORMSTD de generare a vectorului normal standard.
for to do begin 1:=i kGenerează ( )1,0~ NT → ; Ia TZi =: ;
end.
Pentru a simula ( )Σμ→ ,~ NX vom folosi mai întâi teorema: Teorema 2. Dacă ( )Σ→ ,0~ NY şi este o matrice C kk × iar este un vector,
atunci μ
1×k ( )CCNCYX ′Σμ→+μ= ,~ ,
Fiind dată matricea (pozitiv definită) se ştie că există o matrice 0Σ0ijcC = inferior triunghiulară (adică =ijc , ji < ), astfel încât $ . Matricea
este matricea Choleski. De aici deducem teorema: Σ=′CC
C Teorema 3. Dacă ( )Σ→ ,0~ NZ , μ este un vector 1×k şi este matricea Choleski a lui , atunci
CΣ ( )Σμ→+μ= ~CZ ,NX .
Matricea C se determină cu ajutorul matricii ijσ=Σ cu formulele lui Choleski:
11
11 σ
σ= i
ic , ki ≤≤1
∑−
=
−σ=1
1
2i
ririiii cc , ki ≤<1
jj
j
rjririj
ij c
ccc
∑−
=
−σ=
1
1 , kij ≤<<1
0=ijc , kji ≤<<1 .
(4.10.)
Simularea vectorului ( )Σμ→ ,~ NX , după ce se calculează într-un pas pregătitor matricea C a lui Choleski, se realizează cu următorul algoritm: Algoritmul NORMMULT de simulare a vectorului normal
Uneori este necesar să simulăm un subvector al unui vector normal când se dă
celălalt subvector. Mai precis dacă ( ) ( )( ) ( Σμ→ )′= ,~, 21 NXXX , cu ( )1X vector
şi ldimensionar − ( )2X vector ldimensionaq − , kqr =+ , se cere să se simuleze ( )1X condiţionat de . ( ) )2 (2x=X
Să considerăm scrierea celulară a lui μ şi Σ , adică ( )
( ) ⎟⎟⎠
⎞⎜⎜⎝
⎛
μμ
=μ 2
1
, ⎟⎟⎠
⎞⎜⎜⎝
⎛ΣΣΣΣ
=Σ2221
1211 (4.11.)
Dacă μ şi sunt cunoscute atunci este adevărată următoarea teoremă. Σ Teorema 4. Repartiţia condiţionată a lui ( )1X când se dă este normală
( ) ( )22 xX =( )( )*
111
* ,ΣμN unde ( ) ( ) ( ) ( )( )221
221211
* μ−ΣΣ+μ=μ − x
211
221211*11 ΣΣΣ−Σ=Σ −
(4.12.)
De aici se deduce algoritmul de simulare a lui ( )1X condiţionat de ( ) ( )22 xX = şi anume: Algoritmul NORMCOND de simulare condiţionată.
Intrare , , ( )2x iμ ijΣ , 2,1, =ji ;
Calculează , ; ( )1*μ
*11Σ
Generează ( ) ( )rrINZ ×→ ,0~1̀ ;
Calculează astfel încât *C +Σ=′
11**CC ;
Calculează ( ) ( ) ( )1*1*
1* ZCX +μ= .
Rezultatul este ( ) ( )( )*
111
*1
* ,~ Σμ→ NX .
IV.4. Simularea repartiţiei Cauchy multidimensionale Repartiţia Cauchy standard care are densitatea:
( )211)(
xxf
+π= (4.13.)
se poate simula cu ajutorul metodei inverse. Dacă , sunt variabile normale independente, atunci variabila
1Z 2Z )1,0(N
2
1
ZZX = (4.14.)
are repartiţia Cauchy standard. Din această proprietate rezultă deci o metodă simplă de simulare a variabilei Cauchy standard unidimensionale ca raport de variabile normale
independente. )1,0(N O variabilă Y are repartiţia ( )σ,cCauchy dacă densitatea sa este
( )2
2
1
11)(
σ−
+⋅
σ=
cxxg , Rc∈ , 0<σ ,
(4.15.)
Simularea lui Y se va face cu formula XcY σ+= . În mod asemănător se poate introduce repartiţia Cauchy standard multidimensională. Vectorul aleator X are repartiţia Cauchy standard multidimensională dacă densitatea sa de repartiăie e te de forma s
( )( ) 2
12
12
1
1
1)( ++
+
′+⋅
π
Γ= kk
k
xxxf , . kRx∈ (4.14'.)
Dacă Z este un vector ldimensionak − normal şi este o variabilă normală ) independentă de
),0( IN η1,0(N Z , atunci X se poate simula cu formula
η=
ZX . (4.14''.)
Se poate considera şi repartiţia multidimensională care are densitatea de forma
( Σ,cCauchy )
( )( ) ( ) 2
1
121
21
1
1
det)( +
−
+
+
⎥⎦⎤
⎢⎣⎡ −Σ′−+
⋅Σπ
Γ= kk
k
cxcx
xg (4.15'.)
Relaţia dintre vectorul Cauchy standard X şi vectorul este asemănătoare relaţiei dintre vectorul
( Σ→ ,~ cCauchyY )Z normal şi vectorul ),0 I(N Y normal ( )Σ,cN ,
adică dacă X are repartiţie Cauchy standard şi ( )Σ,cCauchy→~Y şi considerăm matricea inferior triunghiulară astfel ca S SS ′=Σ , atunci Y se poate genera cu formula
cSXY += (4.15''.) Combinând (4.14') cu (4.15'') rezultă că ( )Σ→ ,~ cCauchyY se simulează cu formula
η=
*ZY (4.16.)
unde , independent de ( Σ→ ,0~ NZ ) )1,0(~ N→η .
IV.5. Simularea repartiţiei multinomiale
Această repartiţie este generalizarea multidimensională a repartiţiei binomiale.
Vectorul aleator cu componente întregi nenegative ( )′= kXXXX ,...,, 21 , are repartiţia multinomială nXXX k =+++ ...21 ( )kpp ,...,2pnMultinom ,, 1 dacă
( ) knk
nn
kkk ppp
nnnnnXnXnXP ...
!!...!!,...,, 2
21
121
2211 ==== ,
, , nnnn k =+++ ..21 0>ip ki ≤≤1 , 1..21 =+++ kppp (4.17.)
Vectorul X are o interpretare asemănătoare variabilei binomiale în cazul unui experiment cu urnă după cum urmează: să presupunem că într-o urnă se găsesc bile de culoarea i , ,
iNk≤1 i≤ kNNNN +++= ...21
iX. Se presupune că se extrag bile din
urnă cu întoarcere, din care sunt bile de culoarea i , n
kXXXn +++= ...21 . Atunci
vectorul aleator ( )′= X kXX ,...,, 21X are repartiţia ( )kppnMultinom ,...,, 1 , iNip = , . Această interpretare ne conduce la următorul algoritm pentru simularea lui 1=∑i ip
X .
Algoritmul MULTINOM de simulare a repartiţiei multinomiale
Intrare şi calculează , kpppkn ,...,,,, 21 [ ] ∑α
=
=α1i
ipF ki ≤≤1 ;
(Acesta este un pas pregătitor); Iniţializează [ ] 0:1 =X , …, [ ] 0=kX ;
for to do begin 1:=i nSimulează randomU := ; Ia 1: += ii ; while [ ]iFU ≥ do 1+= ii ; Insumează [ ] [ ] 1: += iXiX ; end.
Testarea algoritmului se face în mod simplu astfel:
- se generează selecţia de volum T , ( ) ( ) ( )( )′= αααkXXX ,...,1 , T≤α≤1 ;
- se calculează mediile aritmetice şi dispersiile empirice ( )
( )[ ] 2
1
2
1
11i
T
ii
T
ii XXT
XT
X −== ∑∑=
α
=α
α
Pentru validarea algoritmului, trebuie ca pentru un T mare să avem ii Xm ≈ , , unde şi sunt date de (4.17’.). 22
ii s≈σ im 2iσ
IV.6. Simularea repartiţiei Dirichlet.
Un vector aleator X are repartiţia ( )121 ,...,, +υυυ kDirichlet dacă densitatea sa de
repartiţie este
( )( )( ) ( ) ( )
⎪⎩
⎪⎨
⎧
∉
∈−−−υΓυΓυ++υΓ
=+υυυ
+
+
k
k1
11
111
11
1
S xdacã
S xdacã
0
...1.........
,...,k
kk
kk
k
k
xxxxxxf (4.18.)
dac , ( )ii GammaY υ→ ,1,0~ 11 +≤≤ ki sunt variabile gama independente, atunci vectorul X ale cărui componente sunt de forma
121 ... ++++=
k
ii YYY
YX , ki ≤≤1 (4.19.)
are o repartiţie ( 121 ,...,,Dirichlet + )υυυ k . Se observă deci că repartiţia Dirichlet este o extensie la cazul multidimensional a repartiţiei Beta. Formulele (4.15) simulează componentele unui vector Dirichlet, presupunându-se desigur că parametrii , 1υ ki ≤≤1 sunt cunoscuţi.
LISTA DE SUBIECTE
1. Modelul matematic al unui model de simulare 2. Clasificari ale modelelor matematice 3. Sistem de asteptare- model mathematic 4. Structura unui model de simulare 5. Agenda simularii 6. Realizari iterative ale sistemului 7. Proprietati ale unei masini secventiale 8. Sistem cu evenimente discrete (SED) 9. Proprietati ale extensiei unei functii de tranzitie autonoma a unui SED 10. Proprietati ale functiei ce da timpul total necesar unui SED pentru a face n tranzitii 11. Proprietati ale numarului maxim de tranzitii facute intr-un SED intr-un interval de
timp dat 12. Proprietati ale unui SED viabil 13. Etapele realizarii unui experiment de simulare 14. Generalitati despre limbajul GPSS 15. Instructiuni GPSS 16. Program GPSS pentru simularea unui sistem de asteptare cu o statie 17. Program GPSS pentru simularea unui sistem de asteptare cu statii paralele 18. Model de simulate cu statii paralele si preferinte 19. Model de simulare a activitatilor de reparatie a unui automobil intr-un atelier service 20. Repartitia uniforma 21. Necesitatea simularii numerelor aleatoare 22. Metod congruentiale liniare 23. Teorema ce da perioada maxima a unui generator congruential liniar 24. Generatori de numere uniforme 25. Amestecarea de generatori si generarea unui indice aleator 26. Metoda inversa pentru simularea variabilelor neuniforme 27. Metoda compunerii pentru simularea variabilelor neuniforme 28. Metoda respingerii pentru simularea variabilelor neuniforme 29. Teoreme ce fundamenteaza algoritmii de respingere pentru simularea variabilelor
neuniforme 30. Teorema lui Forsythe (a sirului descendent) 31. Algoritmii Compdiscr si Compcont 32. Algoritmii Lapcomp si Algres 33. Algoritmii Gamresm1, Resp34 si Respsird 34. Metode de simulare care se obtin prin transformarea unor variabile uniforme 0-1 35. Stabilitatea familiilor de variabile aleatoare gamma standard si normale 36. Simularea repartitiilor χ2 si Student 37. Simularea repartitiilor Snedecor si lognormala 38. Simularea repartitiei exponentiale 39. Algoritmul Exrj 40. Simularea repartitiei gamma 41. Algoritmul Gamcomnl1, Gamresnl1 si Gamresnl2 42. Metode pentru simularea variabilei gamma cu parametru supraunitar 43. Algoritmul Gamcauchy 44. Metode pentru simularea repartitiei beta 45. Algoritmul Beta13 46. O metoda de compunere-respingere pentru simularea repartitiei normale 47. Algoritmul Rjnorm 48. Metoda polara pentru simularea repartitiei normale 49. Simularea unor repartitii bazate pe probe Bernoulli 50. Simularea repartitiei binomiale 51. Algoritmii Bern si Binnorm
52. Simularea repartitiei Pascal si algoritmul Pascal 53. Simularea repartitiei geometrice 54. Simularea repartitiei hipergeometrice 55. Algoritmul Hipergeom 56. Simularea repartitiei Poisson 57. Validarea generatorilor 58. Constructia histogramei 59. Algoritmul Histogramei 60. Testul de concordanta χ2 61. Testarea unui generator care simuleaza o variabila neuniforma 62. Vectori aleatori (generalitati) 63. Teorema ce reprezinta metoda inversa pentru vectori aleatori 64. Simularea repartitiei Gumbel bidimensionale 65. Algoritmul Gumbel si Aimult 66. Simularea vectorilor uniformi 67. Algoritmii Vunifint si Vectunifd 68. Simularea vectorilor normali 69. Algoritmii Vectnormst si Normmult 70. Algoritmul Normcond de simulare conditionata 71. Simularea repartitiei Cauchy multidimensionale 72. Simularea repartitiei multinomiale 73. Algoritmul Multinorm 74. Simularea repartitiei Dirichlet
