Graafin erillisten kaariparien laskeminen

Olkoon $G=(V, E)$ yksinkertainen suuntaamaton, äärellinen graafi. Haluamme laskea kuinka monella tapaa $E$ :stä voidaan valita kaksi kaarta $e_1$ ja $e_2$ siten että näiden päätepisteiden joukot $\{v_{11}, v_{12}\}$ ja $\{v_{21}, v_{22}\}$ ovat pistevieraat. (Suuntaamaton kaarihan itse asiassa on päätepisteidensä joukko, joten asian voisi ilmaista myös niin, että kaaret ovat pistevieraat.)

Halutunlaisten kaariparien laskeminen voidaan tehdä käymällä läpi kaikki kaaret ja summaamalla kullekin kaarelle sen pariksi sopivien kaarien lukumäärä. Tämä lukumäärä kaarelle $e=\{v_1, v_2\}$ on $|E|-\deg(v_1) - \deg(v_2) + 1$ .

Summan sieventämiseksi voidaan yhdistää saman lukumäärän antavat kaaret. Tätä varten muodostetaan graafista polynomi $p_E(x)$ asettamalla jokaista kaarta $e = \{v_1, v_2\} \in E$ vastaamaan termi $x^{\deg(v_1)+\deg(v_2)-1}$ , jossa eksponentti siis laskee kaaren päätepisteisiin liittyvien kaarien lukumäärän. Nämä summataan yhteen. Esimerkki:

Merkitään $p_E(x) = a_m x^m + \dots +a_1 x$ (huom. vakiotermi on $0$ , sillä jokaiseen kaareen liittyy ainakin yksi kaari nimittäin se itse). Nyt haluttu lukumäärä saadaan kätevästi:

$\begin{aligned} \sum_{k=1}^{m} a_k (|E|-k). \end{aligned}$

Ylläolevalle esimerkille saadaan ( $|E|=9$ )

$\begin{aligned} 4(9-5) + 4(9-4) + 1(9-3) = 42. \end{aligned}$

Huomiota

Tämä erillisten kaariparien lukumäärä on itse asiassa verkon ”2-paritusten” lukumäärä https://en.wikipedia.org/wiki/Matching_(graph_theory)
Olisi voitu määritellä myös polynomi $q_E$ , asettamalla jokaiseen kaareen suoraan termi $x^{|E|-\deg(v_1)-\deg(v_2)+1}$ eli siihen liittymättömien kaarien lukumäärä. Tällöin haluttu summa saa helposti ilmaistavan muodon: $q_E' (1)$ .
Ja koska yleensä vielä halutaan että derivaatta pitää ottaa nollassa, olisi voitu tehdä siirto $x\mapsto x+1$ . Tällöin esimerkin polynomi olisi ollut $\begin{aligned} q(x) = 4(x+1)^{9-5} + 4(x+1)^{9-4} + (x+1)^{9-3} = x^6 + 10 x^5 + 39 x^4 + 76 x^3 + 79 x^2 + 42 x + 9 \end{aligned}$ ja kuten huomataan, tämän derivaatta nollassa eli ensimmäisen asteen kerroin on $42$ .

Vangin pasianssin voittotodennäköisyys

Vangin pasianssi on korttipeli, jossa tavallisesta 52 kortin pakasta valitaan ensin pöydälle 13 satunnaista korttia. Näistä samanarvoiset yhdistellään pinoiksi. Tämän jälkeen pakasta aletaan nostamaan kortteja, joiden avulla pöydän kortteja voi poistaa. Tämä tehdään nostamalla kerralla kolme korttia, joista vain viimeinen otetaan huomioon. Jos pöydällä on tämän kortin osoittaman arvon pino, se voidaan poistaa. Näin jatketaan, kunnes pakka on tyhjä (tai pöytä tyhjä). Peli menee läpi, mikäli pöytä saadaan tyhjäksi. Kysymys kuuluukin: Mikä on pelin läpimenemisen todennäköisyys?

Esimerkki pelitilanteesta, jossa pöydälle valittu kortit ja ensimmäinen kolmikko pakasta nostettu. Tässä katsotaan korttia risti neljä. Koska pöydällä ei ole nelosia, niin pöydältä ei poisteta mitään.

Kokeile peliä täällä

