Kahden vääntyneen kolikon peli

Pelissä on kaksi kolikkoa, kutsutaan niitä ensimmäinen ja toinen kolikko. Todennäköisyydet että näillä tulee kruuna (H) ovat vastaavasti x ja y.

Peli etenee seuraavasti: Pelaaja heittää kolikkoja toistuvasti ja pistää merkille minkä tuloksen hän sai (HT, TH, HH tai TT). Jos hän jossain vaiheessa saa kaksi kertaa putkeen tuloksen HH tai kaksi kertaa putkeen tuloksen TT, niin peli päättyy. Jos taas hänen viime tuloksensa on HH ja hän saa heti tämän jälkeen tuloksen TT, niin hän pääsee bonuspeliin. Bonuspelissä hän heittää molempia kolikkoja (erikseen) kaksi kertaa. Jos ensimmäisellä kolikolla tulee ensimmäisellä heitolla H, pelaaja voittaa 5 pistettä, mutta jos tämän jälkeen toisella heitolla tulee H, pelaaja menettää 3 pistettä (eli voittaa vain 2 pistettä yhteensä). Mikäli ensimmäisellä heitolla tulee T, pelaaja ei voita eikä häviä mitään. Vastaavat luvut toiselle kolikolle ovat 3 ja -4 eli HT antaa 3 pistettä ja HH 3-4 = -1 pistettä. Bonuspelin jälkeen peli jatkuu normaalisti (bonuspelin heittoja ei rekisteröidä, eli pelissä viimeksi on tullut TT).

Miten vääntyneet kolikot peliin kannattaa ottaa (eli kuinka muuttujat x ja y väliltä [0, 1] valita), kun halutaan maksimoida pelaajan voittojen odotusarvo?

Ratkaisu

Lasketaan todennäköisyydet

  • a = \mathbb{P}(HT \wedge TH) = x+y-2xy
  • b = \mathbb{P}(HH) = xy
  • c = \mathbb{P}(TT) = 1-x-y+xy

Muodostetaan Markovin ketju tiloina ”HT tai TH”, ”HH”, ”TT” ja ”End”.

Pelin Markovin ketju, sen siirtymätodennäköisyysmatriisi ja fundamentaalimatriisi. Olemme kiinnostuneet fundamentaalimatriisin ensimmäisen rivin toisesta alkiosta, joka kertoo odotusarvon kuinka monta kertaa ketju käy tilassa HH ennen absorboitumista, kun se alkaa tilasta ”HT tai TH”. (Aloitetaan tästä tilasta, koska silloin ei ole väliä mitä on tullut viimeksi.)

Siirtymätodennäköisyydet ovat muuten tilasta riippumattomat, paitsi tiloista ”HH” ja ”TT” ei mennäkkään vastaavilla heitoilla niihin itseihinsä vaan tilaan ”End”. Tila ”End” on absorboiva.

Pelin voiton odotusarvo saadaan (iteroidun odotusarvon kaavan nojalla) kertomalla yhden bonuspelin voiton odotusarvo (5x -3x^2 + 3y - 4y^2) sillä kuinka monta kertaa siihen odotusarvoisesti päästään. Tämä odotusarvo taas saadaan kertomalla tilassa ”HH” vierailujen odotusarvo siirtymätodennäköisyydellä ”HH” -> ”TT”, joka on c. Näin ollen pelin voiton odotusarvo on

\begin{aligned} \frac{( 5x -3x^2 + 3y - 4y^2 )cb(c+1)}{1-a(b+1)(c+1)-bc} = \frac{ (5x -3x^2 + 3y - 4y^2)xy(1-x-y+xy)(2-x-y+xy) }{1-(x+y-2xy)(1+xy)(2-x-y+xy)-xy(1-x-y+xy)} \end{aligned}

Kuvasta nähdään, että tämän maksimikohta joukossa [0,1]^2 on jossain (0.7, 0.3):n kohdilla:

Odotusarvo riippuu muuttujista x ja y. Kuvaaja ja tasa-arvo käyrä-kuvaaja.

Tässä javascript-funktiot odotusarvon ja sen gradientin laskemiseksi:

