Felsőbb Mathematica házi feladat

Készítette: Temesváry András (U6QW38)

2006/2007 I.félév, 2006.12.06.

Diszkrét idejű Markov láncok tranziens-és határeloszlásának vizsgálata z-transzformáció segítségével

Bevezetés

A kommunikációs hálózatok megjelenésével egyidőben fontos kérdéssé vált ezen hálózatok teljesítményének elemzése. Számos kommunikációs hálózatbeli protokoll modellezhető diszkrét idejű Markov láncok segítségével, így fontossá vált ezen modellek minél könnyebb számítását lehetővé tevő módszerek bevezetése. Különösen fontos lehet, hogy eg diszkrét idejű Markov lánc tranziens eloszlását könnyedén, minél kisebb számítási igénnyel határozzuk meg, lehetővé téve a rendszer érdekes állapotainak valószínűségének meghatározását. Szintén fontos lehet a rendszert leíró modell határeloszlásának meghatározása is, hogy képünk legyen a rendszer stabilitását illetően. 
      A következőkben ismertetem a diszkrét idejű Markov láncok tranziens és egyensúlyi eloszlásának analíziséhez használt z-transzformációt, ennek néhány fontos tulajdonságát, használhatóságát a diszkrét eloszlásokkal kapcsolatban, valamint egy példán keresztül megoldást mutatok a z-transzformációt használó eljárás használatára a Mathematica segítségével. Elmondható, hogy már a bemutatott példa is komoly számítást igényel, így a matematikai segédprogramok használata ezen problémára szinte elkerülhetetlen általános esetben. Ahol lehet a levezetésekben, illetve a példa bemutatása során természetesen minden számításnál a Mathematica segítségét veszem igénybe.

1. fejezet - Definiciók, tulajdonságok

\(A\ diszkrét\ idejű\ z - 
      transzormáció, \ 
    vagy\ más\ elnevezéssel\ a\ generátor\ függvény\ általános\ definíciója\ a\ \(következő : 
f \((z)\)\)\  = \ ∑\+\(k = \ 0\)\%∞\( f\_k\) 
          z\^k\)
\(ahol\ {f\_i} - k\ egy\ valós\ sorozat\ elemei, \ 
    f\_i ∈ 
      R\ minden\ i - re . 
A\ továbbiakban\ a\ z - 
        transzformációt\ diszkrét\ valószínűségi\ változókra\ alkalmazzuk . \ 
            A\ diszkrét\ valószínűségi\ változókat\ az\ eloszlásfüggvényükkel\ adjuk\ meg, \ 
    melyre\ igaz\ a\ \(következő : 
∑\+\(k = \ \(-∞\)\)\%∞ p\_k\) = \ 
      1, \ ahol\ p\_k = P \((X = k)\), \ és\ 0 ≤ 
      p\_k ≤ 1\)
\(tehát\ \(p\_k\) 
      annak\ a\ valószínűsége, \ 
    hogy\ az\ X\ valószínűségi\ változó\ éppen\ a\ "k"\ értéket\ veszi\ fel\ egy\ adott\ kísérlet\ \(\(eredményeképpen\)\(.\)\(\ \)\)\)

Valószínűségi változók összegének kiszámítása:

\(X\ és\ Y\ egy - 
      egy\ valószínűségi\ változó, \ 
    az\ összegüket\ szeretnénk\ meghatározni, \ 
    ezt\ jelöljük\ Z - vel\ \((Z\  = \ X + \ Y)\), \ 
    ekkor : 
P \((Z = k)\) = \(∑\+\(n = \ 0\)\%k 
            P \((X = n)\) P \((Y = k - n)\) = 
        X*Y\)\)
\(a\ két\ valószínűségi\ változó\ konvolúciója, \ 
    amely\ a\ legtöbb\ esetben\ nehézkesen\ kezelhető\ kifejezés . 
Azonban\ ha\ vesszük\ az\ egyenlet\ mindkét\ oldalának\ a\ z - 
      transzformáltját, \ úgy\ a\ következő\ kifejezéshez\ \(jutunk : 
∑\+\(k = \ 0\)\%∞\( z\^k\) 
              P \((Z = 
                  k)\)\) = ∑\+\(k = \ 0\)\%∞\(∑\+\(n = \ 0\)\%k P \((X = n)\) P \((Y = k - n)\) \(z\^n\) 
            z\^\(k - n\)\)\ \)
