Imagine that you have a friend with 7 coins: some combination of pennies, nickels, dimes, or quarters. Your friend tells you that if you can guess the denomination of each of the coins, he will give them to you. He then proceeds to flip them onto the table, and to add the values of all the coins that land heads up, and tells you their sum. He agrees to repeat this process until you're confident in your answer.

#number of coins

physical_N = 7

#denominations of available coins, in cents

denominations = [1, 5, 10, 25]

#choose the value of the actual coins

physical_coins = np.random.choice(denominations, size = physical_N)

#probability that the coin will land heads (in this case, fair)

physical_pheads = np.ones(shape=physical_N)/2

#model of the physical system

Number_of_samples = 300

physical_measurement = np.array([np.sum(np.random.binomial(n=1,p=physical_pheads,

size=physical_N)*physical_coins)

for i in range(Number_of_samples)])

Now we can use our simulation to generate some test data:

samples: [12 41 11 45 22 27 36 17 41 21 ... ]

With seven coins, each able to take on four values, there are a total of 4^7 = 16384 different permutations. Luckily for us, symmetry dictates that not all of these solutions are unique: we don't care which coins are nickels and which are dimes, as long as we guess the correct number of each. In this case, we're interested in the combinations:

A large, but more manageable number. Each of these combinations represents a hypothesis (ie. two dimes, two nickels, two pennies, and a quarter) and together they form a 'hypothesis space'. For example, if we simplify the problem down to two coins, our probability space looks like a two dimensional matrix, with each coin represented along an axis, and the probability of the combination of those two coins making up the body of the table.

We translate the problem of guessing the coins to be a problem of filling in the question marks in this probability space with the probability that each particular combination represents the true value of the coins. As we have already seen, because we don't care about the order of the coins, our hypothesis space is redundant. In the picture above, the probability that the coins are worth 5 and 1 cent will be the same as the probability that the coins are worth 1 and 5 cents. To reduce the hypothesis space, we want to eliminate one of the symmetric hypotheses (5,1) and (1,5). One way we can do this is by saying that we will only look at coordinates that are in ascending order. This stipulation breaks the inherent symmetry of the problem and lets us look at the region of the hypothesis space with unique values. For the two dimensional case this makes our matrix:

We can generate this hypothesis space as a numpy array with seven dimensions (one for each coin) and each dimension of length four, one for each coin denomination.

potential_num_coins = physical_N

potential_coin_values = denominations

hypothesis_space_shape = [len(potential_coin_values)

for i in range(potential_num_coins)]

hypothesis_space = np.zeros(shape=hypothesis_space_shape)

it = np.nditer(hypothesis_space, flags=['multi_index'], op_flags=['readwrite'])

while not it.finished:

#are coordinates sorted?

if (sorted(list(it.multi_index)) == list(it.multi_index)):

it[0] = 1

it.iternext()

#normalize the state space

count = np.sum(hypothesis_space)

hypothesis_space /= count

We iterate over all elements of the array, setting all to zero except the ones whose coordinates come in ascending order. These we set such that their sum equals one, and they form a uniform prior hypothesis over the state space for us to apply Bayes rule and our data to.

In order to apply Bayes rule we need to define a likelihood function. The likelihood function gives us the probability of getting a certain value in each toss, given a specific hypothesis for the coins' individual values. To do this we can calculate all the different permutations of the coins falling heads up, looking for combinations that sum to the measured value. In python:

def likelihood(sum, coin_set_hypothesis):

p_heads = np.ones(len(coin_set_hypothesis))/2

p = 0

#the highest index we want is 2**len(bulbs)-1

for state in range(0,2**len(coin_set_hypothesis)):

#convert the state number to an array of individual coin states

st = np.array(map(int, bin(state)[2:].zfill(len(coin_set_hypothesis))))

value = np.sum(st*coin_set_hypothesis)

#accumulate the probabilities of permutations that result in the measured value

if value == sum:

p += np.product(st*p_heads + (1-st)*(1-p_heads))

return p

Now that we have data, a hypothesis space (with prior), and a likelihood function, we are ready to apply Bayes rule. For our problem, this takes the form:

So to apply the formula, when we get a new datapoint, we crawl over the entire hypothesis space, updating the probabilities of each hypothesis, and then we normalize the probabilities so that they sum to zero. In python, this becomes:

for datapoint in physical_measurement[:]:

it = np.nditer(hypothesis_space, flags=['multi_index'], op_flags=['readwrite'])

while not it.finished:

#only work on non-zero probabilities

if it[0]>0:

it[0] *= likelihood(datapoint, [potential_coin_values[i] for i in it.multi_index])

it.iternext()

#normalize

hypothesis_space /= np.sum(hypothesis_space)

