<jupyter_start><jupyter_text>[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mravanba/comp551-notebooks/blob/master/MLE.ipynb) Maximum Likelihood<jupyter_code>import numpy as np
#%matplotlib notebook
%matplotlib inline
import matplotlib.pyplot as plt
from IPython.core.debugger import set_trace
np.random.seed(1234)<jupyter_output><empty_output><jupyter_text>Lets assume that a thumbtack has a probability $0 < \theta^* < 1$ of _heads_ .We use a **Bernoulli** distribution [$\mathrm{Bernoulli}(x | \theta) = \theta^{x} (1-\theta)^{1-x}$] to model this<jupyter_code>#Function to compute parametric probability mass function
#If you pass arrays or broadcastable arrays it computes the elementwise bernoulli probability mass
Bernoulli = lambda theta,x: theta**x * (1-theta)**(1-x)
theta_star = .4
Bernoulli(theta_star, 0)<jupyter_output><empty_output><jupyter_text>The **likelihood** function for this dataset is$$L(\theta) = \prod_{n=1}^{10} \mathrm{Bernoulli}(\theta, x^{(n)})$$Note that this is a function of $\theta$, it is not a probability distribution, so it doesn't sum to 1.Next we plot this function for $n$ randomly generated samples using $\theta^*$.<jupyter_code>n = 10                                                              #number of random samples you want to consider 
xn = np.random.rand(n) < theta_star                                 #Generates n element boolean array where elements are True with probability theta_star and otherwise False
xn = xn.astype(int)                                                 #to change the boolean array to intergers [0:False, 1:True]
print("observation {}".format(xn))
#Function to compute the log likelihood
#Note that you can either pass this function a scalar(always broadcastable) theta or a broadcastable(in data axis) theta to get likelihood value or values
#Also note that we added an extra dimension in xn to broadcast it along theta dimension
L = lambda theta: np.prod(Bernoulli(theta, xn[None,:]), axis=-1)
#we generate 100 evenly placed values of theta from 0 to 1
theta_vals = np.linspace(0,1,100)[:, None]                         #Note that we made the array broadcastable by adding an extra dimension for data
plt.plot(theta_vals, L(theta_vals), '-')
plt.xlabel(r"$\theta$")
plt.ylabel(r"Likelihood $L(\theta)$")
#plt.title(r"likelihood function Bernoulli("+str(xn.astype(int))+r"$|\theta)$")
plt.show()<jupyter_output>observation [1 0 0 0 0 1 1 0 0 0]<jupyter_text>**maximum likelihood** method for learning a probabilistic model picks $\theta$ that maximizes $L(\theta)$ for a given set of observations $\mathcal{D} = \{x^{(1)}, \ldots, x^{(N)}\}$. For _Bernoulli_ dist. max. likelihood estimate is the $\widehat{\theta} = \frac{1}{N} \sum_n x^{(n)}$, corresponding to the peak of the likelihood plot. The value of the likelihood shrinks exponentially as we increase the number of observations $N$ -- try `N=100`. A more numerically stable value with the same _maximum_ is the **log-likelihood** $\ell(\theta) = \log(L(\theta))$. It is also customary to minimize the _negative log-likelihood_ (NLL).Let's plot NLL for different $N$ -- as we increase our data-points the ML estimate often gets better.<jupyter_code>#Generates 2^12 element boolean array where elements are True with probability theta_star and otherwise False
xn_max = np.random.rand(2**12) < theta_star
for r in range(1,6):
    n = 4**r              #number of data samples for r-th iteration
    xn = xn_max[:n]       #slice them from the total samples generated
    #Function to compute the log likelihood (Implementation exactly similar to the likelihood function)
    ll = lambda theta: np.sum(np.log(Bernoulli(theta, xn[None,:])), axis=-1)
    theta_vals = np.linspace(.01,.99,100)[:, None]
    ll_vals = -ll(theta_vals)
    #Plot the log likelihood values
    plt.plot(theta_vals, ll_vals, label="n="+str(n))
    max_ind = np.argmin(ll_vals)                            #Stores the theta corresponding to minimum log likelihood
    plt.plot(theta_vals[max_ind], ll_vals[max_ind], '*')
#to get the horizontal line for theta
plt.plot([theta_star,theta_star], [0,ll_vals.max()], '--', label=r"$\theta^*$")
plt.xlabel(r"$\theta$")
plt.ylabel(r"Negative Log-Likelihood $-\ell(\theta)$")
plt.yscale("log")
plt.title("ML solution with increasing data")
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12372535690006799, 'fraction_numerical': 0.015182415590301382, 'mean_word_length': 4.132558139534884, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 16, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '63144', 'license': 'MIT', 'n_tokens_mistral': 1429, 'n_tokens_neox': 1328, 'n_words': 457, 'path': 'MLE.ipynb', 'repo_name': 'mravanba/comp551-notebooks', 'size': '4413'}