\(amely\ kifejezés\ kicsit\ átrendezve\ és\ az\ összegzési\ határokat\ \(átalakítva : 
∑\+\(k = \ 0\)\%∞\( z\^k\) P \((Z = k)\)\) = \((∑\+\(n = \ 0\)\%∞ 
              P \((X = n)\) z\^n)\) \((∑\+\(k = \ n\)\%∞ 
              P \((Y = k - n)\)\ z\^\(k - n\))\)\ \)
\(ahonnan\ látszik, \ 
    hogy\ a\ diszkrét\ valószínűségi\ változók\ z - 
        transzformáltjaira\ az\ összegképzésnél\ a\ konvolúció\ helyett\ egy\ nagyon\ egyszerű\ szorzat\ alakot\ \(kapunk : 
\(p\_Z\) \((z)\)\) = \(p\_X\) \((z)\) \(p\_Y\) \((z)\)\)

\(\(\(
\)\(\(Határértéktétel, \ 
    mely\ közvetlen\ alkalmazható\ lesz\ a\ diszkrét\ idejű\ Markov\ láncok\ egyensúlyi\ eloszlásának\ \(meghatározására : 
p\_∞\) = 
      lim\_\(z -> \ 1\)\((1 - z)\) p \((z)\)\)\(
\)
    \)\)\)

2. fejezet - Diszkrét valószínűségi változók z-transzformáltjai, momentumainak kiszámítása

Általános tulajdonságok:

A z-transzformáltak segítségével könnyen kiszámíthatjuk egy valószínűségi változó momentumait. A legfontosabb momentum a vizsgálatok során az első momentum, azaz a valószínűségi változó várható értéke. Kiszámítására a következő formulát lehet alkalmazni:

ExpectedValue[distributionZtransform_] := 
      D[distributionZtransform, z] /. z -> 1

A valószínűségi változó várható értékét tehát úgy kaphatjuk meg, hogy kiszámítjuk az eloszlásfüggvény z-transzformáljának deriváltját, és ezt a kifejezést a z=1 helyen kiértékeljük, hiszen

p[z_] := (∑_(k = 0)^∞) (Subscript[p, k]) (z^k)

derivative1 = D[p[z], z]

\(És\ ha\ a\ fenti\ kifejezésbe\ z = 
      1 - et\ helyettesítünk, \ úgy\ éppen\ \(\(\(\ \)\(∑\)\)\+\(k = 0\)\%∞\) k\ p\_k - 
      t\ kapunk, \ 
    ami\ éppen\ megegyezik\ a\ várható\ érték\ \(\(definíciójával\)\(:\)\)\)

derivative1 /. z -> 1

Érdekes lehet még továbbá a diszkrét valószínűségi változók második momentumának kiszámítása, hiszen számos modell arra a feltételezésre épít, hogy ez az érték ismert (pédául az általános M/G/1-es rendszerek analízisénél). A számítás menete a következő:

derivative2 = D[p[z], {z, 2}]

Hogyha a fenti kifejezést a z=1 helyen vesszük, akkor a következő alakot kapjuk:

derivative2 /. z -> 1

Amely kifejezés éppen a valószínűségi változó várható értékével kisebb a valószínűségi változó második momentumánál. Így tehát a következő formulát kaptuk:

E[(X^2)]=((d^2)/(dz^2)) p (z)\( | \_\(z = 1\)\)+E[X]

momentum2 = 
      FullSimplify[(D[p[z], {z, 2}] /. z -> 1) + (D[p[z], z] /. 
              z -> 1)]

A második momentum kiszámítására való formula tehát, mely később is használható lesz:

SecondMomentum[
        distributionZTransformal_] := (D[distributionZTransformal, {z, 2}] /. 
            z -> 1) + (D[distributionZTransformal, z] /. 
            z -> 1)

A következő diszkrét valószínűségi változók z-transzformáltjait fogjuk vizsgálni:

A binomiális eloszlás:

dBinomial[n, k] := 
      Binomial[n, k] (p^k) (((1 - p))^(n - k))

A binomiális eloszlás z-transzformáltja:

dBinomialZtrans = (∑_(k =  0)^∞) dBinomial[n, k] 
          (z^k)