/*
(x(5x - 3x^2 + 3y-4y^2)y(1 - x - y + xy)(2 - x - y + xy))
/
(1 - xy(1 - x - y + xy)
- (x + y - 2xy) (1 + xy) (2 - x - y + xy))
, x\in[0,1], y\in [0,1]
*/

const f = (x, y) => {
    let g = 5*x - 3*x**2 + 3*y-4*y**2;
    let num = g*x*y*(1 - x - y + x*y)*(2 - x - y + x*y);
    let denom = 1 - x*y*(1 - x - y + x*y)
        - (x + y - 2*x*y)*(1 + x*y)*(2 - x - y + x*y);
    
    return num/denom;
};

const gradF = (x, y) => {
    let p1 = 5*x - 3*x**2 + 3*y-4*y**2;
    let p2 = x*y*(1-x-y+x*y); //= xy - x^2y - xy^2 + x^2y^2
    let p3 = 2 - x - y + x*y;
    let num = p1*p2*p3;
    let denom = 1 - x*y*(1 - x - y + x*y)
            - (x + y - 2*x*y)*(1 + x*y)*(2 - x - y + x*y);
    
    let dxp1 = 5 - 6*x;
    let dxp2 = y - 2*x*y - y**2 + 2*x*y**2;
    let dxp3 = -1+y;
    let dxnum = dxp1*p2*p3 + p1*dxp2*p3 + p1*p2*dxp3;
    let dxdenom = (y-2)*(y-1)**2 + 3*x**2*y*(1-3*y+2*y**2)
                + x*(2-8*y+14*y**2-6*y**3);
    
    let dyp1 = 3-8*y;
    let dyp2 = x - x**2 - 2*x*y + 2*x**2*y;
    let dyp3 = -1 + x;
    let dynum = dyp1*p2*p3 + p1*dyp2*p3 + p1*p2*dyp3;
    let dydenom = 2*(y-1) + x**2*(-4+14*y-9*y**2)
                + x*(5-8*y+3*y**2) + x**3*(1-6*y+6*y**2)
    
    return [
        (dxnum*denom-num*dxdenom)/denom**2,
        (dynum*denom-num*dydenom)/denom**2
    ];
    
};

Jyrkimmän nousun menetelmällä (gradient descent) tai muuten löydetään maksimi

(x,y) =  (0.71046, 0.28314)

ja odotusarvo

1.2863.

Värisuorien todennäköisyyksistä

Jos tavallisesta 52 kortin korttipakasta nostetaan satunnainen 7 kortin käsi, mikä on todennäköisyys, että käsi sisältää (ainakin) r peräkkäistä saman maan korttia?

Merkitään A_{i, j} = ”käsi sisältää kortit (i, j), (i+1, j), \dots, (i+r-1, j)”, arvoille i=1,2,\dots, 13-r+1 ja j= ♠, ♥, ♣, ♦. Merkitään näiden kaikkien tapahtumien joukkoa \cal{A}:lla. Kysytty todennäköisyys on tällöin

\begin{aligned} \mathbb{P} \left( \bigcup_{A\in\cal{A}} A \right) \end{aligned} .

Inkluusio-eksluusion periaatteen nojalla tämä voidaan laskea

\begin{aligned}   \mathbb{P} \left( \bigcup_{A\in\cal{A}} A \right)    =  \sum_{j=1}^{|\cal{A}|} (-1)^{j+1} \sum_{ \substack{ I\subset \cal{A} \\ |I| = j}} \mathbb{P} \left( \bigcap_{A\in I} A \right)  \end{aligned}

Merkitään tapahtumalle A \in \cal{A} symbolilla \bar {A} niiden korttien joukkoa, jotka siinä vaaditaan kädessä esiintymään. Huomataan, että tämä relaatio on siinä mielessä käänteinen, että kun tapahtumia A leikataan, vastaa tämä joukkojen \bar{A} yhdistettä. (Leikkautapahtuma tapahtuu jos ja vain jos käsi sisältää kaikki vaaditut kortit.)

Todennäköisyys \begin{aligned} \mathbb{P} \left( \bigcap_{A\in I} A \right)  \end{aligned} riippuu vain vaadittujen korttien unionin koosta t = \left(   \bigcup_{A\in I} \bar{A}  \right) sillä siinä vaaditaan leikkauksen korttien esiintyvän kädessä eikä sillä ole väliä mitkä kortit ne ovat. Hypergeometrisen jakauman mukaan tämä todennäköisyys on

