Search code examples
algorithmgraphprobabilityfractionsfractals

total path probability of weighted circular graph


Let's say there is a game where for each move there are probable paths, depending on the throw of fancy dice. Depending on the results there might be transitions forward, backward or staying on one place. Eventually (even after infinite amount of throws) the graph leads to final states. Each edge is weighted with probability .

For the case where there are no loops I can simply sum+multiply and re-normalize the probabilities for each outcome if I start at the same vertex(cell).

However, if I have loops it starts getting confusing. For example, let's say we have same probability on each edge:

  start0
   /\   ^
  /  \  |
end1  tr2
      /
     end2

the graph starts at start0 and has 50% chance of terminating at end1 or going to transition tr2. From tr2 there is again 50% chance terminating at end2 or going back to start0.

How could I calculate the total probability for reaching each stop end1 and end2. If I try using a converging series like this:

pEnd1=1/2 + 1/2*1/2+1/8+.. ->lim->1. which makes no sense since end2 is getting no probability. Obviously I have a mistake there.

So my question is how can I calculate the probabilities of reaching the final nodes if I have the probabilities of each edge but I could have loops.

example 1) simple fork with a loop all edges are 50% probable

start0-> p=50% ->end1
start0-> p=50% ->tr1
tr2-> p=50% ->start0
tr2-> p=50% ->end2

example 2) more loops

start0-> p=1/3 ->e1
start0-> p=1/3 ->tr1
start0-> p=1/3 ->start0
tr1-> p=1/3 ->tr2
tr1-> p=2/3 ->start0
tr2-> p=7/9 ->start0
tr2-> p=2/9 ->end2

example 3) - degenerate test case - since all paths end at e1 - it should end up having 100% probability

start0-> p=1/3 ->e1
start0-> p=2/3 ->tr1
tr1-> p=3/4 ->start0
tr2-> p=1/4 ->e1

Solution

  • This is not really a programming problem, although you could write a simulation and perform it 100000 times to see what the distribution would be.

    You wrote:

    pEnd1=1/2 + 1/2*1/2+1/8+.. ->lim->1. which makes no sense since end2 is getting no probability. Obviously I have a mistake there.

    Indeed, there is a mistake. You did not take into account the probability to go from tr2 to start0 (50%). The probability that the path will cycle once to start0 and then end up in end1 is 1/2 (going to tr2) * 1/2 (going back to start0) * 1/2 (going to end1). The number of decisions (of 50%) is always odd when ending up in end1. And it is even when ending up in end2. So the formula would be:

    pEnd1 = 2-1 + 2-3 + 2-5 + ... -> lim = 2/3

    pEnd2 = 2-2 + 2-4 + 2-6 + ... -> lim = 1/3

    Simulation

    To make this a programming question, here is a simulation in JavaScript:

    function run(transitions, state) {
        while (transitions[state][state] != 1) { // not in sink
            let probs = transitions[state];
            let rnd = Math.random(); // in range [0, 1)
            for (let i = 0; i < probs.length; i++) {
                rnd -= probs[i];
                if (rnd < 0) {
                    state = i; // transition
                    break;
                }
            }
        }
        return state;
    }
    
    // Define graph
    let names = ["start0", "end1", "tr2", "end2"]
    let transitions = [
        [0.0, 0.5, 0.5, 0.0],
        [0.0, 1.0, 0.0, 0.0], // sink
        [0.5, 0.0, 0.0, 0.5],
        [0.0, 0.0, 0.0, 1.0]  // sink
    ];
    
    // Start sampling
    let numResults = [0, 0, 0, 0];
    let numSamples = 0;
    setInterval(function () {
        let endstate = run(transitions, 0);
        numSamples++;
        numResults[endstate]++;
        document.querySelector("#" + names[endstate]).textContent = (100 * numResults[endstate] / numSamples).toFixed(4) + "%";
    }, 1);
    <div>Arriving in end1: <span id="end1"></span></div>
    <div>Arriving in end2: <span id="end2"></span></div>

    You may also like to read about Absorbing Markov chains. From that we learn that the "absorbing probabilities" matrix B can be calculated with the formula:

    B = NR

    Where:

    • N is the "fundamental matrix" (I - Q)⁻¹     
    • I is the identity matrix of the same shape as Q
    • Q is the probability matrix for transitions between non-final states
    • R is the probability matrix for transitions to final states

    Here is a script (including the relevant matrix functions) to calculate B for the example problem you depicted:

    // Probabilities to go from one non-final state to another
    let Q = [
        [0.0, 0.5], // from start0 to [start0, tr2]
        [0.5, 0.0]  // from tr2    to [tr2, start0]
    ];
    // Probabilities to go from one non-final state to a final one
    let R = [
        [0.5, 0.0], // from start0 to [end1, end2]
        [0.0, 0.5]  // from tr2    to [end1, end2]
    ];
    // See https://en.wikipedia.org/wiki/Absorbing_Markov_chain#Absorbing_probabilities
    let N = inversed(sum(identity(Q.length), scalarProduct(Q, -1)));
    let B = product(N, R);
    console.log("B = (I-Q)⁻¹R:\n" + str(B));
    
    // Generic matrix utility functions:
    // cofactor is copy of given matrix without given column and given row 
    function cofactor(a, y, x) {
        return a.slice(0, y).concat(a.slice(y+1)).map(row => row.slice(0, x).concat(row.slice(x+1)));
    } 
    
    function determinant(a) {
        return a.length == 1 ? a[0][0] : a.reduceRight((sum, row, y) =>
            a[y][0] * determinant(cofactor(a, y, 0)) - sum
        , 0);
    } 
    
    function adjoint(a) {
        return a.length == 1 ? [[1]] : transposed(a.map((row, y) => 
            row.map((_, x) => ((x + y) % 2 ? -1 : 1) * determinant(cofactor(a, y, x)))
        ));
    } 
    
    function transposed(a) {
        return a[0].map((_, x) => a.map((_, y) => a[y][x]));
    }
    
    function scalarProduct(a, coeff) {
        return a.map((row, y) => row.map((val, x) => val * coeff));
    }
    
    function inversed(a) {
        return scalarProduct(adjoint(a), 1 / determinant(a));
    }
    
    function product(a, b) {
        return a.map((rowa) =>
            b[0].map((_, x) =>
                b.reduce((sum, rowb, z) =>
                    sum + rowa[z] * rowb[x]
                , 0)
            )
        );
    }
    
    function sum(a, b) {
        return a.map((row, y) => row.map((val, x) => val + b[y][x]));
    }
    
    function identity(length) {
        return Array.from({length}, (_, y) => 
            Array.from({length}, (_, x) => +(y == x))
        );
    }
    
    function str(a) {
        return a.map(row => JSON.stringify(row)).join("\n");
    }

    The output is:

    [
        [2/3, 1/3] // probabilities when starting in start0 and ending in [end1, end2]
        [1/3, 2/3] // probabilities when starting in tr2 and ending in [end1, end2]
    ]