2005/06 õszi félévi témalabor prezentációs segédanyag

Az alábbiakban a témalabor elõadásomhoz készítettem szép, látványos példát. Leírom, hogy miként valósítottam ezt meg. Jópár grafikont felhasználtam, és 100 képet kiexportáltam, majd VirtualDub-bal (http://www.virtualdub.org/) videóvá fűztem őket.

A grafikonoknál megadtam, hogy milyen méretben, milyen tartományon legyenek ábrázolva. Ezek gyakran hasznosak, hiszen a prezentáció során fontos, hogy az azonos grafikonok azonos beosztással, mérettel készüljenek, a videónál pedig alapkövetelmény.

Alapeloszlás

Keressünk egy megfelelõ sûrûségfüggvényt. Mivel ilyen púpos lett, a felújítási függvénynél majd lehet látni, hogy a pontok között lévõ távolság hol nagyobb, hol kisebb.

(∫_0^1)-x (x - 1) (1 + Cos[2 π x] + x - 
            0.49) \[DifferentialD]x

f[x] -szel jelölöm a sûrûségfüggvényt. Természetesen normálok, hogy sûrûségfüggvényt kapjak.

f[x_] := 
      If[x < 0 || x > 1, 
        0, -x (x - 
                1) (1 + Cos[2 π x] + x - 
                  0.49)/(0.11767274151216445) // N]

Grafikonon:

Plot[f[x], {x, 0, 1}, ImageSize -> {640, 480}]

Realizáció

Készítsük realizációt! Ha F az eloszlásfüggvény, akkor (F^(-1))[RND] megfelelõ eloszlású lesz. Ilyeneket adok össze, és elmentem a részletösszegeket.
q tartalmazza a véletlen számot 0 és 1 között. r egy lista, amely folyamatosan nő egy taggal, és tartalmazza a pontok helyeit ((Subscript[S, n]) értékeit). Közelítõ megoldást alkalmazok az inverz kiszámolására. Legenerálom az elsõ 30 pontot, és kirajzolom. Az eljárás egy percig is futhat!

r = {0};

Do[q = Random[]; 
      r = Join[r, {r[[-1]] + FindRoot[(∫_0^y)-x (x - 
                            1) (1 + Cos[2 π x] + x - 
                              0.49)/(0.11767274151216445) \[DifferentialD]x == q, {y, 0.2, 0, 1}][[1, 
                2]]}], {30}]

ListPlot[{Drop[r, 1], Table[0.01, {Length[r] - 1}]} // 
        Transpose, {Axes -> {True, False}, 
        PlotStyle -> {Hue[0], PointSize[ .015]}, 
        ImageSize -> {640, 480}, 
        PlotRange -> {{0, 15}, {-1, 1}}, 
        Ticks -> {Range[16] - 1, Automatic}}]

Ugyanezt megcsinálom úgy, hogy az alapeloszlás exponenciális. Így az ábrán Poisson pontfolyamatot kapunk.

r = {0};

Do[q = Random[]; 
      r = Join[r, {r[[-1]] + FindRoot[(∫_0^y) (E^(-x)) \[DifferentialD]x == q, {y, 0.1, 0, 100}][[1, 
                2]]}], {20}]

r

ListPlot[{Drop[r, 1], Table[0.01, {Length[r] - 1}]} // 
        Transpose, {Axes -> {True, False}, 
        PlotStyle -> {Hue[0], PointSize[ .015]}, 
        ImageSize -> {640, 480}, 
        PlotRange -> {{0, 15}, {-1, 1}}, 
        Ticks -> {Range[16] - 1, Automatic}}]

Konvolúció számolása

Az alábbiakban az első 20 konvolúció sûrûségfüggvényeit fogjuk közelítõleg meghatározni. Azaz (Subscript[S, n])eloszlás sûrûségfüggvényét az x pontban s[n, Floor[100*x]] fogja megmondani. Értelemszerûen s[1, . ] alakúak meg fognak egyezni f függvény megfelelõ értékeivel, míg s[n, . ] értékei s[ n-1, . ] értékekbõl és f-bõl származtatom. A 100-as osztó fontos, hiszen minden egész között 100 különbözõ értéket tárolok el. Valójában az összegzésnél közelítõ integrálást végzünk.

A program sokáig fut, ezért íratom ki a számokat, hogy lássam előbb-utóbb véget ér.

s = Table[-1, {j, 20}, {i, 2000}];

For[i = 1, i ≤ 2000, i++, 
      s[[1, i]] = f[i/100]]

For[j = 2, 
      j ≤ 20, j++, 
For[i = 1, 
        i ≤ 2000, i++, 
        s[[j, 
            i]] = (1/100) (∑_(k = 1)^i)(s[[j - 1, 
                      i - k + 1]] s[[1, 
                      k]] // N)
]; 
      Print[j]]

Kirajzolom az első pár konvolúció sűrűségfüggvényét.

Table[
      Plot[s[[j, Ceiling[x*100]]], {x, 
          0, 4}, PlotRange -> {{0, 4}, {0, 2.3}}, 
        ImageSize -> {640, 480}], {j, 1, 4}]

m(t) számolása

Az n. konvolúció (Subscript[F, n]) eloszlásfüggvénye t helyen könnyen számolható, összegezzük az Floor[100*t]-nél kisebb értékekre s[[ n, . ]] -et. m(t) meghatározásához pedig összegezni kell (Subscript[F, n])(t)-ket n-re.

ListPlot[
      Table[{i/100, (∑_(j = 1)^20)(∑_(k = 1)^i) 
                  s[[j, k]]/100}, {i, 
          250}], {PlotJoined -> True, 
        ImageSize -> {640, 480}}]

A várható érték reciproka:

1/(∫_0^1) 
          y ((-y (y - 1) (1 + Cos[2 π y] + y - 
                    0.49`))/0.11767274151216445)` \[DifferentialD]y

ListPlot[
      Table[{i/100, (∑_(j = 1)^20)(∑_(k = 1)^i) 
                    s[[j, k]]/
                  100}, {i, 250}] ⋃ 
        Table[{i/200, 1.7518722630595462 i/200}, {i, 
            500}], {PlotJoined -> False, 
        ImageSize -> {640, 480}}]

A különbség 10-nél:

(∑_(j = 1)^20)(∑_(k = 1)^250) s[[j, 
              k]]/100 - 
      2.5 1.7518722630595462

Az elméletileg számolt asszimtotikus különbség: ((E((X^2)) - 2 (((EX))^2))/(2 (((EX))^2)))

(((∫_0^1) (y^2) ((-y (y - 1) (1 + Cos[2 π y] + y - 0.49`))/0.11767274151216445)` \[DifferentialD]y - 2 ((( (∫_0^1) y ((-y (y - 1) (1 + Cos[2 π y] + y - 0.49`))/0.11767274151216445)` \[DifferentialD]y))^2))/(2 ((( (∫_0^1) y ((-y (y - 1) (1 + Cos[2 π y] + y - 0.49`))/0.11767274151216445)` \[DifferentialD]y))^2)))

Nagyobb értékekre már nem pontos a közelítés.

Következő pont távolságának eloszlása

Felejtsük el egy kicsit az eddigi feladatot! Képzeljük el azt, hogy van egy diszkét eloszlásunk, amely k/100 alakú számokra van koncentrálva, amelyeknek f(k/100)/100 súlyt ad. Ennek a konvolúciójának adatait tartalmazza s. Hozzunk létre egy Z táblázatot. Z[[a, b]] azt mondja meg, hogy "a" idő eltelte után ha tekintjük a következő pont távolságának eloszlását, akkor az eloszlának a sűrűségfüggvénye (a diszkét közelítésben a kérdéses valószínűség) mennyi b egység (vagyis b század) múlva.

Azt a valószínűséget kell kiszámolni, amikor "a" időpont után nem történt esemény egészen (a+b)-ig, majd (a+b)-nél bekövetkezett egy. 1. eset: történt már esetmény "a" előtt: felbontom az eseteket aszerint, hogy mikor volt az utolsó esemény, amelynek idejét jelöljük k-val. Ez történhetett 1-től (a-1) -ig bármikor, tehát utána az (a+b-k+1) időintervallumban 1 esemény történt, előtte pedig bármennyi (ezt számolom a j-vel való szummázásnál). 2. eset: az első esemény a+b-nél következik be.

Z[a_, 
        b_] := ((∑_(k = 1)^(a - 1))(s[[1, 
                      a + b - k + 1]]/
                    100 ((∑_(j = 1)^19) s[[j, 
                          k]]/100)) + 
            s[[1, a + b]]/
              100)*100

A képeket exportálom az export parancsal. A képeknevei "XY"+(sorszám-1)+".bmp" lesz. A többi csak formázás kérdése. Az exportált képek a következő pont eloszlását mutatják, mikozben telik az idő.

Foor[l = 1, l ≤ 200, l++, 
        ListPlot[
          Table[{k/100 + l 2/100 - 2/100, Z[2*l, k]}, {k, 
              100}], {PlotJoined -> False, 
            PlotRange -> {{2 l/100, 1 + 2 l/100}, {0, 
                  3}}, {ImageSize -> {640, 480}}, 
            Ticks -> {(Range[35]/4 - 1/2 // N), Automatic}, 
            AxesOrigin -> {2 l/100 // N, 0}}]];

Foor[l = 1, l ≤ 200, l++, 
        Export["Z" <> ToString[l - 1] <> ".bmp", 
          ListPlot[
            Table[{k/100 + l 2/100 - 2/100, Z[2*l, k]}, {k, 
                100}], {PlotJoined -> False, 
              PlotRange -> {{2 l/100, 1 + 2 l/100}, {0, 
                    3}}, {ImageSize -> {640, 480}}, 
              Ticks -> {(Range[35]/4 - 1/2 // N), Automatic}, 
              AxesOrigin -> {2 l/100 // N, 0}}]]];

Határesetben a következő függvényt kapjuk (ez az elméletből jött ki):

g = ((1 - (∫_0^x)-((y (y - 1) (1 + Cos[2 π y] + y - 0.49`))/0.11767274151216445)` \[DifferentialD]y)/((∫_0^1) y ((-y (y - 1) (1 + Cos[2 π y] + y - 0.49`))/0.11767274151216445)` \[DifferentialD]y))

Plot[
      g, {x, 0, 1}, {PlotRange -> {0, 3}, 
        ImageSize -> {640, 480}}]

Előző pont távolságának eloszlása

Hozzunk létre egy Y táblázatot. Y[[a, b]] azt mondja meg, hogy "a" idő eltelte után ha tekintjük az előző pont távolságának eloszlását, akkor az eloszlának a sűrűségfüggvénye (a diszkét közelítésben a kérdéses valószínűség) mennyi b egységgel (vagyis b századdal) ezelőtt. Természetesen kis "a"-ra b=a esetben nagy értéket kapunk, hiszen pozitív valószínűséggel nem esett még pont a (0,a) intervallumba, de ezzel most nem törődünk.

Azt a valószínűséget kell kiszámolni, amikor "a" időpont előtt nem történt esemény egészen (a-b)-től, majd (a-b)-nél bekövetkezett egy. 1. eset: történt már esetmény "a-b" előtt.

Y[a_, b_] := 
      If[b == a, 0, 
        If[b > a, 
          0, ((∑_(k = 1)^(2000 - b))s[[1, 
                    b + k]]/
                  100) ((∑_(j = 1)^19) s[[j, 
                    a - b]]/100)*100]]

Foor[l = 1, l ≤ 200, l++, 
        ListPlot[
          Table[{k/100 - 1/100, Y[2*l, k]}, {k, 100}], {PlotJoined -> 
              False, PlotRange -> {{0, 1}, {0, 
                  3}}, {ImageSize -> {640, 480}}, 
            Ticks -> {(Range[35]/4 - 1/2 // N), Automatic}, 
            AxesOrigin -> {0, 0}}]];

Foor[l = 1, l ≤ 200, l++, 
        Export["Y" <> ToString[l - 1] <> ".bmp", 
          ListPlot[
            Table[{k/100 - 1/100, Y[2*l, k]}, {k, 100}], {PlotJoined -> 
                False, PlotRange -> {{0, 1}, {0, 
                    3}}, {ImageSize -> {640, 480}}, 
              Ticks -> {(Range[35]/4 - 1/2 // N), Automatic}, 
              AxesOrigin -> {0, 0}}]]];

Intervallum hosszának eloszlása

Hozzunk létre egy T táblázatot. T[[a, b]] azt mondja meg, hogy "a" idő eltelte után ha tekintjük az előző esemény és a következő esemény között lévő intervallum hosszának eloszlását, mennyi az eloszlának a sűrűségfüggvénye (a diszkét közelítésben a kérdéses valószínűség).

Azt a valószínűséget kell kiszámolni, amikor "a" időpont előtt nem történt esemény egészen (a-b)-től, majd (a-b)-nél bekövetkezett egy. 1. eset: történt már esetmény "a-b" előtt.

T[a_, b_] := 
      s[[1, b]]/
          100 ((∑_(k = 1)^b)(If[k > a, 0, 
                If[k == a, 
                  1, ((∑_(j = 1)^19) s[[j, 
                          a - k]]/
                        100)]]))*100

Foor[l = 1, l ≤ 200, l++, 
        ListPlot[
          Table[{k/100 - 1/100, T[2*l, k]}, {k, 100}], {PlotJoined -> 
              False, PlotRange -> {{0, 1}, {0, 
                  5}}, {ImageSize -> {640, 480}}, 
            Ticks -> {(Range[35]/4 - 1/2 // N), Automatic}, 
            AxesOrigin -> {0, 0}}]];

Foor[l = 1, l ≤ 200, l++, 
        Export["T" <> ToString[l - 1] <> ".bmp", 
          ListPlot[
            Table[{k/100 - 1/100, T[2*l, k]}, {k, 100}], {PlotJoined -> 
                False, PlotRange -> {{0, 1}, {0, 
                    5}}, {ImageSize -> {640, 480}}, 
              Ticks -> {(Range[35]/4 - 1/2 // N), Automatic}, 
              AxesOrigin -> {0, 0}}]]];

Az elméleti úton számolt határérték:

h = (((∫_0^x) y (( y (y - 1) (1 + Cos[2 π y] + y - 0.49`))/0.11767274151216445)` \[DifferentialD]y)/((∫_0^1) y ((-y (y - 1) (1 + Cos[2 π y] + y - 0.49`))/0.11767274151216445)` \[DifferentialD]y))

(∫_0^1) x f[x] \[DifferentialD]x // N

Plot[((f[x] x)/0.5708178735894571), {x, 0, 
        1}, {ImageSize -> {640, 480}, 
        PlotRange -> {{0, 1}, {0, 5}}}]