p_t := \begin{aligned} \frac{ {{t} \choose {t}}  {{52-t} \choose {7-t}} }{ {{52} \choose {7}} } =  \frac{  {{52-t} \choose {7-t}} }{ {{52} \choose {7}} }  \end{aligned}

Täytyy siis selvittää kuinka monta t:n kokoista unionia \bigcup_{I\subset \cal{A}} \bar{A} tasolla |I|=j on. Merkitään tätä c_{j, t}.

Koska p_t = 0 , kun t>7, niin tarvitaan vain arvot, joissa t\leq 7. Lisäksi, koska tason j unioni koostuu j joukosta ja ensimmäinen tuo r alkiota ja seuraavat jokainen ainakin yhden lisää (joukot \bar{A} ovat pötköjä, jotka menevät maksimissaan r-1:n verran päällekkäin), niin saadaan myös yläraja j\leq 7-r+1.

Lukujen c_{j, t} laskeminen voidaan tehdä brute-forcena. Tapaus r=2 on melko raskahko, sillä siinä täytyy laskea arvoon j=6 asti ja {{4\times 12}\choose{6}} = 12 271 512. Tapaus r=2 on kuitenkin toisella tapaa yksinkertaisempi ja sen voi ratkaista tutkimalla kuinka monella tapaa neljän 13-ketjun graafista voidaan valita aligraafi, joka on kokoelma ketjuja (nämä aliketjujen pituudet vastaavat sitä tyyppiä millainen yhdiste A’ista muodostuu). Tässä joukot A vastaavat graafin sivuja ja se että joukko on yhdisteessä mukana vastaa sitä, että sivu on aligraafissa mukana.

Saadaan seuraavat tulokset:

r P(väh. r pituinen värisuora) approks
2 686898/1194505 0.575
3 1074209/16723070 0.0642
4 151/30940 0.00488
5 407/1454180 0.000280
6 361/33446140 0.0000108
7 1/4778020 2.093e-7

Javascript-koodi:

function binCoeff(n, k) {
    if (k < 0 || k > n) return 0;
    if (k === 0 || k === n) return 1;
    let ret = 1;
    for (let i=0; i<Math.min(k, n-k); i++){
        ret = ret*(n-i)/(i+1);
    }
    return ret;
}

function* generCombinations(arr, r) {
    const result = [];
    result.length = r;
    function* makeThem(len, start) {
        if(len===0) {
            yield result.slice();
        } else {
            for (var i=start; i<=arr.length-len; i++) {
                result[result.length-len] = arr[i];
                yield* makeThem(len-1, i+1);
            }
        }
    }
    yield* makeThem(r, 0);
}

function linkBruteGeneral(j, k, m, linkLen) {
    let ret = {};
    let runStarts = [];
    for (let i=0; i<k; i++) {
        runStarts = [...runStarts,
...new Array(m-(linkLen-1)).fill(null).map((_,z)=>[z, i])];
    }
    for (let comb of generCombinations(runStarts, j)) {
        let cardSet = new Set();
        for (let card of comb) {
            for (let follow=0; follow<linkLen; follow++) {
                cardSet.add((card[0]+follow)+","+card[1]);
            }
        }
        let key = cardSet.size;
        if (!(key in ret)) ret[key] = 0;
        ret[key] += 1;
    }
    return ret;
}

function calcProb(ds) {
    const calcPart = (j) => {
        let s = 0;
        let d = ds[j];
        for (let k in d) {
            if (k>7) continue;
            s += d[k]*binCoeff(52-k, 7-k);
        }
        return s;
    };
    
    let ret = 0;
    for (let i=1; i<=ds.length; i++) {
        ret += (i%2?1:-1)*calcPart(i);
    }
    return ret;
}