A z-tranformált alapján a várható érték:

FullSimplify[ExpectedValue[dBinomialZtrans]]

Ahogy láthatjuk a binomiális eloszlás várható értéke éppen (n p).

A binomiális eloszlás második momentuma:

FullSimplify[SecondMomentum[dBinomialZtrans]]

A geometriai eloszlás:

dGeometry[k_] := (((1 - p))^(k - 1)) p

A geometriai eloszlás z-transzformáltja:

dGeometryZtrans = (∑_(k =  
              0)^∞)((dGeometry[k]) 
            (z^k))

A geometriai eloszlás várható értékének kiszámítása a z-transzformált segítségével:

ExpectedValue[dGeometryZtrans]

\(A\ várható\ érték\ tehát\ \(\(1\/p\)\(.\)\)\)

A geometriai eloszlás második momentuma:

FullSimplify[SecondMomentum[dGeometryZtrans]]

A Poisson eloszlás:

dPoisson[λ, 
        k] := ((λ^k)/Factorial)[
              k] (E^(-λ))

A Poisson eloszlás z-transzformáltja:

dPoissonZtrans = (∑_(k =  0)^∞) dPoisson[λ, 
            k] (z^k)

A z-transzformált alapján a várható érték kiszámítása:

ExpectedValue[dPoissonZtrans]

A várható érték tehát megegyezik a Poisson eloszlás paraméterével, λ-val.

A Poisson eloszlás második momentuma, amelyre sok valóságban használt modell támaszkodik:

FullSimplify[SecondMomentum[dPoissonZtrans]]

Megmutatható továbbá, hogy a binomiális eloszlás határértéke (n->∞,p->0 esetén)éppen a Poisson eloszlással egyezik meg:

FullSimplify[
      Limit[Binomial[n, 
            k] ((((λ/n)))^k) (((1 - (λ/n)))^(n - k)), 
        n -> ∞]]

Sajnos a Mathematica a nevezőben valamiért k! helyett a Gamma[k+1] függényt számította ki, egyébként a számítás helyes lenne.

3. fejezet - Az inverz z-transzformáció

\(Inverz\ z - 
      transzfomáción\ azt\ a\ módszert\ értjük\ melynek\ során\ a\ z - 
      transzformáltból\ visszaszámítjuk\ az\ egyes\ \(p\_k\) együtthatókat . 
A\ z - 
      transzformált\ általános\ alakja\ már\ szerepelt, \ 
    most\ \(megismételjük : 
p[
            z]\) := ∑\+\(k = 0\)\%∞\( p\_k\) 
          z\^k\)
\(Ez\ alapján\ a\ \(p\_k\) 
      együtthatók\ kiszámításának\ \(\(módja\)\(:\)\)\)

p0 = FullSimplify[p[z] /. z -> 0]

p1 = FullSimplify[D[p[z], z] /. z -> 0]

p2 = FullSimplify[D[p[z], {z, 2}]/2 /. z -> 0]

\(Látszik\ tehát, \ 
    hogy\ általánosan\ a\ következő\ összefüggés\ adható\ az\ \(együtthatókra : 
p\_k\) = \(\(\(1\/\(k!\)\) \(d\^k\/dz\^k\) p[z]\)\( | \_\(z = 0\)\)\)\)

Definiáljuk az általános képletet, hogy használható legyen bármikor:

pcoeff[k_] := (1/Factorial)[k] D[p[z], {z, k}] /. 
        z -> 0

\(\(\(Így\ ellenőrzésképpen\ a\ 0. \ együttható\ \((p\_0)\)\)\(:\)\)\)

pcoeff[0]

Az 1.  együttható ((Subscript[p, 1])):

pcoeff[1]

És a 2. együttható:

pcoeff[2]

Megjegyzés: ezek az eredmények ugyan helyesek, csak sajnos nincsenek kiértékelve. Sajnálatos módon ezen a ponton nem találtam meg a Mathematicában azt a lehetőséget, ami ezeket a végtelen összegzéseket szépen kiértékelné.

4. fejezet - Diszkrét idejű Markov láncok tranziens-és határeloszlásának vizsgálata

A példa definiálása:

Példaként tekintsünk egy két állapotú diszkrét idejű Markov láncot, melyet a következő egylépéses átmenetvalószínűség-mátrix ír le:

