.md

Question

Consider a simple model for whether a person has the flu or not. Let F=1 indicate that a person has the flu and F=0 indicate that they don't have the flu. Let C=1 indicate that the person has a cough and C=0 indicate that they don't have a cough. Let M=1 indicate that the person has muscle pain and M=0 indicate that they don't have muscle pain. Assume that C and M are conditionally independent given F so that the probability model is P(C=c,M=m,F=f)=P(C=c|F=f)P(M=m|F=f)P(F=f).

Suppose that we ask two different doctors to supply probabilities for this model and we obtain the following results: Doctor 1: P(F=1)=0.4 P(C=1|F=0)=0.2, P(C=1|F=1)=0.8 P(M=1|F=0)=0.3, P(M=1|F=1)=0.9 Doctor 2: P(F=1)=0.5 P(C=1|F=0)=0.1, P(C=1|F=1)=0.7 P(M=1|F=0)=0.2, P(M=1|F=1)=0.8

Suppose we also have access to 10 patient records recording both the symptoms and whether each patient was diagnosed with flu or not as shown below:

Patient C M F 1 1 0 0 2 1 1 1 3 1 1 1 4 0 1 0 5 0 1 1 6 1 0 1 7 0 0 0 8 1 0 0 9 0 0 0 10 1 1 1

What is the likelihood of the data under the model specified by doctor 1? Give your answer in scientific notation (ie: 1.2345E-10). Hint: you may want to write a short program to compute this value.

Work

I started by doing some math to simplify shit

let D = patient records
let D_i_j =  ith patient record, where i∈ {1, 2, 3, ... 10}
let D_i_j =  j attribute of ith patient record, where i ∈ {1, 2, 3, ... 10}, j ∈ {c, m, f}
let M_i = model specified by doctor i, where i ∈ {1, 2}
let X = a patient record

P(D|M_1)
=Product{P(X=D_i|M_1), for all i}

P(X=D_i|M_1)
=P(F=D_i_f, C=D_i_c, M=D_i_m,|M_1)
=P(C=D_i_c|F=D_i_f)P(M=D_i_m|F=D_i_f)P(F=D_i_f,)

So I wanna compute P(D|M_1)...

I know how to do this, but I wonder if I can write it elegantly.

from operator import mul

product = lambda l: reduce(mul, l, 1)


def prob_data_given_model(data, model):
    '''
    let D = `data`
    let M_i = `model`, for some i ∈ {1, 2}
    
    Returns P(D|M_i) = Product{P(X=D_i|M_1), for all i}
    '''
    return product(prob_record_given_model(data[i], model))

def prob_record_given_model(record, model):
    '''
    let D_i = `record`, for some i ∈ {1, 2, 3, ... 10} 
    let M_i = `model`, for some i ∈ {1, 2}
    
    Returns P(X=D_i|M_1)=P(C=D_i_c|F=D_i_f)P(M=D_i_m|F=D_i_f)P(F=D_i_f)
    '''
    c, f, m = record
    return model.c_giv_f(c, f) * model.m_giv_f(m, f) * model.f(f)


class Model(object):
    def __init__(self, f_1, c_1_giv_f_0, c_1_giv_f_1, m_1_giv_f_0, m_1_giv_f_1):
        self.f_1 = f_1
        self.c_1_giv_f_0 = c_1_giv_f_0
        self.c_1_giv_f_1 = c_1_giv_f_1
        self.m_1_giv_f_0 = m_1_giv_f_0
        self.m_1_giv_f_1 = m_1_giv_f_1
        
    def c_giv_f(self, c, f):
        if f == 0:
            if c == 1:
                return self.c_1_giv_f_0
            else:
                return 1 - self.c_1_giv_f_0
		# wow this is boring and dumb, I'm stopping now        

doctor_1 = Model(0.4, 0.2, 0.8, 0.3, 0.9)
doctor_2 = Model(0.5, 0.1, 0.7, 0.2, 0.8)