12 Homework on Likelihoods and MLEs - M2L5HW1

Caution

Section omitted to comply with the Honor Code

--- title: "Homework on Likelihoods and MLEs - M2L5HW1" subtitle: "Bayesian Statistics: From Concept to Data Analysis" categories: - Bayesian Statistics keywords: - Likelihoods - MLEs - Homework --- ::::: {.content-visible unless-profile="HC"} ::: {.callout-caution} Section omitted to comply with the Honor Code ::: ::::: ::::: {.content-hidden unless-profile="HC"} [Political preferences]{.column-margin} You are trying to ascertain your American colleague's *political preferences*. To do so, you design a *questionnaire* with five yes/no questions relating to current issues. The questions are all worded so that a "yes" response indicates a conservative viewpoint. [$X_i \sim Bernoulli(\theta)\quad i\in 1...5$]{.column-margin} Let $\theta$ be the unknown political viewpoint of your colleague, which we will assume can only take values 'Conservative' or 'Liberal' [$\theta \in \{\text{C}, \text{L}\}$]{.column-margin}. You have no reason to believe that your colleague leans one way or the other, so you assign the **prior** [$\mathbb{P}r(\theta = \text{C}) = 0.5$]{.column-margin}. Assume the five questions are independent and let $Y$ count the number of "yes" responses. [$Y\stackrel{iid}\sim \sum_{i=0}^5 X_i(x=y)$]{.column-margin} - If your colleague is 'conservative', then the probability of a "yes" response to any given question is 0.8. [$\mathbb{P}r(X_i=y \mid \theta=\text{C}) = 0.8$]{.column-margin} - If your colleague is 'liberal', the probability of a "no" response on any given question is 0.7.[$\mathbb{P}r(X_i = n \mid \theta=\text{L}) = 0.7$]{.column-margin} ::: {#exr-likelihood-1} [Political preferences]{.column-margin} What is an appropriate likelihood for this scenario? ::: ::: {.solution .callout-tip collapse="true"} #### Solution: $$ \begin{aligned} f(\theta \mid y) &= {5 \choose y} \cdot 0.8^y \cdot 0.2^{5-y}\cdot \mathbb{I}_{(\theta=C)} \\ &+ {5 \choose y} \cdot 0.3^y \cdot 0.7^{5-y} \cdot \mathbb{I}_{(\theta=L)} \end{aligned} $$ If your colleague is conservative, the number of "yes" responses will follow a $Binomial(5, 0.8)$. If liberal, the number of "yes" responses will follow a $Binomial(5, 0.3)$. ::: ::: {#exr-likelihood-2} [Political preferences]{.column-margin} Suppose you ask your colleague the five questions and he answers "no" to all of them. What is the MLE for $\theta$? ::: ::: {.solution .callout-tip collapse="true"} #### Solution: since all answers are the same the MLE is easy to simplify $$ \begin{aligned} {5 \choose 0} 0.8^0 0.2^{5} &\mathbb{I}_{\{\theta=C\}} + {5 \choose 0} 0.3^0 0.7^{5} &\mathbb{I}_{\{\theta=L\}} = \\ 0.2^{5} &\mathbb{I}_{\{\theta=C\}} + 0.7^{5} &\mathbb{I}_{\{\theta=L\}} = \end{aligned} $$ I implemented this expression as a function using python ```{python} #| label: C1-L05-Ex1-1 from math import comb def l5q2(y:int=0, N=5, Py1C=0.8, Pn1L=0.7, theta='C'): """ estimates the MLE of parameter theta using Args: y (int): number of yes answers. Defaults to 0. N (int, optional): number of questions in the poll. Defaults to 5. Py1C (float, optional): Probability of answering yes conditioned on conservative. Defaults to 0.8. Pn1L (float, optional): Probability of answering no conditioned on liberal. Defaults to 0.7. theta (string): 'C' if conservative 'L' if liberal. Defaults to 'C' Returns: float: the MLE estimate """ # indicator functions I_c = lambda theta : theta == 'C' I_l = lambda theta : theta == 'L' y_comp = N-y # complement for y for use in power Pn1C = 1 - Py1C Py1L = 1 - Pn1L binomial = comb(5,y) mle = \ binomial * (Py1C ** y) * (Pn1C ** y_comp) * I_c(theta)+\ binomial * (Py1L ** y) * (Pn1L ** y_comp) * I_l(theta) return round(mle,4) print(f"MLE of theta=comservative: {l5q2(y=0, theta='C')}") print(f"MLE of theta=liberal: {l5q2(y=0, theta='L')}") ``` since the expected value for $\mathbb{E}(\text{L}) = 0.17 > \mathbb{E}(\text{C})=0.00$ is greater, the MLE is: `liberal` This result is intuitive because if your colleague answered "no" to all conservative-leaning questions, he is unlikely to be a conservative. ::: ::: {#exr-likelihood-3} [Political preferences]{.column-margin} Recall that Bayes' theorem gives $$ f(\theta \mid y)=\frac{f(y \mid \theta )\color{red}{f (\theta)} }{\sum_{\theta} f(y \mid \theta )f( \theta )} $$ What is the corresponding expression for this problem? ::: ::: {.solution .callout-tip collapse="true"} #### Solution: $$ f(\theta \mid y)= \frac{ {5 \choose y} 0.8^y 0.2^{5-y} {\color{red}(0.5)} \mathbb{I}_{(\theta=C)} + {5 \choose y} 0.3^y 0.7^{5-y} {\color{red}(0.5)}\mathbb{I}_{(\theta=L) } }{ {5 \choose y} 0.8^y 0.2^{5-y} {\color{red}(0.5)} + {5 \choose y} 0.3^y 0.7^{5-y}{\color{red}(0.5)}} $$ - The prior probability was 0.5 for both values of $\theta$. - We have summed over all possibilities of $\theta$ to get the denominator. - The denominator is the **marginal probability** of observing $y$, which gives us a *normalizing constant*, which does not contains $\theta$ and evaluates to a number when we plug in y. ::: ::: {#exr-likelihood-4} [Political preferences]{.column-margin} Evaluate the expression in @exr-likelihood-3 for y=0 and report the posterior probability that your colleague is conservative, given that he responded "no" to all of the questions. ::: ::: {.solution .callout-tip collapse="true"} #### Solution: and this is code to implement Bayes' rule ```{python} #| label: C1-L05-Ex1-2 from math import comb def l5q4(y:int=0, Py1C:float=0.8, Pn1L:float=0.7, prior={'C': 0.5,'L': 0.5}, theta='C')->float: """ The posterior probability for theta, the political preference is given y 'yes' answers in a poll, using Bayes rule with: Args: y (int): number of yes answers. Defaults to 0. N (int, optional): number of questions in the poll. Defaults to 5. Py1C (float, optional): Probability of answering yes conditioned on conservative. Defaults to 0.8. Pn1L (float, optional): Probability of answering no conditioned on liberal. Defaults to 0.7. prior (dict, optional): Prior for theta the political preference. Defaults to {'C': 0.5,'L': 0.5}. theta (string): 'C' if conservative 'L' if liberal. Defaults to 'C' Returns: float: the posterior probability for theta given y. """ i_c = lambda theta : theta=='C' i_l = lambda theta : theta=='L' y_comp=5-y Pn1C=1-Py1C Py1L=1-Pn1L binomial=comb(5,y) numerator = \ binomial * (Py1C**y) * (Pn1C**y_comp) * prior['C'] * i_c (theta)+\ binomial * (Py1L**y) * (Pn1L**y_comp) * prior['L'] * i_l(theta) denominator = \ binomial*(Py1C**y)*(Pn1C**y_comp)*prior['C'] +\ binomial*(Py1L**y)*(Pn1L**y_comp)*prior['L'] posterior = numerator / denominator return round(posterior,3) print(f"{l5q4(y=0, theta='C')=}") ``` ::: ::: {#exr-likelihood-5} [Political preferences]{.column-margin} Evaluate the expression in @exr-likelihood-3 for $y=0$ and report the posterior probability that your colleague is liberal, given that he responded "no" to all of the questions. ::: ::: {.solution .callout-tip collapse="true"} #### Solution: we can reuse the function ```{python} #| label: C1-L05-Ex1-3 print(f"{l5q4(y=0, theta='L')=}") ``` ::: ::: {#exr-likelihood-6} [Loaded Coins]{.column-margin} Recall in [@exm-two-coins], your brother has a fair coin which comes up heads 50% of the time and a loaded coin which comes up heads 70% of the time. Suppose now that he has a third coin that comes up tails 70% of the time. Again, you don't know which coin your brother has brought you, so you are going to test it by flipping it 4 times, where X counts the number of heads. Let $\theta$ identify the coin so that there are three possibilities $\theta =fair$, $\theta =loaded$ favoring heads and $\theta =loaded$ favoring heads, and \theta =loaded favoring tails Suppose the prior is now $\mathbb{P}r( \theta = fair)=0.4$, $\mathbb{P}r( \theta =loaded heads)=0.3$, and loaded tails $\mathbb{P}r( \theta =loaded tails) = 0.3$. Our prior probability that the coin is loaded is still 0.6, but we do not know which loaded coin it is, so we split the probability evenly between the two options. - What is the form of the likelihood now that we have three options? ::: ::: {.solution .callout-tip collapse="true"} #### Solution: $$ \begin{aligned} f(x|\theta) =& {4 \choose x} 0.5^x 0.5^{4-x} \mathbb{I}_{\{\theta=fair\}} + \\ & {4 \choose x} 0.3^x 0.7^{4-x} \mathbb{I}_{\{\theta=loaded heads\}} +\\ & {4 \choose x} 0.7^x 0.3^{4-x} \mathbb{I}_{\{\theta=loaded tails\}} \end{aligned} $$ ::: ::: {#exr-likelihood-7} [Loaded Coins]{.column-margin} Suppose you flip the coin four times and it comes up heads twice. What is the MLE for $\theta$? ::: ::: {.solution .callout-tip collapse="true"} #### Solution: $$ {4 \choose x} 0.5^4 $$ is the highest value among the three when $X=2$ ::: ::: {#exr-likelihood-8} [Loaded Coins]{.column-margin} Suppose you flip the coin four times and it comes up heads twice. What is the posterior probability that this is a fair coin? ::: ::: {.solution .callout-tip collapse="true"} #### Solution: ```{python} #| label: C1-L05-Ex1-4 # key # 'F' for fair, 'LH' for Loaded Heads, 'LT' for Loaded Tails # PH1F is the conditional probability of heads given Fair def l5q8(x, PH1F=0.5, PH1LH=0.7, PT1LT=0.7, prior={'F':0.4,'LH':0.3,'LT':0.3}, theta='F'): # indicator functions i_F = lambda theta : theta == 'F' i_LH = lambda theta : theta == 'LH' i_LT = lambda theta : theta == 'LT' c_comp=4-x PT1F =1-PH1F PT1LH=1-PH1LH PH1LT=1-PT1LT binomial=comb(4,x) numerator=binomial*(PH1F**x)*(PT1F**c_comp)*prior['F']*i_F(theta) + \ binomial*(PH1LH**x)*(PT1LH**c_comp)*prior['LH']*i_LH(theta) +\ binomial*(PH1LT**x)*(PT1LT**c_comp)*prior['LT']*i_LT(theta) denominator = binomial*(PH1F**x)*(PT1F**c_comp)*prior['F'] + \ binomial*(PH1LH**x)*(PT1LH**c_comp)*prior['LH'] +\ binomial*(PH1LT**x)*(PT1LT**c_comp)*prior['LT'] posterior=numerator/denominator return round(posterior,2) print(f"{l5q8(x=2,theta='F')=}") ``` ::: ::: {#exr-likelihood-9} [Loaded Coins]{.column-margin} Suppose you flip the coin four times and it comes up heads twice. What is the posterior probability that this is a loaded coin (favoring either heads or tails)? Round your answer to two decimal places. ::: ::: {.solution .callout-tip collapse="true"} #### Solution: ```{python} #| label: C1-L05-Ex1-5 print(f"{1-l5q8(x=2,theta='F')=}") print(f"{l5q8(x=2,theta='LH')+l5q8(x=2,theta='LT')=}") ``` ::: ::: {.solution .callout-tip collapse="true"} #### Solution: I think that the above is pretty neat for a quick calculation but that one could do better. The conditional probabilities and the priors can be specified using a map from keys to probabilities. ```{python} #| label: C1-L05-Ex1-6 from math import comb # key # 'F' for fair, 'LH' for Loaded Heads, 'LT' for Loaded Tails # PH1F is the conditional probability of heads given Fair def l5bayes(N: int, x:int, cond: dict[str, float], prior: dict[str, float],theta:str) -> float: """ estimates the posterior probability for a Bernoulli experiment with three types of coins using Bayes' law and a binomial likelihood. Args: N (int): total trials. x (int): count of successful trials (Heads). cond (dict[str, float]): conditional probability of "H|F" expressed as a map from keys to their probability. prior (dict[str, float]): prior probability expressed as a map from keys {F,LT,LH} to their to the probability. theta (str): the value of the parameter F for fair LH and LT for loaded heads and tails. Returns: int: the posterior probability for the given value of the parameter theta """ thetas = list(prior.keys()) # we get the parameter's keys from the prior. binomial=comb(N,x) indicator = {} likelihood = {} joint_terms = {} marginal_terms = {} joint_probability,marginal = 0,0 for theta_key in thetas: indicator[theta_key] = lambda theta : theta == theta_key likelihood[theta_key] = binomial * cond[f'H|{theta_key}'] ** x * cond[f'T|{theta_key}'] ** (N-x) joint_terms[theta_key] = likelihood[theta_key] * prior[theta_key] * indicator[theta_key](theta) marginal_terms[theta_key] = likelihood[theta_key] * prior[theta_key] joint_probability += joint_terms[theta_key] marginal += marginal_terms[theta_key] posterior = joint_probability / marginal return round(posterior,3) cond_p_map: dict[str, float] = {'H|F': 0.5,'T|F': 0.5, 'H|LH': 0.7,'T|LH': 0.3, 'H|LT': 0.3,'T|LT': 0.7} prior_map: dict[str, float] = {'F': 0.4, 'LH': 0.3, 'LT': 0.3} N: int=5 print(f"{l5bayes(x=2, N=N, cond=cond_p_map, prior=prior_map, theta='F')=}") print(f"{1 - l5bayes(x=2, N=N, cond=cond_p_map, prior=prior_map, theta='F')=}") print(f"{l5bayes(x=2, N=N, cond=cond_p_map, prior=prior_map, theta='LH') + l5bayes(x=2,N=N, cond=cond_p_map, prior=prior_map, theta='LT')=}") ``` ::: I still think I could do better: - this has the binomial encoded into the likelihood, perhaps when we look at other problems it can use for the distribution and prior actual functions, possibly from SCIPY. - I like that the values of theta are collected from the prior. - the H, T events are hard coded, it would be better if they too were derived from the `cond_p_map` event-values set. :::::