Voimme olettaa, että jälkimmäisessä (pakan 3 kerrallaan nostelussa) itse asiassa vain valitaan 13 satunnaista korttia. (Koska jokainen joka kolmas kortti on yhtä satunnainen kuin 13 päällimmäistäkin; sillä oletuksella että pakan alkuperäinen järjestys on satunnainen.) Läpimeno on yhtäpitävää sen kanssa, että jälkimmäinen valinta sisältää ainakin yhden kerran jokaisen pöydälle tulleen arvon.

Laskentaa hankaloittaa se, että pöydälle valitut 13 korttia vaikuttavat jäljellä olevan pakan valintojen todennäköisyyksiin. Tämän vuoksi meidän täytyy huomioida minkä tyyppinen valinta pöydälle on tullut. Määritellään tämä ”tyyppi” tarkasti ottaen pöydän pinojen kokojen (laskevasti) järjestettynä vektorina. Ylläolevassa esimerkkikuvassa on tullut tyyppi $(3, 2, 1, 1, 1, 1, 1, 1, 1, 1)$ : 3 vitosta, 2 rouvaa, 1 kasi, 1 ässä, …, 1 kutonen. Tyypissä ei siis huomioida mitä kortteja on tullut, ainoastaan näiden lukumäärät. Tämä rajoittaakin mahdollisten tyyppien määrä huomattavasti: ne ovat itseasiassa luvun 13 ositukset osiin, joista jokainen on korkeintaan 4 (koska maita on 4, niin suurin mahdollinen pinon koko on 4). Ks. OEIS: http://oeis.org/A001400. Näitä on 39 kappaletta ja ne voi generoida seuraavalla Sage-koodilla:

#list of partitions of n to at most maxPartsN parts that are at most maxSize
def getPartitions(n, maxPartsN, maxSize):
    if n==0: return [[]]
    canMakeMax = maxPartsN*maxSize
    if n>canMakeMax: return []
    if n==canMakeMax: return [[maxSize]*maxPartsN]
    ret = []
    for x in xrange(min(n, maxSize), -1, -1):
        ret += [[x] + sP for sP in getPartitions(n-x, maxPartsN-1, x)]
    return ret

def getAllPossibleTypes(boardN, suits, vals):
    return getPartitions(boardN, min(boardN, vals), suits)

a = getAllPossibleTypes(13, 4, 13)
print "{} types:\n{}".format(len(a), a)

Seuraavaksi meidän täytyy laskea mikä on kullekin tyypille $t=(t_1, t_2, \dots, t_r)$ todennäköisyys tulla pöydälle ja sitten mikä valita jäljellä olevasta pakasta 13 korttia siten, että jokaista tyypin korttia tulee ainakin yksi. Tässä onkin kätevää, että jälkimmäisessä valinnassa pakasta on poistettu tyypin mukaan juuri niitä kortteja, joita halutaan tulevan, joten sillä mitkä kortit ne olivat, ei ole väliä, vaan ainoastaan niiden tyypillä $t$ . Merkitään näitä kysyttyjä todennäköisyyksiä

$\begin{aligned} A(t) =& \mathbb{P}(\text{saadaan 1. valinnassa tyyppi } t) \\ B(t) =& \mathbb{P}(\text{kun pakasta poistettu tyypin } t \text{ mukaan, saadaan 2. valinnassa jokaista poistettua arvoa ainakin } 1) \end{aligned}$

Lasketaan ensin $A(t)$ . Kuinka moni pakan 13-kombinaatioista on tyyppiä $t =(t_1, t_2, \dots, t_r)$ ? Lasketaan tulo

$\begin{aligned} \prod_{j=1}^r (13-j){{4}\choose{t_j}} \end{aligned}$

tässä valitaan jokaiselle $t$ :n komponentille arvo ja kerrotaan sillä kuinka monella tavalla tuon arvoiset kortit voidaan valita maiden joukosta. Huomaa: koska tyyppi on järjestetty vektori, arvojen valintoja ei jaeta $r!$ :llä, mutta tulos täytyy jakaa sillä kuinka monella tavalla tyypin luvut voivat permutoida tyypin muuttumatta. Se on tyypin ”itsensä tyypin” yli laskettujen kertomien tulo $\prod_{j} {y_j}!$ , kun tyypin tyyppi on $(y_j)_j$ .

Sage-koodi näiden laskemiseen:

import fractions
from collections import Counter

