Kombinatoryka i prawdopodobieństwo

Generowanie podzbiorów

Kombinacje, permutacje, wariacje, podzbiory vs Tuples, Permutations, Subsets

Mathematica posiada zestaw komend pozwalających bardzo łatwo generować rozmaite podzbiory używane w rachunku prawdopodobieństwa. Niestety, związek tych komend z terminologią używaną w nauczaniu matematyki w Polsce jest dosyć luźny. Po prostu nie ma odpowiedników naszych określeń typu ,,wariacje bez powtórzeń’’ czy ,,kombinacje’’. Mathematica posiada komendy Tuples, Permutations (permutacje) oraz Subsets (podzbiory).

Dla przykładu rozważmy zbiór 4-elementowy

zbior={"amelinium", "tego", "nie", "pomalujesz"}

Polecenie Tuples generuje wszyskie możliwe sposoby wybrania k=3 elementów powyższego zbioru. Możemy sobie to wyobrazić jako losowanie k=3 razy z urny zawierającej 4 karteczki ale zwracając karteczkę za każdym razem. Po naszemu nazywamy to wariacją z powtórzeniami:

Tuples[zbior,3]

Polecenie Permutations generuje permutacje, czyli wszystkie możliwe przestawienia:

Permutations[zbior]

Permutations ma drugi parametr. Jeżeli podamy do, wygenerowane zostaną wszystkie permutacje podzbiorów co najwyżej k=3 elementowych:

Permutations[zbior,3]

Jeżeli dodatkowo podamy parametr k w nawiasach wąsatych, to wygenerowane zostaną wszystkie permutacje podzbiorów dokładnie k=3 elementowych, W nomenklaturze polskiej są to wariacje bez powtórzeń:

Permutations[zbior,{3}]

Polecenie Subsets generuje podzbiory. Istotna różnica polage na tym, że teraz nie rozróżniamy elementów ani ich kolejności. Wszystkie podzbiory:

Subsets[zbior]

Wszystkie możliwe podzbiory zawierające conajwyżej k=3 elementy:

Subsets[zbior,3]

Wszystkie podzbiory zawierające dokładnie k=3 elementy. Po polsku są to kombinacje bez powtórzeń:

Subsets[zbior,{3}]

Poniżej przydatna tabelka pokazująca jaką komendę należy użyć do wygenerowania znanych ze szkolnego kursu prawdopodobieństwa zbiorów:

{{"Nazwa polska", "Polecenie\ngenerujące", "Liczba"},{},
{"Wariacje z powtórzeniami", "Tuples[zbior,k]", "n^k"},{},
{"Wariacje bez powtórzeń","Permutations[zbior,{k}]","n!/(n-k)!"},{},
{"Permutacje", "Permutations[zbior]", "n!"},
{"Kombinacje","Subsets[zbior, {k}]", "(Binomial[n, k]) = Binomial[n,k]"}
}//TableForm

Przykład: wygenerować wszystkie możliwe n-krotne złożenia funkcji f1[x]=x+1 oraz f2[x]=1/x i obliczyć wynik złożenia dla x=1

Zastosowanie metod kombinatorycznych na szerszą skalę jest możliwe w praktyce tylko przy użyciu komputera, stąd ich mała popularność w nauczaniu tradycyjnym. Można oczywiście wykonać poniższe rachunki ręcznie, ale będzie to trwało nieproporcjonalnie długo w porównaniu do osiągniętych korzyści. Zastosowanie komputera odwraca tą sytuację. Otóż często kombinatoryczne  rozwiązanie okazuje się być tym wymagającym najmniejszej wiedzy i wysiłku.

Generowanie wszytkich kombinacji na raz

Definiujemy 2 funkcje:

f1=#+1& (* funkja dodająca 1 *)

f2=1/#& (* funkcja obliczająca odwrotność *)

f1[3]

f1[2]

f2[3]

Pytamy się, jaki może być rezultat składania tych funkcji:

f1@f1@f1@f1@x

f2@f2@f2@f2@f2@x

