This lab is built around the process of identifying the fault with a coffee machine.

Your task is to:

Given the structure of a graphical model for the state of a coffee machine learn the distributions from data.
Implement belief propagation, so you can evaluate the probability of each failure given the available evidence.
Identify the most probable failure for a set of broken coffee machines.

Marking and Submission

These lab exercises are marked, and contribute to your final grade. For this lab exercise there are 3 places where you are expected to enter your own code, for 15 marks overall. Every place you have to add code is indicated by

# **************************************************************** n marks

with instructions above the code block.

Please submit your completed workbook using Moodle before 5pm on the 28th November 2018. The workbook you submit must be an .ipynb file, which is saved into the directory you’re running Jupyter; alternatively you can download it from the menu above using File -> Download As -> Notebook (.ipynb). Remember to save your work regularly (Save and checkpoint in the File menu, the icon of a floppy disk, or Ctrl-S); the version you submit should have all code blocks showing the results (if any) of execution below them. You will normally receive feedback within a week.

%matplotlib inline

Coffee Machine Data Set

You can download from moodle a zip file that contains the data as a csv file. The below code will load the data directly from the zip file, so you don’t have to unzip it.

Each row of the file contains a fully observed coffee machine, with the state of every random variable. The random variables are all binary, with False represented by 0 and True represented by 1. The variables are:

Failures (you’re trying to detect these):

1. he – No electricity
1. fp – Fried power supply unit
1. fc – Fried circuit board
1. wr – Water reservoir empty
1. gs – Group head gasket forms seal
1. dp – Dead pump
1. fh – Fried heating element

Mechanism (these are unobservable):

1. pw – Power supply unit works
1. cb – Circuit board works
1. gw – Get water out of group head

Diagnostic (these are the tests the mechanic can run – observable):

1. ls – Room lights switch on
1. vp – A voltage is measured across power supply unit
1. lo – Power light switches on
1. wv – Water visible in reservoir
1. hp – Can hear pump
1. me – Makes espresso
1. ta – Makes a hot, tasty espresso

For the above the number is the column number of the provided data (dm) and the two letter code a suggested variable name.

For if you are unfamiliar with an espresso coffee machine here is a brief description of how one works (you can ignore this if you want):

The user puts ground coffee into a portafilter (round container with a handle and two spouts at the bottom), tamps it (compacts the coffee down), and clamps the portafilter into the group head at the front of the machine. A gasket (rubber ring) forms a seal between the portafilter and group head. A button is pressed. Water is drawn from a reservoir by a pump into a boiler. In the boiler a heating element raises the waters temperature, before the pump pushes it through the group head and into the portafilter at high pressure. The water passes through the coffee grinds and makes a tasty espresso.

The graphical model showing how the variables are related is included on moodle as coffee machine.pdf; here it is given as conditional probabilities:

Failures:

P_he: P(𝚗𝚘 𝚎𝚕𝚎𝚌𝚝𝚛𝚒𝚌𝚒𝚝𝚢)
P_fp: P(𝚏𝚛𝚒𝚎𝚍 𝚙𝚜𝚞)
P_fc: P(𝚏𝚛𝚒𝚎𝚍 𝚌𝚒𝚛𝚌𝚞𝚒𝚝 𝚋𝚘𝚊𝚛𝚍)
P_wr: P(𝚠𝚊𝚝𝚎𝚛 𝚛𝚎𝚜𝚎𝚛𝚟𝚘𝚒𝚛 𝚎𝚖𝚙𝚝𝚢)
P_gs: P(𝚐𝚛𝚘𝚞𝚙 𝚑𝚎𝚊𝚍 𝚐𝚊𝚜𝚔𝚎𝚝 𝚜𝚎𝚊𝚕 𝚋𝚛𝚘𝚔𝚎𝚗)
P_dp: P(𝚍𝚎𝚊𝚍 𝚙𝚞𝚖𝚙)
P_fh: P(𝚏𝚛𝚒𝚎𝚍 𝚑𝚎𝚊𝚝𝚒𝚗𝚐 𝚎𝚕𝚎𝚖𝚎𝚗𝚝)

Mechanism:

P_pw_he_fp: P(𝚙𝚜𝚞 𝚠𝚘𝚛𝚔𝚜|𝚗𝚘 𝚎𝚕𝚎𝚌𝚝𝚛𝚒𝚌𝚒𝚝𝚢,𝚏𝚛𝚒𝚎𝚍 𝚙𝚜𝚞)
P_cb_pw_fc: P(𝚌𝚒𝚛𝚌𝚞𝚒𝚝 𝚋𝚘𝚊𝚛𝚍 𝚠𝚘𝚛𝚔𝚜|𝚙𝚜𝚞 𝚠𝚘𝚛𝚔𝚜,𝚏𝚛𝚒𝚎𝚍 𝚌𝚒𝚛𝚌𝚞𝚒𝚝 𝚋𝚘𝚊𝚛𝚍)
P_gw_cb_wr_dp: P(𝚐𝚎𝚝 𝚠𝚊𝚝𝚎𝚛|𝚌𝚒𝚛𝚌𝚞𝚒𝚝 𝚋𝚘𝚊𝚛𝚍 𝚠𝚘𝚛𝚔𝚜,𝚠𝚊𝚝𝚎𝚛 𝚛𝚎𝚜𝚎𝚛𝚟𝚘𝚒𝚛 𝚎𝚖𝚙𝚝𝚢,𝚍𝚎𝚊𝚍 𝚙𝚞𝚖𝚙)

Diagnostic:

P_ls_he: P(𝚕𝚒𝚐𝚑𝚝𝚜 𝚜𝚠𝚒𝚝𝚌𝚑 𝚘𝚗|𝚗𝚘 𝚎𝚕𝚎𝚌𝚝𝚛𝚒𝚌𝚒𝚝𝚢)
P_vp_pw: P(𝚟𝚘𝚕𝚝𝚊𝚐𝚎 𝚊𝚌𝚛𝚘𝚜𝚜 𝚙𝚜𝚞|𝚙𝚜𝚞 𝚠𝚘𝚛𝚔𝚜)
P_lo_cb: P(𝚙𝚘𝚠𝚎𝚛 𝚕𝚒𝚐𝚑𝚝 𝚘𝚗|𝚌𝚒𝚛𝚌𝚞𝚒𝚝 𝚋𝚘𝚊𝚛𝚍 𝚠𝚘𝚛𝚔𝚜)
P_wv_wr: P(𝚠𝚊𝚝𝚎𝚛 𝚟𝚒𝚜𝚒𝚋𝚕𝚎|𝚠𝚊𝚝𝚎𝚛 𝚛𝚎𝚜𝚎𝚛𝚟𝚘𝚒𝚛 𝚎𝚖𝚙𝚝𝚢)
P_hp_dp: P(𝚌𝚊𝚗 𝚑𝚎𝚊𝚛 𝚙𝚞𝚖𝚙|𝚍𝚎𝚊𝚍 𝚙𝚞𝚖𝚙)
P_me_gw_gs: P(𝚖𝚊𝚔𝚎𝚜 𝚎𝚜𝚙𝚛𝚎𝚜𝚜𝚘|𝚐𝚎𝚝 𝚠𝚊𝚝𝚎𝚛,𝚐𝚛𝚘𝚞𝚙 𝚑𝚎𝚊𝚍 𝚐𝚊𝚜𝚔𝚎𝚝 𝚜𝚎𝚊𝚕 𝚋𝚛𝚘𝚔𝚎𝚗)
P_ta_me_fh: P(𝚝𝚊𝚜𝚝𝚢|𝚖𝚊𝚔𝚎𝚜 𝚎𝚜𝚙𝚛𝚎𝚜𝚜𝚘,𝚏𝚛𝚒𝚎𝚍 𝚑𝚎𝚊𝚝𝚒𝚗𝚐 𝚎𝚕𝚎𝚖𝚎𝚗𝚝)

Note that while the model is close to what you may guess the probabilities are not absolute, to account for mistakes and unknown failures. For instance, the mechanic may make a mistake while brewing an espresso and erroneously conclude that the machine is broken when it is in fact awesome. The probabilities associated with each failure are not uniform.

# It may prove helpful to have a mapping between the suggested variable names and

Data: 262144 exemplars, 17 features
     Broken machines  = 122603
     Working machines = 139541

​

1. Learn Model

