Véletlen partíciók generálása

Programozási Feladat 3

A kínai étterem folyamat

A következő program a kínai étterem folyamatot reprezentálja. A beérkezők számát megadva az egyes asztaloknál ülő emberek számát kapjuk eredményül.

ChinesePartition[hanyan_]:= Module[{},lista={1};
Table[m=RandomChoice[Append[lista,1] -> Range[Length[lista]+1]];
If[ m>Length[lista], AppendTo[lista,1],  lista[[m]]= lista[[m]]+1 ],{hanyan-1}];lista]

ChinesePartition[100]

Ennek megjelenítésére az alábbi program használható:

Graphics[{Array[Circle[{3*(#-1),0}, 0.7]&,Length[lista]],
Flatten[Table[Array[Circle[{(3*(q-1))-Cos[10*#1*Pi/200],Sin[10*#1*Pi/200]}, 0.05]&,lista[[q]]],{q,Length[lista]}]]}]

Adott listához megjeleníthető az érkezési sorrend:

ChinesePartitionuj[hanyan_] := Module[{}, lista = {1};
  Table[m = RandomChoice[Append[lista, 1] -> Range[Length[lista] + 1]];
   If[ m > Length[lista], AppendTo[lista, 1],  lista[[m]] = lista[[m]] + 1 ];lista, {hanyan - 1}]]

listauj=ChinesePartitionuj[100]

Manipulate[Graphics[{Array[Circle[{3*(#-1),0}, 0.7]&,Length[listauj[[i]]]],
Flatten[Table[Array[Circle[{(3*(q-1))-Cos[10*#1*Pi/200],Sin[10*#1*Pi/200]}, 0.05]&,listauj[[i]][[q]]],{q,Length[listauj[[i]]]}]]}],{i,1,99,1}]

Legnagyobb ciklus

A ciklusokat rendezhetjük nagyság szerinti sorrendbe is:

Sort[lista,Greater]

Vizsgálhatjuk a legnagyobb keletkezett ciklus méretét. A RandomPartition paranccsal összehasonlítva lényeges eltérés látható, a ChinesePartition legnagyobb ciklusa majdnem minden esetben nagyobb a RandomPartition által generált legnagyobb ciklus méreténél.

<< Combinatorica`

legnagyobbciklus[hanyan_,hanyszor_]:=Table[Sort[ChinesePartition[hanyan],Greater][[1]],{hanyszor}]

RPlegnagyobbciklus[hany_,hanyszor_]:=Table[Sort[RandomPartition[hany],Greater][[1]],{hanyszor}]

LegNagyobbCiklusAbra[hany_,hanyszor_]:=Show[ListPlot[legnagyobbciklus[hany,hanyszor]],ListPlot[RPlegnagyobbciklus[hany,hanyszor],PlotStyle->{Red}]]

LegNagyobbCiklusAbra[100,100]

Az ábrán látszik, hogy a ChinesePartition legnagyobb ciklusainak kékkel jelölt mérete általában lényegesen nagyobb a beépített programénál.

A ciklusok száma

A (Subscript[K, n]) jelöli az n érkező után elfoglalt asztalok számát, vagyis a ciklusok számát. Erre a következő teljesül: ((Subscript[K, n])/logn)->1, ha n->∞:

A következő ábrán látható, hogy a ciklusok száma logn-hez tart. Érdemes nagyobb lépésekben vizsgálni, hogy nagyobb számokra is fusson.

KNAbra[meddig_,hanyasaval_]:=Show[ListPlot[Transpose[{Range[1,meddig,hanyasaval],Map[Length,Flatten[Table[ChinesePartition[x],{x,1,meddig,hanyasaval},{1}],1]]}]],Plot[Log[x],{x,1,meddig}]]

KNAbra[10000,500]

A következő több próbára számolja ki a ciklusok számának átlagát. Itt is haladhatunk nagyobb lépésekben, hogy  nagyobb számokra is fusson. Jól látszik, hogy minél többször végezzük el a próbát, annál jobban simul a logaritmus görbéjére.

Atlag[mennyi_,hanyszor_,hanyasaval_]:=Module[{},t={};
Array[AppendTo[t,(Sum[Map[Length,Flatten[Table[ChinesePartition[x],{x,1,mennyi,hanyasaval},{hanyszor}],1]][[i]],{i,(#-1)*hanyszor+1,#*hanyszor}])/hanyszor]&,{Floor[mennyi/hanyasaval]}];t]

Atlag[1000,3,100]

AtlagAbra[meddig_,hanyszor_,hanyasaval_]:=Show[ListPlot[Transpose[{Range[1,meddig,hanyasaval],Atlag[meddig,hanyszor,hanyasaval]}]],Plot[Log[x],{x,1,meddig}]]

AtlagAbra[500,5,50]

A következő ábrán látható, hogy ((Subscript[K, n])/logn)->1, ha n->∞:

Hatar[meddig_]:=Show[ListPlot[Table[(Length[ChinesePartition[x]])/(Log[x]),{x,2,meddig}]],Plot[(x^0),{x,1,meddig},PlotStyle -> {Red,Thick}]]

Hatar[1000]

Jól látszik, hogy a ciklusok számában is jelentős eltérés van a RandomPartition és a ChinesePartition között.

RPKNAbra[meddig_,hanyasaval_]:=Show[ListPlot[Transpose[{Range[1,meddig,hanyasaval],Map[Length,Flatten[Table[RandomPartition[x],{x,1,meddig,hanyasaval},{1}],1]]}],PlotStyle->{Red}],ListPlot[Transpose[{Range[1,meddig,hanyasaval],Map[Length,Flatten[Table[ChinesePartition[x],{x,1,meddig,hanyasaval},{1}],1]]}]]]

RPKNAbra[1000,10]

Az ábrán kékkel jelölt ChinesePartition ciklusainak száma lényegesen kisebb a pirossal jelölt RandomPartitionénál, valamint láthatóan a számok növelésével nem nő egyenes arányban a ciklusok száma.

A ciklusok hossza

Itt a ciklusok hoszzát vizsgálom. A program több próba esetén összeszámolja, hogy milyen hosszú ciklusból hány keletkezett.

ciklus[hanyan_,hanyszor_]:=Module[{},listauj=Sort[Flatten[Table[ChinesePartition[hanyan],{hanyszor}],1],Greater];Print[listauj];
k={};
Table[
If[
MemberQ[listauj,i],AppendTo[k,{i,Count[listauj,i]}],Cases[Range[hanyan],Except[i]]
],{i,1,hanyan}];k]

ciklus[1000,10]

Az 1 hosszú ciklusok száma

A keletkezett permutációknál érdemes lehet az 1 hosszú ciklusokat számolni. Ezek várható értéke a ChinesePartition esetében 1.
Ebben is szépen látszik az eltérés a RandomPartitionhoz képest, ahol általában több 1 hosszú ciklus keletkezik.

ciklus1[hanyan_,hanyszor_]:=Module[{},
h={};
Table[egyekszama=Count[ChinesePartition[hanyan],1];
AppendTo[h,egyekszama],{i,1,hanyszor}];h]

ciklus1[100,10]

RPciklus1[hanyan_,hanyszor_]:=Module[{},
h={};
Table[egyekszama=Count[RandomPartition[hanyan],1];
AppendTo[h,egyekszama],{i,1,hanyszor}];h]

RPciklus1[100,10]

Abra1ciklus[hanyan_,hanyszor_]:=Show[ListPlot[ciklus1[hanyan,hanyszor]],ListPlot[RPciklus1[hanyan,hanyszor],PlotStyle->{Red}]]

Abra1ciklus[100,500]

A következő esetben a legfeljebb "hanyhosszu" hosszúságú ciklusokat számoljuk meg. A jegyzet szerinti jelölést alkalmaza: (Subscript[C, j]^n), ahol n a "hanyan", j pedig a "hanyhosszu", vagyis n számra a j hosszú ciklusokat számoljuk meg, több próba esetén. A ciklusval eredménye a következő: a listán belüli i. lista az i méretű ciklusokat számolja. Az i. listán belül pedig azt kapjuk, hogy a j. próbára hány darab i hosszú ciklust kaptunk.

ciklusval[hanyan_,hanyszor_,hanyhosszu_]:=Transpose[Module[{},
h=Table[{},{hanyhosszu}];
Table[
Table[szama=Count[ChinesePartition[hanyan],j];
AppendTo[h[[j]],szama],{i,1,hanyszor}],{j,1,hanyhosszu}];h]]

ciklusval[20,10,2]

A j hosszú ciklusok száma független 1/j várható értékű valószínűségi változók.

VarhatoErtek[hanyan_,hanyszor_,hanyhosszu_]:=Table[(Total[Transpose[ciklusval[hanyan,hanyszor,hanyhosszu]][[j]]])/hanyszor,{j,1,hanyhosszu}]

N[VarhatoErtek[10,5000,4],3]

Most megszámoljuk, hogy a több próba során hányszor hány darab adott hosszúságú ciklust talált. (Pl. 10 próbából 2-szer volt 1 darab 1-hosszú, 3-szor 2 db 1 hosszú ciklus, stb...). A jegyzet szerinti jelöléssel a kapott listán belüli i. lista j. eleme azt adja meg, hogy "hanyszor" próba után (Subscript[C, hanyhosszu]^hanyan) hányszor vett fel j értéket.

ciklusszam[hanyan_,hanyszor_,hanyhosszu_]:=Module[{},
h=Table[{},{hanyhosszu}];
Table[
Table[szama=Count[ChinesePartition[hanyan],j];
AppendTo[h[[j]],szama],{i,1,hanyszor}],{j,1,hanyhosszu}];
o=Table[Table[Count[h[[q]],w],{w,1,hanyan}],{q,1,hanyhosszu}];o]

ciklusszam[10,10,2]

Véletlen ülésrend adott valószínűségekhez

Véletlen ültetési rend parciálisan felcserélhető partícióhoz. Adott (Subscript[P, i]) valószínűségekhez generálunk ülésrendet. A (Subscript[P, i]) valószínűségeknek teljesíteiük kell a következőket: 0 ≤ (Subscript[P, i]) ∀i-re, illetve (∑_(i=1)^k)(Subscript[P, i]) < 1. Afeltételeknek megfelelő valószínűségeket generálunk (a lényeg az, hogy bármely ilyen adott valószínűségekre működik).

val=Table[RandomReal[],{100}];
val2=Map[#/(Sum[val[[i]],{i,Length[val]}]+1)&,val]

Ehhez pedig a következő módon adhatjuk meg az ülésrendet:

adottsulyokhoz[hanyan_]:= Module[{},lista={1};
Table[m = RandomChoice[ Append[val2[[1;;Length[lista]]],1-(Sum[val2[[p]],{p,Length[lista]}])]-> Range[Length[lista] + 1]];
If[ m>Length[lista], AppendTo[lista,1],  lista[[m]]= lista[[m]]+1 ],{hanyan-1}];lista]

adottsulyokhoz[96]

Graphics[{Array[Circle[{3*(#-1),0}, 0.7]&,Length[lista]],
Flatten[Table[Array[Circle[{(3*(q-1))-Cos[10*#1*Pi/200],Sin[10*#1*Pi/200]}, 0.05]&,lista[[q]]],{q,Length[lista]}]]}]

A kétparaméteres modell

A kétparaméteres modell esetén az előző modellt módosítjuk. Az új érkezők az (α,θ) értékek függvénye alapján megválasztott valószínűséggel foglal helyet.

ketparameteres[alfa_,theta_,hanyan_]:= Module[{a=alfa,t=theta,h = hanyan} ,lista={1};
Table[m=RandomChoice[lista2=Append[lista,(t + a*(Length[lista]+1))];
Function[x,(x-a)/(Sum[lista[[i]],{i,1,Length[lista]}]+ t)] /@lista2 -> Range[Length[lista]+1]];
If[ m>Length[lista], AppendTo[lista,1],  lista[[m]]= lista[[m]]+1 ],{h-1}];lista]

ketparameteres[0,8,100]

Az így kapott listát szintén ábrázolhatjuk:

Graphics[{Array[Circle[{10*(#-1),0}, 2]&,Length[lista]],Array[Inset[(((Subscript[n, #])-α)/(n+ θ)),{10*(#-1),0}]&,Length[lista]-1],Inset[((θ+ kα)/(n+θ)),{10*(Length[lista]-1),0}],
Flatten[Table[Array[Circle[{(10*(q-1))-3*Cos[10*#1*Pi/200],3*Sin[10*#1*Pi/200]}, 0.05]&,lista[[q]]],{q,Length[lista]}]]}]

Ekkor az α és θ paraméterek függvényében a következőket figyelhetjük meg:
Ha α < 0, akkor (Subscript[K, n]) θlogn szerint növekszik.

ParAbra[hanyan_,hanyasaval_,alpha_,theta_]:=ListPlot[{Table[Length[ketparameteres[alpha,theta,x]],{x,1,hanyan,hanyasaval}],Table[2*Log[x],{x,1,hanyan,hanyasaval}]}]

ParAbra[1000,20,-0.001,2]

Prímfelbontások

Itt egy általunk kívánt határig véletlen választott szám típusát keressük.

randomtipus[meddig_]:=Module[{},
m=RandomInteger[meddig];
Print[m];
primek=FactorInteger[m];
legnagyobb= First[Last[primek]];
c=Table[0,{PrimePi[legnagyobb]}];
van=Flatten[primek,1];
csak=Delete[van, Table[{2*u},{u,1,(Length[van])/2}]];
Table[If[MemberQ[csak,Prime[i]],c[[i]] = primek[[Flatten[Position[csak,Prime[i]],1][[1]]]][[2]]],{i,1,Flatten[Position[Table[Prime[n],{n,1,m}],legnagyobb],1][[1]]}];c]

randomtipus[10000]

A következőkben egy adott szám c típusát keressük:

tipus[szam_]:= Module[{},
primek=FactorInteger[szam];
legnagyobb= First[Last[primek]];
c=Table[0,{PrimePi[legnagyobb]}];
van=Flatten[primek,1];
csak=Delete[van, Table[{2*u},{u,1,(Length[van])/2}]];
Table[If[MemberQ[csak,Prime[i]],c[[i]] = primek[[Flatten[Position[csak,Prime[i]],1][[1]]]][[2]]],{i,1,Flatten[Position[Table[Prime[n],{n,1,m}],legnagyobb],1][[1]]}];c]

tipus[4232345]

primekszama[szam_]:=Total[tipus[szam]]

A következő program nagy számokra nézi meg a prímtényezők számát (multiplicitással). Véletlen számonként nézi, mert ha adott számokkal haladunk, akkor az erősen befolyásolja a prímtényezők számát (pl. 1000-től haladva 10000-ével nem ad túl általános képet). Hibát ír ki, de kiszámolja.

PrimTenyezok[meddig_,hanyasaval_]:=Module[{},l=Table[Sum[Table[RandomInteger[hanyasaval],{meddig/hanyasaval}][[i]],{i,1,j}],{j,1,meddig/hanyasaval}];Table[primekszama[l[[i]]],{i,1,meddig/hanyasaval}]
]

plist=PrimTenyezok[4298364,213222]

A prímtényezők multiplicitása

A multipl[] kiszámolja, hogy egy adott számig véletlen választott számban a keresett prím milyen kitevőn szerepel.

multipl[meddig_,prim_]:=Module[{}, e=Select[FactorInteger[RandomInteger[meddig]],#[[1]]==prim&];If[e=={},0,e[[1,2]]]]

multipl[10000,2]

A következő pedig egy határig véletlen választott számok prímfelbontásában keresi a kérdéses prím kitevőjét. Ez nyilván a szám prím alapú logaritmusához tart.

Table[multipl[1000000,2],{1000}]

multiplhiszt[meddig_,prim_,hany_]:=Module[{},Show[ListPlot[a=Table[multipl[meddig,prim],{hany}];Table[Count[a,i],{i,0,Log[prim,meddig]}]],Plot[hany*(1-1/prim)(((1/prim))^(x-1)),{x,0,Log[prim,meddig]}]]]

multiplhiszt[100000,2,1000]

A témához kapcsolódó további számítások

A RandomPartitionnal szemben jól látszik, hogy a futásidőben is jelentős különbségek vannak. Míg a  ChinesePartition[] akár (10^6)-os nagyságrendre is használható időben fut, a RandomPartition[] már a 1000-es nagyságrendet is rendkívül lassan számolja.

AbsoluteTiming[RandomPartition[2000]]

AbsoluteTiming[ChinesePartition[100000]]