function makeTheFormula(ds) {
    const makePart = (j) => {
        let f = "(";
        let d = ds[j];
        let isFirst = true;
        for (let k in d) {
            if (k>7) continue;
            if (isFirst) {
                isFirst = false;
            } else {
                f += "+";
            }
            f += d[k]+`*binomial(${52-k}, ${7-k})`;
        }
        return f+")";
    };
    
    let ret = "";
    for (let i=1; i<=ds.length; i++) {
        ret += (i===1?'':(i%2?'+':'-')) + makePart(i);
    }
    return ret;
}

let testK = 4;
let testM = 13;
let testR = 3;

let maxNeedJ = 7-testR+1;

let ds = [];
for (let j=0; j<=maxNeedJ; j++) {
    ds.push(linkBruteGeneral(j, testK, testM, testR));
}

console.log(calcProb(ds));
console.log(makeTheFormula(ds));

Koodi toiselle tavalle tapauksessa r=2:

/** list of partitions of n to at most maxPartsN parts that are at most maxSize */
function getPartitions(n, maxPartsN, maxSize) {
    if (n===0) return [[]];
    let canMakeMax = maxPartsN*maxSize;
    if (n>canMakeMax) return [];
    if (n===canMakeMax) return [new Array(maxPartsN).fill(maxSize)];
    let ret = [];
    for (let x=Math.min(n, maxSize); x>=0; x--) {
        let toConcat = getPartitions(n-x, maxPartsN-1, x).map(p=>[x].concat(p));
        ret = ret.concat(toConcat);
    }
    return ret
}

function factorial(n) {
    let ret = 1;
    for (let i=2; i<=n; i++) ret *= i;
    return ret;
}

function binCoeff(n, k) {
    if (k < 0 || k > n) return 0;
    if (k === 0 || k === n) return 1;
    let ret = 1;
    for (let i=0; i<Math.min(k, n-k); i++){
        ret = ret*(n-i)/(i+1);
    }
    return ret;
}

/** In how many ways  can subchains be fitted into chains.
* Chain length is the number of edges in it.
*/
const linkageMemo = {};
function countLinkages(chains, subchains) {
    //debugger;
    let memoKey = chains.join(",")+";"+subchains.join(",");
    if (memoKey in linkageMemo) return linkageMemo[memoKey];
    if (subchains.length===0) return 1;
    let s = 0;
    let currSub = subchains[0];
    let restSubs = subchains.slice(1);
    for (let cI=0; cI<chains.length; cI++) {
        let cLen = chains[cI];
        let maxInd = cLen-currSub;
        for (let ind=0; ind<=maxInd; ind++) {
            let restChains = chains.slice(0, cI).concat(chains.slice(cI+1));
            if (ind>1) restChains.push(ind-1);
            if (ind<maxInd-1) restChains.push(maxInd-ind-1);
            restChains.sort((a,b)=>b-a);
            s += countLinkages(restChains, restSubs);
        }
    }
    linkageMemo[memoKey] = s;
    return s;
}


function permuFactor(arr) {
    let counts = {};
    for (let x of arr) {
        if (!(x in counts)) counts[x] = 0;
        counts[x] += 1;
    }
    let ret = 1;
    for (let val of Object.values(counts)) ret *= factorial(val);
    return ret;
}

function countAllLinkages(j, k, m) {
    let chains = new Array(k).fill(m-1);
    let linksN = k*(m-1);
    let partis = getPartitions(j, linksN, linksN);
    let ret = {};
    for (let parti of partis) {
        let key = parti.reduce((a,b)=>a+b+1,0);
        if (!(key in ret)) ret[key] = 0;
        ret[key] += countLinkages(chains, parti)/permuFactor(parti);
    }
    return ret;
}


/* Given the amounts of linkages calculate the numerator of probability of having a 2-run
* ds = [
    {
      "0": 1
     },
     {
      "2": 48
     },
     {
      "3": 44,
      "4": 1084
     },
     {
      "4": 40,
      "5": 1944,
      "6": 15312
     }, ...
]
*/
function calcProb(ds) {
    const calcPart = (j) => {
        let s = 0;
        let d = ds[j];
        for (let k in d) {
            if (k>7) continue;
            s += d[k]*binCoeff(52-k, 7-k);
        }
        return s;
    };
    
    let ret = 0;
    for (let i=1; i<=6; i++) {
        ret += (i%2?1:-1)*calcPart(i);
    }
    return ret;
}