PiMatrix = {{0.5, 0.5}, {0.75, 0.25}}

MatrixForm[PiMatrix]

Ez a példa egyszerű,de kézi számítása már komoly nehézségeket is okozhat, általános esetben  pedig a probléma megoldása számítógépes segítség nélkül nagyon nehéz volna.

Definiáljuk továbbá a 2x 2-es egységmátrixot is, mert a számítások során szükség lesz rá:

IdMatrix = IdentityMatrix[2]

MatrixForm[IdMatrix]

\(A\ rendszer\ állapotát\ a\ n . \ 
        diszkrét\ időpillanatban\ a\ \(P\^\((n)\)\) 
      eloszlásvektor\ jellemzi, 
    melynek\ elemszáma\ a\ rendszer\ állapotainak\ számával\ egyezik\ meg . \ 
        Tranziens\ eloszlás\ vizsgálatáról\ beszélhetünk, \ 
    ha\ ezt\ a\ valószínűségi\ eloszlás\ vektort\ szeretnénk\ egy\ általunk\ választott\ \(n\_i\) 
      időpillanatban\ meghatározni, \ 
    illetve\ határeloszlásról\ beszélünk, \ 
    ha\ a\ rendszer\ állapotát\ a\ végtelenben, \ 
    n -> ∞\ esetén\ kívánjuk\ meghatározni\)
\(A\ diszkrét\ rendszer\ állapotát\ minden\ időpontban\ meghatározhatjuk\ a\ kezdeti\ eloszlás\ \((P\^\((0)\))\)\ és\ az\ egylépéses\ átmenetvalószínűség - 
        mátrix\ segítségével\), hiszen:
(P^((n)))=(P^((0))) (PiMatrix^n)

A tranziens-és határeloszlás meghatározása az egylépéses átmenetvalószínűség-mátrix n. hatványának segítségével:

Meghatározhatjuk tehát a tranziens -és határeloszlást az egylépéses átmenetvalószínűség-mátrix n. hatványának segítségével:

PiMatrixOnN = FullSimplify[MatrixPower[PiMatrix, n]]

Ebből a tranziens eloszlást könnyedén megkaphatjuk egy egyszerű vektor-mátrix szorzással (ehhez azonban meg kellett határozni az átmenetvalószínűség-mátrix n. hatványát):

{Pn0, Pn1} = {P00, P01} . PiMatrixOnN

A 0. állapot valószínűsége tehát az n. időpontban:

Pn0

Az 1. állapot valószínűsége az n. időpontban:

Pn1

A fenti kifejezésben {Pn0,Pn1} jelöli az n. időpillanatban a rendszert leíró eloszlásvektort, tehát az egyes állapotok valószínűségét ebben a időpontban. A {P00,P01} vektor pedig a rendszer kezdeti eloszlása,a rendszer indításakor az egyes állapotok valószínűsége.

A határeloszlás meghatározása ({Pinf0,Pinf1}):

Pn0

Pinf0 értékét Pn0 kifejezés n->∞-ben felvett értéke adja meg:

Pinf0 = 
      Limit[(0.4 (((-0.25))^n) + 
                0.6 (1^n)) P00 + (-0.6 (((-0.25))^n) + 
                0.6 (1^n)) P01, n -> ∞]

Pn1

Hasolóan Pinf1 értékét a Pn1 kifejezés n->∞-ben felvett értéke határozza meg:

Pinf1 = 
      Limit[(-0.4 (((-0.25))^n) + 
                0.4 (1^n)) P00 + (0.6 (((-0.25))^n) + 
                0.4 (1^n)) P01, n -> ∞]

Most tehát megmutattuk, hogy hogyan függ a rendszer állapotainak eloszlása a kezdeti eloszlástól. A határeloszlás létezik, a rendszer stabil.

A fent bemutatott eljárásnál azonban a vizsgálatokhoz egy könnyebben kezelhető, egyszerűbb számítási módszerre van szükség, amit a z-transzformáltakkal való számítás tesz lehetővé.

A tranziens eloszlás meghatározása z-transzformáltak segítségével:

Kiszámíthatjuk a rendszer állapotát leíró eloszlásvektor z-transzformáltját, amelyet P(z)-vel jelölünk:

P (z) = (∑_(k = 0)^∞) (P^((k))) 
            (z^k) = (P^((0))) + (∑_(k = 1)^∞) 
                    (P^((k - 1))) PiMatrix (z^k) = (P^((0))) + 
              PiMatrix z (∑_(k = 1)^∞) (P^((k - 1))) 
                    PiMatrix (z^(k - 1)) = 
            (P^((0))) + 
              P (z) PiMatrix z
\(ahonnan\ P \((z)\) - 
        t\ kifejezve\ a\ következő\ formulát\ \(kapjuk : 
P \((z)\)\) = \(P\^\((0)\)\) \((IdMatrix - PiMatrix\ z)\)\^\(-1\)\)

\(A\ fenti\ számításokban\ kiszámolt\ inverz\ mátrixos\ kifejezés, \(\(\((IdMatrix - 
              PiMatrix\ z)\)\^\(-1\)\)\(:\)\)\)

InvMatrixExpr = Inverse[IdMatrix - z PiMatrix]

MatrixForm[InvMatrixExpr]

Jól jöhet,ha tudjuk ennek a mátrixnak a determinánsát, például a gyökök meghatározásához.

DetInvExpr = Det[(IdMatrix - z PiMatrix)]

A gyököket ebből könnyeb meghatározhatjuk általános esetben is, melyek kellhetnek az egyes számításokhoz:

Solve[DetInvExpr == 0, z]

Ha elvégzünk a fent bemutatott formulákon néhány átalakítást, akkor a következő jól használható alakot kapjuk:

\(P \((z)\) = \(\(P\^\((0)\)\) \(∑\+\(n = 0\)\%∞\( 
                  PiMatrix\^\((n)\)\) 
              z\^n\) = \(\(P\^\((0)\)\) \((IdMatrix - PiMatrix\ z)\)\^\(-1\) = \(P\^\((0)\)\) InvMatrixExpr\)\)\)

\(Azaz\ ez\ azt\ jelenti, \ 
    hogy\ amennyiben\ sorba\ fejtjük\ az\ InvMatrixExpr - 
      nek\ nevezett\ \((IdMatrix - PiMatrix\ z)\)\^\(-1\)\ mátrixot, \ úgy\ éppen\ az\ egylépéses\ átmenetvalószínűség - 
      mátrix\ n . \ 
          haványának\ alakját\ kaphatjuk\ meg, 
    ezzel\ nagy\ mértékben\ csökkentve\ a\ módszer\ számításigényét . \ 
          A\ fenti\ képletben\ k\ helyett\ n - 
      et\ használtam\ az\ összegzés\ indexeként, \ hogy\ jobban\ kitűnjön, \ 
    hogy\ mit\ is\ szeretnénk\ elérni\ ezzel\ az\ \(\(átalakítással\)\(.\)\)\)

Az inverz mátrixos kifejezés egyszerűbb alakra hozva:

FullSimplify[InvMatrixExpr]

Az inverz mátrixos kifejezés részlettörtekre bontva:

Apart[InvMatrixExpr]

InvMatrixExpr1 = Apart[FullSimplify[InvMatrixExpr]]

Naiv módszer az együtthatók kiszámítására, ha a kiszámolt gyökök alapján felírjuk a következő alakot:
(((IdMatrix - PiMatrix z))^(-1)) = (A/(1 - z)) + (B/(1 + 0.25  z)) , hiszen a gyökök éppen {z = 1} és {z = -4} voltak.

Ezek alapján például az A mátrix a11 elemének a kiszámítása a következő lehet:

Part1 = Simplify[Part[InvMatrixExpr1[[1, 1]], 1]]

Solve[Part1 == (ElementA11/(1 - z)), ElementA11]

Az A11-es elem valóban 0.6, mint azt már az egylépéses átmenetvalószínűség-mátrix n. hatványának segítségével megállapítottuk. A többi mátrixelem meghatározása a naiv módszerrel hasonlóan történik.

Clear[A, B]

Az A és B mátrix kiszámítására általános esetben lineáris egyenletrendszert kell megoldanunk, melynek menete a következő:

A = {{a11, a12}, {a21, a22}}

B = {{b11, b12}, {b21, b22}}

Felírva a megfelelő egyenleteket, és az egyenletrendszert megoldva:

FullSimplify[
      Solve[{InvMatrixExpr[[1, 1]] == 
            (a11/(1 - z)) + (b11/(1 + 0.25  z)), 
          InvMatrixExpr[[1, 2]] == 
            (a12/(1 - z)) + (b12/(1 + 0.25  z)), 
          InvMatrixExpr[[2, 1]] == 
            (a21/(1 - z)) + (b21/(1 + 0.25  z)), 
          InvMatrixExpr[[2, 2]] == 
            (a22/(1 - z)) + (b22/(1 + 0.25  z)), a11 + b11 == 1, 
          a12 + b12 == 0, a21 + b21 == 0, 
          a22 + b22 == 1}, {a11, a12, a21, a22, b11, b12, b21, 
          b22}]]

FullSimplify[
      Solve[{InvMatrixExpr[[1, 1]] == 
            (a11/(1 - z)) + (b11/(1 + 0.25  z)), 
          InvMatrixExpr[[1, 2]] == 
            (a12/(1 - z)) + (b12/(1 + 0.25  z)), 
          InvMatrixExpr[[2, 1]] == 
            (a21/(1 - z)) + (b21/(1 + 0.25  z)), 
          InvMatrixExpr[[2, 2]] == 
            (a22/(1 - z)) + (b22/(1 + 0.25  z)), a11 + b11 == 1, 
          a12 + b12 == 0, a21 + b21 == 0, a22 + b22 == 1, 
          a11 + a12 == 1, a21 + a22 == 1}, {a11, a12, a21, a22, 
          b11, b12, b21, b22}]]

A = {{0.6, 0.4}, {0.6, 0.4}}

MatrixForm[A]

B = {{0.4, -0.4}, {-0.6, 0.6}}

MatrixForm[B]

Innen a tranziens eloszlás megállapítása már ugyanúgy történhet, mintha ismernénk az átementvalószínűség-mátrix n. hatványát, amelyet az előzőekben mutattunk be.

Pn0

Pn1

A határeloszlás számítása z-tranformáció segítéségével:

A diszkrét rendszer határeloszlását az első fejezetben már bevezetett határértéktétel segítségével rögtön megkaphatjuk:

{Pinf0, Pinf1} = 
      Limit[(1 - z) {P00, P01} . InvMatrixExpr, z -> 1]

Ahol {Pinf0,Pinf0} a  rendszer határeloszlását adja meg.

Így tehát a 0. állapotban tartózkodás valószínűsége a végtelenben:

Pinf0

Illetve az 1. állapotban való tartózkodás valószínűsége a végtelenben:

Pinf1

Az így meghatározott határeloszláshoz nem kellett ismerni az egylépéses átmenetvalószínűség-mátrix n. hatványát, így ez a módszer nagyon hatékonyan vezet az eredményhez.

Összefoglalás:

A kommunikációs hálózatok tervezése és analízise során sokszor diszkrét idejű Markov láncok írják le a rendszer működését, állapotát. Két lehetőségünk van ha ezen rendszerek tranziens-és határeloszlását szeretnénk meghatározni. Költséges módon meghatározhatjuk a diszkrét rendszert leíró egylépéses átmenetvalószínűség-mátrix n. hatványát,és erre alapozhatjuk a számításokat, illetve segítségül hívhatjuk a z-tranformált aparátust, mellyel a probléma jóval könnyebben, kisebb komplexitással kezelhető. A dolgozatban ismertettem a z-transzformáltak alapvető tulajdonságait, alkalmazhatóságát valószínűségi változók tulajdonságainak kiszámítására, valamint egy példán keresztül megmutattam a z-transzformáltakat használó módszert a rendszer tranziens-és határeloszlásának meghatározására.
A továbbiakban szeretnék még egy olyan általános módszert adni, amely egy mátrixot sorba tud fejteni beépített eljárások segítségével, a Mathematica használatával, ezzel egyszerűsítve az egylépéses átmenetvalószínűség-mátrix n. hatványának meghatározását.

Irodalomjegyzék:

dr .  Jereb László  -  
        Teljesítményelemzés, egyetemi jegyzet,előadásjegyzet, BME 2006. 
\(Mathematica\ dokumentáció\  - \ beépített\ segédletek, \ 
    dokumentációk\)