def aProb(t, suits, vals):
    num = 1
    for i in xrange(len(t)):
        num *= (vals-i)*binomial(suits, t[i])
    num = num/reduce(lambda x,y: x*factorial(y), Counter(t).values(), 1)
    denom = binomial(suits*vals, sum(t))
    return fractions.Fraction(int(num), int(denom))

Sitten lukujen $B(t)$ laskemiseen. Olkoon taas tyyppi $t=(t_1, t_2, \dots, t_r)$ . Nyt jäljellä oleva pakka on multijoukko $X = \{(4-t_j)\times x_j\}$ , missä $x_j$ on tyypin komponenttia $j$ vastaavan kortin arvo ja $t_j=0, \text{ kun } j>r$ . Tällaisen valinnan tekemisen todennäköisyyttä, että tiettyjä arvoja saadaan ainakin yksi, käsittelimme viime kerralla: https://membolicsythod.home.blog/2019/05/19/multijoukon-valintojen-todennakoisyyksia/ . Nyt vähintään kerran valituksi vaadittavat alkiot ovat $x_1, x_2, \dots, x_r$ . Kuten sanottu, sillä ei ole väliä mitkä nuo arvot ovat, vain niiden tyypillä.

Lopulta kysytty läpimenotodennäköisyys saadaan laskettua osina (ehdollistetaan sillä millainen tyyppi tulee, summa on kaikkien tyyppien $t$ yli):

$\begin{aligned} \mathbb{P}(\text{voitto}) = \sum_{t} A(t)B(t) \end{aligned}$

Koottu, lopullisen vastauksen antava koodi:

from fractions import Fraction
from collections import Counter

R.<z> = QQ['z']
def waysToSelect(multiplicities, mustHave=[]):
    g = R([1])
    for j,nj in enumerate(multiplicities):
        coeffs = [binomial(nj, k) for k in range(nj+1)]
        if j in mustHave: coeffs[0] = 0 #subtracts the 1
        g *= R(coeffs)
    return [g[k] for k in range(g.degree()+1)]

#list of partitions of n to at most maxPartsN parts that are at most maxSize
def getPartitions(n, maxPartsN, maxSize):
    if n==0: return [[]]
    canMakeMax = maxPartsN*maxSize
    if n>canMakeMax: return []
    if n==canMakeMax: return [[maxSize]*maxPartsN]
    ret = []
    for x in xrange(min(n, maxSize), -1, -1):
        ret += [[x] + sP for sP in getPartitions(n-x, maxPartsN-1, x)]
    return ret

def getAllPossibleTypes(boardN, suits, vals):
    return getPartitions(boardN, min(boardN, vals), suits)


def aProb(t, suits, vals):
    num = 1
    for i in xrange(len(t)):
        num *= (vals-i)*binomial(suits, t[i])
    num = num/reduce(lambda x,y: x*factorial(y), Counter(t).values(), 1)
    denom = binomial(suits*vals, sum(t))
    return Fraction(int(num), int(denom))

def bProb(t, suits, vals, pickN):
    n = suits*vals
    if len(t)==0 and pickN<=n: return Fraction(int(1), int(1))
    if (len(t) and pickN==0) or pickN>n-sum(t): return Fraction(int(0), int(1))
    deckMults = [suits-(t[i] if i<len(t) else 0) for i in range(vals)]
    w = waysToSelect(deckMults, range(len(t)))
    num = w[pickN] if len(w)>pickN else 0
    denom = binomial(suits*vals-sum(t), pickN)
    return Fraction(int(num), int(denom))

def winProb(suits, vals, boardN, pickN):
    ret = Fraction(int(0), int(1))
    n = suits*vals
    if boardN+pickN > n: return ret
    
    #for a type t, check if for some value all suits of that
    #value have come and, as a result, we can't win
    def someValueMaxed(t):
        for x in t:
            if x>=suits: return True
        return False
    
    ts = getAllPossibleTypes(boardN, suits, vals)
    for t in ts:
        if someValueMaxed(t): continue
        ret += aProb(t, suits, vals)*bProb(t, suits, vals, pickN)
    return ret

#the usual 4x13 deck, pick 13 in each phase of the game
p = winProb(4, 13, 13, 13)
print "{} \n= {}".format(str(p), float(p))

Tämä tulostaa läpimenotodennäköisyyden

$\begin{aligned} \frac{964444044208}{262190765217675} = 0.00367840584853 \end{aligned}$

Multijoukon valintojen todennäköisyyksiä