# if we've reached a sufficiently high confidence, we can stop taking measurements

if np.max(hypothesis_space) > .999: break

measurement: 11 potential solutions: 72 prob of best guess: 0.0431965442765

measurement: 37 potential solutions: 30 prob of best guess: 0.0939334637965

`measurement: 15 potential solutions: 18 prob of best guess: 0.212076583211`

measurement: 31 potential solutions: 17 prob of best guess: 0.265437788018

measurement: 22 potential solutions: 12 prob of best guess: 0.246006006006

measurement: 26 potential solutions: 12 prob of best guess: 0.341038082771

measurement: 26 potential solutions: 12 prob of best guess: 0.400122469488

measurement: 36 potential solutions: 12 prob of best guess: 0.577365934659

measurement: 47 potential solutions: 9 prob of best guess: 0.64290385683

measurement: 40 potential solutions: 9 prob of best guess: 0.350104310773

measurement: 16 potential solutions: 9 prob of best guess: 0.379965178792

measurement: 21 potential solutions: 9 prob of best guess: 0.571889851959

measurement: 55 potential solutions: 4 prob of best guess: 0.518260091103

measurement: 35 potential solutions: 4 prob of best guess: 0.529247615951

measurement: 16 potential solutions: 4 prob of best guess: 0.521924746104

measurement: 21 potential solutions: 4 prob of best guess: 0.526892338442

measurement: 37 potential solutions: 4 prob of best guess: 0.539400203991

measurement: 30 potential solutions: 4 prob of best guess: 0.554541922563

measurement: 26 potential solutions: 4 prob of best guess: 0.510603887618

measurement: 55 potential solutions: 4 prob of best guess: 0.50891404123

measurement: 36 potential solutions: 4 prob of best guess: 0.578816353773

measurement: 36 potential solutions: 4 prob of best guess: 0.645292356054

measurement: 16 potential solutions: 4 prob of best guess: 0.580580572545

measurement: 47 potential solutions: 4 prob of best guess: 0.581762705335

measurement: 31 potential solutions: 4 prob of best guess: 0.513238821848

measurement: 47 potential solutions: 4 prob of best guess: 0.514291829821

measurement: 12 potential solutions: 4 prob of best guess: 0.514996238922

measurement: 16 potential solutions: 4 prob of best guess: 0.586513177343

measurement: 16 potential solutions: 4 prob of best guess: 0.654416532032

measurement: 46 potential solutions: 4 prob of best guess: 0.654606418732

measurement: 26 potential solutions: 4 prob of best guess: 0.587215897594

measurement: 42 potential solutions: 4 prob of best guess: 0.654720985541

measurement: 35 potential solutions: 4 prob of best guess: 0.587144733857

measurement: 26 potential solutions: 4 prob of best guess: 0.516254101845

measurement: 31 potential solutions: 4 prob of best guess: 0.615383141046

measurement: 30 potential solutions: 4 prob of best guess: 0.705670642705

measurement: 41 potential solutions: 4 prob of best guess: 0.761589968294

measurement: 46 potential solutions: 4 prob of best guess: 0.76176304029

measurement: 16 potential solutions: 4 prob of best guess: 0.810104284176

measurement: 16 potential solutions: 4 prob of best guess: 0.8505338556

measurement: 21 potential solutions: 4 prob of best guess: 0.85060579791

measurement: 15 potential solutions: 4 prob of best guess: 0.883616952898

measurement: 47 potential solutions: 4 prob of best guess: 0.883623423049

measurement: 40 potential solutions: 4 prob of best guess: 0.910098631886

measurement: 21 potential solutions: 4 prob of best guess: 0.910107783616

measurement: 30 potential solutions: 4 prob of best guess: 0.93821763582

measurement: 15 potential solutions: 4 prob of best guess: 0.952937925621

measurement: 25 potential solutions: 4 prob of best guess: 0.938221807153

measurement: 50 potential solutions: 4 prob of best guess: 0.910111271486

measurement: 16 potential solutions: 4 prob of best guess: 0.931033913426

measurement: 41 potential solutions: 4 prob of best guess: 0.947367831568

measurement: 20 potential solutions: 4 prob of best guess: 0.947368224558

measurement: 27 potential solutions: 4 prob of best guess: 0.93103434879

measurement: 21 potential solutions: 4 prob of best guess: 0.931034438103

measurement: 17 potential solutions: 4 prob of best guess: 0.947368397934

`etc...`

`And now we can assess the value we arrived at using bayes law with the actual value from our model:`

`estimated coin values: [1, 1, 5, 5, 10, 10, 25]`

`actual coin values: [1, 1, 5, 5, 10, 10, 25]`

estimated == actual? True