Assigning probabilities to processes, when you know which probabilities exist but not which one belongs to which process

Question

I have multiple stochastic processes that can have 2 different realisations. I think of them as coins that can take the values heads or tails. For the example I assume that there are 3 coins. I know the 3 probabilities of showing heads that these coins have, but I don't know which belongs to which coin. For example I know that the probabilities are 60%, 50%, and 40% and each occurs exactly once i.e. if Coin A has a probability of 60% for heads, Coin B and C cannot have 60%.

I observe realisations for each coin over multiple rounds. With 3 coins there are 3! = 6 different ways the probabilities could be assigned to the coins. For each of these 6 possibilities I want to know probability that it is the true one given the data on the realisations of the 3 coins that I have. For example I want a result like with 20% probability Coin A has 60% probability of heads, Coin B has 50% and Coin C has 40% and with 15% Probability Coin A has 50% probability of heads, Coin B has 60% and Coin C has 40% and so on.

How can I calculate these probabilities?

Answer 1

The idea is to iterate over all possible permutations, compute the posterior probability (likelihood * prior), and then normalize.

Here's a Python program to compute the posterior probabilities:

from math import factorial
from itertools import permutations

# these are our fixed probabilities
P = [0.4, 0.5, 0.6]

# need a prior
PRIOR = 1.0/3.0

# here's the observed data: 5 "rounds", and the number of successes of each coin.
observed_data = (5, (3,2,1))


# some helper functions
def C(n, r):
    return factorial(n) / factorial(r) / factorial(n-r)


def likelihood(observed_data, permutation):
    """
    permutation is a permutation of (0,1,2)

    """

    N, successes = observed_data
    prob = 1.0
    for observed_i, proposed_i in enumerate(permutation):
        p = P[proposed_i]
        r = successes[observed_i]
        prob *= C(N, r) * p**r * (1 - p)**(N - r)
    return prob 


def unnormalized_posterior_probability(observed_data, permutation):
    """
    observed_data is a tuple, first element is the number of rounds
    the second element is a list representing the count of successes
    ex: (5, [2, 0, 3])

    """
    return likelihood(observed_data, permutation) * PRIOR 

# compute the posterior probabilities
running_sum = 0
posterior_probabilities = {}

for possible_order in permutations((0,1,2)):
    posterior_probabilities[possible_order] = unnormalized_posterior_probability(observed_data, possible_order)
    running_sum += posterior_probabilities[possible_order]

for possible_order in permutations((0,1,2)):
    posterior_probabilities[possible_order] = posterior_probabilities[possible_order] / running_sum

The output looks like:

{(0, 1, 2): 0.063,
 (0, 2, 1): 0.094,
 (1, 0, 2): 0.094,
 (1, 2, 0): 0.213,
 (2, 0, 1): 0.213,
 (2, 1, 0): 0.320}

Which suggests that (2, 1, 0) is the most probable order (which makes sense, given P and the observed data.