Below a set of variables to represent conditional probability distributions have been defined. They are a Bernoulli trial for each combination of conditional variables, given as P(𝙵𝚊𝚕𝚜𝚎|...) in [0,...] and P(𝚃𝚛𝚞𝚎|...) in [1,...] (It may be easier to think of them as boring categorical distributions).

To fit a maximum likelihood model you should first use them as counters – loop over the data set and count how many of each combination exist. You must then normalise them so that the sum over axis 0 is always 1. There is an extra mark for doing the right thing and including a prior. You may want to know that the conjugate prior to the Bernoulli trial represented as [P(𝙵𝚊𝚕𝚜𝚎),P(𝚃𝚛𝚞𝚎)] is a Beta distribution; a uniform prior would be a reasonable choice (it can be argued that the expected failure probability is low, and you should adjust the first 7 variables accordingly, but given the quantity of data available it’s not going to matter).

Hint:

The use of 0=False and 1=True both in the dm table and in the conditional probability distributions is very deliberate.
Consider putting all of the variables into a list with extra information about them (such as indices from nti) to make your code neater.
If you make a mistake you could easily end up with NaN or infinity – which would break the next step. Print them out so you can check they are valid!
Do not write unique code for each variable – that will be hundreds of lines of code and very tedious. It’s possible to get all of the marks in 7 lines of code, if you’re particularly sneaky.

(4 marks)

2 marks for counting all conditions
1 mark for normalising distributions
1 mark for including a sensible prior

# A set of variables that will ultimately represent conditional probability distributions.

​

Brute Force

The below code does exactly what you want to implement for step 2, but it brute forces it. Slow and memory inefficient of course, but lets you test it.

​

def brute_marginals(known):

[[0.99766698 0.00233302]
 [0.30730747 0.69269253]
 [0.80499257 0.19500743]
 [0.72705735 0.27294265]
 [0.99285305 0.00714695]
 [0.79354149 0.20645851]
 [0.99776075 0.00223925]
 [0.69478544 0.30521456]
 [0.88219973 0.11780027]
 [0.97686761 0.02313239]
 [0.91244598 0.08755402]
 [0.66892206 0.33107794]
 [0.87904112 0.12095888]
 [0.73941899 0.26058101]
 [0.72389638 0.27610362]
 [0.99157774 0.00842226]
 [0.96820919 0.03179081]]
[[0.99766699 0.00233301]
 [0.307313   0.692687  ]
 [0.80499429 0.19500571]
 [0.72711024 0.27288976]
 [0.99283907 0.00716093]
 [0.79349557 0.20650443]
 [0.99768722 0.00231278]
 [0.69477988 0.30522012]
 [0.88219196 0.11780804]
 [0.9770849  0.0229151 ]
 [0.91244598 0.08755402]
 [0.66891672 0.33108328]
 [0.87903337 0.12096663]
 [0.73940003 0.26059997]
 [0.72391236 0.27608764]
 [0.99373317 0.00626683]
 [1.         0.        ]]

2. Belief Propagation

Your task is to write a function that take known states and calculates the marginal probability for every state of the graphical model. This will require using the sum-product algorithm and passing all messages on the graph.

Here is the wikipedia page for reference: https://en.wikipedia.org/wiki/Belief_propagation (not the greatest, but not horrible)

Hints:

The order of the variables above is such that each variable is dependent only on variables that occur before it. This should save you some hassle in terms of calculating a message passing order.
You may want to add code before the function to prepare. This might involve creating a list of edges in the graphical model, in the order they need to be processed for instance. It can also be done with dictionaries or classes – choose whatever you are comfortable with, but spend time thinking it through first. If you get it wrong your code will get messy!
A good approach to recording messages is a dictionary indexed [from, to]
The easiest way to include a known value in the model is often to add/replace the unary term for that value.
This problem is small enough that you can brute force it – code to do so has been provided above. Do test that your implementation is correct by comparing them.

(10 marks)

3 marks for sensible data structures
5 marks for sending the messages correctly
2 marks for correctly calculating the marginals at the end

def marginals(known):

3. What’s broken?

Simply print out what the most probable broken part of each of the below machines is. Please print it as dictionary key: name of broken part, e.g. A: wr. This will allow automarking to be used.

If you can’t get your belief propagation implementation to work feel free to use brute_marginals instead.

(1 mark)

repair = {}

​