f1@f2@f1@x

Spróbujmy nieco zautomatyzować powyższe czynności. Zaczynamy od losowego wybrania powiedzmy 12 funkcji. Służy do tego RandomChoice (ang: losowy wybór ):

RandomChoice[{f1,f2},12]

Aby je złożyć użyjemy polecenia Composition  (ang: złożenie). Polecenie @@ przekazuje wszystkie elementy listy jako argumenty:

Composition@@%

Teraz możemy obliczyć wynik złożenia dla argumentu x:

%[x]

Przystępujemy do wygenerowania wszystkich możliwości. Zaczynamy od wszystkich możliwych złożeń 6 funkcji:

n=6

funs={f1,f2}

Jeżeli ilość funkcji i liczba złożeń nie jest duża, możemy wygenerować wszystkie kombinacje na raz:

l=Tuples[funs,6]

(Composition@@l[[13]])[x]

(Composition@@l[[13]])[1]

Generujemy wszystkie wyniki:

z=Table[(Composition@@l[[i]])[1],{i,1,Length[l]}]

Jest ich 2^n:

Length[l]

2^n

Niektóre wyniki powtarzają się. Używamy DeleteDuplicates (ang: usuń duplikaty) :

zU=z//DeleteDuplicates

Otrzymaliśmy 33 różne dodatnie liczby wymierne:

Length[zU]

Problem z powyższą metodą ujawnia się dla dużych n. Ilość ,,kombinacji’’ w powyższym przykładzie
to (2^n), natomiast 1 GB pamięci to 2^30 bajtów. Struktury tworzone przez program Mathematica są dosyć ,,pamięciożerne’’:

ByteCount[l]

Zajęta pamięć to znacznie więcej niż 2^n bajtów (2^16 vs 2^6):

Log[2,%]//N

ByteCount[l]/2^n//N

Generowanie kobinacji jedna po drugiej

Niestety, Mathematica nie posiada wbudowanej funkcji generującej Tuples po jednej. Aby to zrobić, zauważmy, że generowanie Tuples można jednoznacznie przypisać liczbie w systemie n-konwym. Np. poniższa komenda jest równoważna wypisaniu wszystkich 2-cyfrowych liczb w systemie 10-tkowym.

Tuples[{0,1,2,3,4,5,6,7,8,9},2]//TableForm

Tuples[{0,1},4]//TableForm

Idea jest więc  następująca: generujemy w pętli liczby naturalne (co jest bardzo proste), zamieniamy je na liczbą w systemie np: dwójkowym, rozbijamy na listę pojedynczych cyfr, a następnie zamieniamy cyfry na elementy zbioru.

Jakie złożenie funkcji odpowiada np: liczbie k=123?

k=123

Listę cyfr w systemie dwójkowym generuje polecenie IntegerDigits:

IntegerDigits[k,2]

Teraz wystarczy zamienić 0 na funkcję f1, a 1 na funkcję f2:

%/.{0->f1,1->f2}

Composition@@%

%[1]

Teraz możemy napisać pętlę, która powtórzy powyższe czynności dla kolejnych wartości k:

Clear[k]

For[k=0,k<128,k++,
x=IntegerDigits[k,2];
f=x/.{0->f1,1->f2};
w=(Composition@@f)[1];
Print[w]
]

Oczywiście generowanie takich złożeń musi czemuś służyć. Np: spróbujemy ustalić, czy wśród powyższych złożeń są wszystkie liczby wymierne o liczniku i mianowniku nie większym od 10:

?Rest

listaR=Flatten[Table[m/n,{m,0,10},{n,1,10}],1]//DeleteDuplicates//Sort//Rest

?Break

Poniższa pętla sprawdza czy wynik złożenia jest na liście. Jeżeli tak, to usuwa go z listy. Pętla jest przerywana, gdy lista stanie się pusta {}. Jak widać, pętla poniżej jest ,,gotowa’’ na obliczenia trwające wiele godzin czy dni, stąd wypisywanie daty i godziny.