Olkoon $X = \{n_j \times x_j \space | \space j=1,2,\dots, m\}$ multijoukko. (Alkio $x_j$ esiintyy $n_j$ kertaa.) Kun siitä valitaan $k$ alkiota, mikä on todennäköisyys, että joitain kiinnitettyjä alkoita, jotka numerointia muuttamalla voidaan olettaa olevan $A = \{x_j \space | \space j=1, 2, \dots, p \}$ , saadaan kutakin vähintään yksi kappale? Valinta tehdään siten, että $X$ ajatellaan korttipakaksi, jossa jokaista $x_j$ -arvoista korttia on $n_j$ kappaletta (eli yhteensä $n = \sum_{j=1}^m n_j$ korttia) ja pakka sekoitetaan (permutaatioiden tasajakauma), jonka jälkeen valitaan $k$ päällimmäistä korttia.

Pohditaanpa seuraavaa polynomien tuloa:

$\begin{aligned} \prod_{j=1}^m (1+z+z^2+z^3+\dots +z^{n_j}) \end{aligned}$

Kun tämän kertoo auki, tulee termin $z^k$ kertoimeksi lukumäärä kuinka monella tavalla voidaan valita $k$ alkiota kaikista $n$ :stä. Eli valitaan vain lukumäärät jokaiselle $x_j$ :lle. Tämä ei kuitenkaan ota huomioon sitä mitkä alkiot $n_j$ :stä $x_j$ :stä valitaan. Tätä varten kuhunkin tulon polynomiin on laitettava kertoimiksi binomikertoimet: $1+{{n_1}\choose{1}}z + {{n_2}\choose{2}}z^2 + \dots + {{n_j}\choose{n_j}}z^{n_j}$ . Mutta tämähän on binomikaavan mukaan $(1+z)^{n_j}$ , joka myös nähdäänkin nyt uudessa valossa: jokaisen yksittäisen $x_j$ :n kohdalla tehdään valinta otetaanko se mukaan ( $z^1$ ) vai ei ( $1$ ). Koko tulo voidaan nyt kirjoittaa:

$\begin{aligned} G_0(z) := \prod_{j=1}^m (1+z)^{n_j} = (1+z)^{n_1+n_2+\dots n_m} = (1+z)^n \end{aligned}$

Eli tämä redusoituu edelleen tavallisiin binomikertoimiin ${{n}\choose{k}}$ , mikä onkin loogista sillä jokaista $x_j$ :tä pidetään erillisenä alkiona. Mutta tästä tullaankin siihen kysymykseen kuinka otamme huomioon, että jokaista $x_j$ , $j=1, 2, \dots, p$ pitää olla valinnassa ainakin yksi. Se tapahtuu jättämällä tulossa $G_0$ vastaavista polynomeista vakiotermit pois (eli ei sallita alkiota valittavan nolla kertaa). Näin saadaan generoivaksi funktioksi

$\begin{aligned} G_1(z) := \left ( \prod_{j=1}^p ((1+z)^{n_j}-1) \right) \left ( \prod_{j=p+1}^m (1+z)^{n_j} \right) = \left ( \prod_{j=1}^p ((1+z)^{n_j}-1) \right) (1+z)^{n- \sum_{j=1}^p n_j} \end{aligned}$

Seuraava Sage-koodi laskee $G_1$ :n kertoimet (multiplicities = luvut $n_j$ , mustHave = lista indekseistä, joita yllä merkittiin $j=1, 2, \dots, p$ ):

R.<z> = QQ['z']

def waysToSelect(multiplicities, mustHave=[]):
    g = R([1])
    #n = sum(multiplicities)
    for j,nj in enumerate(multiplicities):
        coeffs = [binomial(nj, k) for k in range(nj+1)]
        if j in mustHave: coeffs[0] = 0 #subtracts the 1
        g *= R(coeffs)
    return [g[k] for k in range(g.degree()+1)]

#an example: how to select from {0,0,1,1,1,2,3,3}
#so that have at least one 0 and at least one 1
a = waysToSelect([2,3,1,2], range(2))
print a

Tällä, kombinaatiot läpikäyvällä koodilla voi vielä tarkastaa tuloksen:

import itertools

def check(multiplicities, k, mustHave):
    deck = []
    for j, nj in enumerate(multiplicities): deck += [j]*nj
    
    def isGood(a):
        for x in mustHave:
            if x not in a: return False
        return True
    goodCount = 0
    totCount = 0
    for hand in itertools.combinations(deck, k):
        if isGood(hand): goodCount += 1
        totCount += 1
    return goodCount
    