function makeTheFormula(ds) {
    const makePart = (j) => {
        let f = "(";
        let d = ds[j];
        let isFirst = true;
        for (let k in d) {
            if (k>7) continue;
            if (isFirst) {
                isFirst = false;
            } else {
                f += "+";
            }
            f += d[k]+`*binomial(${52-k}, ${7-k})`;
        }
        return f+")";
    };
    
    let ret = "";
    for (let i=1; i<=6; i++) {
        ret += (i===1?'':(i%2?'+':'-')) + makePart(i);
    }
    return ret;
}


let testK = 4;
let testM = 13;
let ds = [];
for (let j=0; j<=6; j++) {
    ds.push((countAllLinkages(j, testK, testM)));
}

console.log(makeTheFormula(ds));
console.log(calcProb(ds));

\begin{aligned} \end{aligned}

Ympyrät tasossa

Kysymys

Olkoon tasossa n toisiaan leikkaamatonta 1-säteistä ympyrää. (Ympyrällä tarkoitetaan sen kaarta.) Mikä on todennäköisyys, että kun valitaan näiltä piste satunnaisesti, tästä pisteestä ei näy muita ympyröitä?

Ympyrän 2 pisteestä A näkyy kaikki muut ympyrät. Pisteestä B taas ei näy mikään toinen ympyrä, sillä ne ovat kaikki B:hen piirretyn tangenttisuoran ”väärällä” puolella.

Kokeile:

https://Invisible-Circles–minkkilaukku2.repl.co

Ratkaisu

Vastaus on \begin{aligned} \frac{1}{n} \end{aligned}.

Ja tämä siis aivan riippumatta ympyröiden sijainneista (paitsi eivät saa leikata toisiaan, kuten sovittiin).

Todistuksen runko

Ympyröiden konveksiverho C. Kaariosuudet keltaisella ja janaosuudet ruskealla. Vihreät pisteet ympyröillä ovat niiden välisten tangenttijanojen päätepisteitä, ”kaikkein uloimmaiset” tangenttijanat ovat mukana C:ssä.
  1. Olkoon C_0 ympyröiden yhdisteen konveksiverho.
  2. Joukon C_0 reuna C on janoista ja ympyränkaarista koostuva derivoituva Jordanin käyrä.
  3. Ympyrän piste P kuuluu C:hen täsmälleen silloin, kun se on sellainen piste, josta ei näy mikään muu ympyrä. (Pois lukien äärellisen monta pistettä, jotka ovat C:n janojen päätepisteitä: janat ovat ympyröiden välisiä tangentteja, joten niiden päätepisteet leikkaavat kahta ympyrää, joten niistä on näköyhteys. Näiden joukko on kuitenkin äärellisenä nollamittainen eikä vaikuta todennäköisyyteen.)
  4. Käyrän C kokonaiskaarevuus on \begin{aligned} \int_{a}^b \kappa (s) ds  = 2\pi \end{aligned}, sillä kierrosluku on 1.
  5. Toisaalta \begin{aligned}   \int_{a}^b \kappa (s) ds  \end{aligned} on käyrällä C olevien ympyränkaarien yhteenlaskettu mitta, sillä janaosuudella kaarevuus \kappa = 0 ja kaariosuudella \kappa = \frac{1}{1} = 1 , sillä jokaisen ympyrän säde on 1 (ks. kaarevuudesta Wikipediassa).
  6. Siis suotuisan osan mitta on 2\pi. Kaikkien ympyröiden mitta on 2\pi n, joten väite seuraa.

Huomioita

  • Yleistyy kolmiulotteiseksi: pallot 3D-avaruudessa. Gauss-Bonnet:n lause.
  • Jos ympyröiden säteet ovat eri suuruisia, ympyröiden paikat vaikuttavat tulokseen. Tämä johtuu siitä, yo. integraalissa kaarevuus ei ole joka ympyrällä vakio vaan tulos riippuu siitä millaisia ympyröitä ulkokaaressa on mukana minkäkin verran. Ks. esimerkkikuva alla.
Erisäteisille ympyröille voi käydä näin. Selvästi jälkimmäisessä tapauksessa suotuisaa mittaa on vähemmän kuin ensimmäisessä.

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}