lista={};

listaOLD=listaR;

listaNEW=listaR;

For[ii=1,ii<10^9,ii++,
listaOLD=listaNEW;
w=(Composition@@(IntegerDigits[ii,2]/.{0->f1,1->f2}))[1];
notpresent=Position[lista,Expand[w]]=={};
If[notpresent,AppendTo[lista,w]];
listaNEW=DeleteCases[listaOLD,w];
(*If[notpresent,Print[Date[]//DateString,"\t",ii,"\t",listaNEW]]*)
If[listaNEW≠listaOLD,Print[Date[]//DateString,"\t",ii,"\t",listaNEW]];
If[listaNEW=={},Break[]]
]

Rzut kostką

Poniższe przykłady ilustrują ciekawe możliwości programu Mathematica. Wykorzystamy fakt, iż symbolem nie musi być literka czy nawet słowo. Może być nim plik graficzny. Rozważamy rzuty kostką. Wgrywamy zdjęcia wszystkich 6 ścianek:

kostki=
{
Import["C:\\temp\\kostki\\jeden_86x86.png"],
Import["C:\\temp\\kostki\\dwa_86x86.png"],
Import["C:\\temp\\kostki\\trzy_86x86.png"],
Import["C:\\temp\\kostki\\cztery_86x86.png"],
Import["C:\\temp\\kostki\\piec_86x86.png"],
Import["C:\\temp\\kostki\\szesc_86x86.png"]
}

dir="https://ribes.if.uj.edu.pl/alsymb/kostki/"

{
Import[dir<>"jeden_86x86.png"],
Import[dir<>"dwa_86x86.png"],
Import[dir<>"trzy_86x86.png"],
Import[dir<>"cztery_86x86.png"],
Import[dir<>"piec_86x86.png"],
Import[dir<>"szesc_86x86.png"]
}

Wynikiem dwukrotnego rzutku kostką musi być jedna z pomizszych sekwencji:

rzuty2=Tuples[kostki,2]

Ilość wszystkich możliwych zdarzeń to:

rzuty2//Length

Ile razy pojawiają się 2 szóśtki? Odpowiedź jest trywialna, ale bez problemu możemy znaleźć ją w poniższy sposób:

Select[rzuty2,#=={,}&]

Powyżej, użyłem polecenia Select (wybierz) do wybrania tych elementów zbioru, które są parą szóstek. Jest to bardzo prosty przykład. Spróbujmy więc czegoś odrobinę trudniejszego. Ile razy pojawiają się w trzykrotnym rzucie conajmniej 2 szóstki?

rzut3=Tuples[kostki,3];

Moglibyśmy je po prostu policzyć, ale zróbmy to porządnie. Na początek przeanalizujemy jeden rzut:

rzut3[[123]]

Chcielibyśmy ustalić (przy pomocy komend programu Mathematica !) czy powyżej są 2 lub więcej szóstek. Można to pytanie sformułować inaczej: czy wśród podzbiorów jest zbiór 2 szóstek?

Subsets[%]

Jeżeli jest, to możemy ustalić jego pozycję:

Position[{,},%]

W powyższym przykładzie nie ma 2 szóstek, więc wynikiem jest zbiór pusty {}. Wybierzemy więc te elementy zbioru rzut3, dla których powyższy zbiór jest niepusty, czyli zawiera 2 szóstki:

Select[rzut3,Position[Subsets[#],{,}]≠{}&]

Teraz wystarczy policzyć znalezione przypadki:

Length[%]

Wszystkich możliwości było:

Length[rzut3]

Prawdopodobieństwo wynosi więc:

P=%%/%

N[p]*100

Student, który sumiennie przerobił zadania z prawdopodobieństwa, potrafi obliczyć powyższe, np: korzystając ze schematu Bernouliego. Powodzenia:

(1/6)^2*(5/6)*Binomial[3,2]+(1/6)^3

p=1/6

q=1-p

n=3

Sum[p^k q^(n-k) Binomial[n,k],{k,2,3}]

Rzut kostką v2

Poprzedni przykład stanowi raczej ciekawostkę, ale w przyszlości: kto wie? W każdym razie dziś raczej operujemy symbolami matematycznymi w tradycyjnym rozumieniu:

kostki={1,2,3,4,5,6}

rzut3=Tuples[kostki,3]

Select[rzut3,Position[Subsets[#],{6,6}]≠{}&]

Length[%]

Teraz możemy bez obawy o zawieszenie komputera spróbować ,,większych’’ pytań. Na przykład, jakie jest prawdopodobieństwo, że suma oczek przy 6-krotnym rzucie kostką jest liczbą pierwszą? Zaczynamy od wygenerowania przestrzeni zdarzeń elementarnych:

rzut6=Tuples[kostki,6];

Poniżej schemat wyszukiwania zdarzeń sprzyjających. Dla każdego zdarzenia, np:

rzut6[[123]]

obliczamy sumę ,,oczek’’:

Total[%]

i sprawdzamy czy jest ona liczbą pierwszą:

PrimeQ[%]

Select[rzut6,PrimeQ[Total[#]]&]

Length[%]

Length[rzut6]

%%/%

N[%]

Jak widać, prawdopodobieństwo jest zaskakująco duże. Powyższą metodą można rozwiązać sporo ,,nietypowych’’ zadań. Np. jakie jest prawdopodobieństwo, że przy 6 rzutach żadna ścianka się nie powtórzy?

Length@Permutations[kostki]/Length[rzut6]

Obowiązuje tu jednak prawo zachowania trudności: musimy napisać funkcję selekcjonującą zdarzenia, działającą poprawnie,  a nie jak selekcjoner reprezentacji w piłce nożnej.

Rzut kostką: metoda Monte Carlo

Prawdopodobieństwo można oszacować, a czasem także wyznaczyć ściśle, naśladując w komputerze proces np: rzucania kostkami. Powtarzając go miliony razy, możemy w miarę dokładnie obliczyć prawdopodobieństwo. Trwa to bardzo krótko, a co więcej jest to metoda, która idealnie nadaje się dla komputerów wieloprocesorowych oraz systemów rozproszonych.

dir="https://ribes.if.uj.edu.pl/alsymb/kostki/"

{
Import[dir<>"jeden_86x86.png"],
Import[dir<>"dwa_86x86.png"],
Import[dir<>"trzy_86x86.png"],
Import[dir<>"cztery_86x86.png"],
Import[dir<>"piec_86x86.png"],
Import[dir<>"szesc_86x86.png"]
}

kostki (* w razie potrzeby wczytać ponownie wcześniej *)

Poleceniem, które ułatwia prowadzenie symulacji losowania jest RandomChoice:

RandomChoice[kostki,3]

Dalej zrezygnuję z obrazków, na rzecz liczb:

kostki={1,2,3,4,5,6}

Ustalmy jakąś liczbę losowań, powiedzmy 128:

num=2^7

Spróbujemy ustalić, jaka jest szansa, że pojawią się 2 szóstki

Table[RandomChoice[kostki,2],{num}]

Select[%,#=={6,6}&]

Length[%]/num//N

1/36//N

Tak oszacowane prawdopodobieństwo jest bardzo przybliżone, co gorsza od czasu do czasu trafiamy na 0. Potrzebjemy więcej losowań. Szybko dojdziemy do momentu, w którym wygenerowane na raz zbiory nie zmieszczą się  w pamięci komputera. Dlatego będziemy teraz losować je po jednym zdarzeniu. Za każdym razem, gdy RandomChoice wypluje 2 szóstki, powiększamy licznik o 1 komendą i++. Teraz możemy użyć naprawdę dużych liczb:

num=2^23

licznik=0

For[i=0,i<num,i++,
If[RandomChoice[kostki,2]=={6,6},licznik++]
]

Prawdopodobieństwo, to po prostu stosunek licznika do całkowitej ilości losowań:

licznik/num//N

Rationalize[%,0.01]

1/36//N