mults = [2,3,1,2]
musts = range(2)
print [check(mults, k, musts) for k in range(sum(mults)+1)]

Koodeja voi suorittaa esim. osoitteessa https://sagecell.sagemath.org/

Todennäköisyydet saa, kun jakaa lukumäärät luvulla ${{n}\choose{k}}$ .
Menetelmää voi myös yleistää valitsemalla kuhunkin tulon polynomiin vain ne termit, joita vastaavat lukumäärät kyseiselle alkiolle sallitaan valinnassa esiintyvän.

Binäärijonot, joissa ei ole parillisen pituisia 0-putkia.

Tutkitaan binäärijonoja $\{0, 1\}^*$ ja näissä olevia yhtenäisiä $0$ -putkia. Esimerkiksi jonossa $011100011$ on kaksi $0$ -putkea: yhden ja kolmen pituiset. Näistä kumpikaan ei ole parillisen pituinen, joten jono on halutun kaltainen. Merkitään tällaisten $n$ :n pituisten jonojen lukumäärää $a_n$ :llä.

Jono $a_n$ alkaa $1, 2, 3, 6, 10, 19, 33, 61, 108, 197, 352, \dots$ ja tämän voi tarkastaa esim. seuraavalla Python-koodilla:

def hasEven0Run(s):
    run = 0
    isRun = False
    isOne = False
    for c in s:
        isOne = c=='1'
        if isRun:
            if isOne and run%2==0: return True
        isRun = not isOne
        if isRun: run += 1
        else: run=0
    return isRun and run>0 and run%2==0

def countBitStrings(n, tester):
    count = 0
    pad = n*'0'
    for i in xrange(2**n):
        s = (pad+(bin(i)[2:]))[-n:]
        if tester(s): count += 1
    return count

def getBitStrings(n, tester):
    ret = []
    pad = n*'0'
    for i in xrange(2**n):
        s = (pad+(bin(i)[2:]))[-n:]
        if tester(s): ret.append(s)
    return ret

runTester = lambda s: not hasEven0Run(s)
print [countBitStrings(n, runTester) for n in range(11)]

Kuinka sitten tämä lukumäärä $a_n$ voitaisiin laskea yleiselle $n$ ? Merkitään kaikkien halutunlaisten (kaikkien äärellisen pituisten) binäärijonojen joukkoa $A$ ja lisäksi $Z_{odd} = \{\cdot, \cdot\cdot\cdot, \}$ sekä triviaalin symbolin yksiötä $E = \{\epsilon\}$ . Tällöin saadaan näiden joukkojen välille (rekursio-)yhtälö

$\begin{aligned} A = (E+Z_{odd})\times (E+zA) \end{aligned}$

joka taas johtaa suoraan $A$ :n generoivan funktion yhtälöön (huom. $Z_{odd}$ :n gen.funktio on $z+z^3+z^5+\dots = \frac{z}{1-z^2}$ )

$\begin{aligned} A(z) = (1+ \frac{z}{1-z^2} )(1+zA(z)) \end{aligned}.$

Tästä voidaan ratkaista

$\begin{aligned} A(z) = \frac{1+z-z^2}{1-z-2z^2+z^3} \end{aligned} .$

Tämä rationaalifunktio voitaisiin jakaa osamurtohajotelmaan ja saada tarkka kaava (käyttämällä geometrisen summan kaavaa $\frac{1}{1-cz} = \sum_{j=0}^\infty (cz)^j$ ) mutta asymptoottista kaavaa varten riittää tutkia $A(z)$ :n lähimpänä origoa olevaa positiivista napaa. Tämä saadaan nimittäjän nollakohtana

$\begin{aligned} \alpha \approx 0.554958132 \end{aligned}$

Koska $\alpha$ on nimittäjän yksinkertainen nollakohta eikä se ole osoittajan nollakohta, niin jonolle $a_n$ saadaan asymptoottinen kaava

$\begin{aligned}a_n &\sim -\frac{1}{\alpha} \frac{1+\alpha-\alpha^2}{-1-4\alpha + 3\alpha^2} \left( \frac{1}{\alpha} \right)^n \end{aligned} \approx a \beta^n,$

missä $a\approx 0.97869358$ ja $\beta = \frac{1}{\alpha} \approx 1.8019377$ .