<jupyter_start><jupyter_text>GS541 phylogenetics homework 2 Estimating likelihoods for sequence modelsAt the end of this segment, I hope you understand:* How we can perform maximum likelihood branch length inference* Modeling substitution with a continuous-time Markov chain<jupyter_code>from scipy.stats import bernoulli, poisson

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import random<jupyter_output><empty_output><jupyter_text>According to the [Poisson PMF](https://en.wikipedia.org/wiki/Poisson_distributionProbability_mass_function) the probability of having **no** mutation in time $t$ is $e^{-\mu t}$. (This is why people say that there are "exponentially distributed waiting times" between Poisson-distributed events.)Let's say that we start in state 0, and when we have a mutation it randomizes the state so that it's 0 or 1 with equal probability._Exercise: fill in the formula below for the probability of being in state 1 after time t with $\mu=1$. Check to make sure it matches (approximately) with the simulations we did for the first homework._<jupyter_code>def probability_of_a_one(t):
    FILL_IN_MISSING_HERE

x_values = np.linspace(0., 2.)
ones_probability = np.array([probability_of_a_one(x) for x in x_values])
pd.DataFrame({"x": x_values, "ones_probability": ones_probability}).plot(x="x", y="ones_probability")<jupyter_output><empty_output><jupyter_text>Before we perfom branch length estimation, I remind you that branch length is defined as the expected number of **substitutions** (not mutations) per site. Refer back to the previous homework if you have forgotten the difference between these two._Exercise: classical maximum-likelihood phylogenetic inference only attempts to infer branch length, which is the product of substitution rate and time. Why is that?_   It turns out that although one can write down probabilities as above for a simple 0/1 model, for more general models one uses matrix exponentiation. In the next few exercises we'll see how that works.We'll formulate the model in terms of a [transition rate matrix](https://en.wikipedia.org/wiki/Transition_rate_matrix) or "Q matrix" for a continuous time Markov chain. This matrix has rates of transition from one state to another on the off-diagonal elements. The diagonal elements are negative of the non-diagonal row sums, which is intuitively the rate with which things are transitioning away from that state.Here's an example, in which states transition from 0 to 1 with rate 0.5, and from 1 to 0 with rate 0.5 (this 0.5 is to make things agree with our previous definition in terms of mutations appearing with rate 1):<jupyter_code>Q = 0.5 * np.matrix([[-1, 1], [1, -1]])
Q
I = np.eye(2)
I<jupyter_output><empty_output><jupyter_text>If we take $I$ and add a small multiple of $Q$, it's like taking a small discrete step as directed by $Q$.For example:<jupyter_code>I + (1/3)*Q<jupyter_output><empty_output><jupyter_text>Recall (from the previous homework) that when we take the power of probability transition matrices, it's like doing that probabilistic transition several times. One of the [definitions of the exponential function](https://en.wikipedia.org/wiki/Characterizations_of_the_exponential_function) is in terms of a limit like$$\lim_{n \rightarrow \infty} \left(1 + \frac1n\right)^n.$$We can do that with matrices too! In that case we think of taking smaller and smaller probabilistic transitions more and more frequently. In terms of DNA, we can think of substitution events becoming rarer as the number of generations gets large.In this case, our "smaller and smaller" transitions are from the matrix$$I + \frac1n Q .$$<jupyter_code>def matrix_exponential_approximation(n):
    return np.linalg.matrix_power(I + (1/n)*Q, n)<jupyter_output><empty_output><jupyter_text>_Exercise: demonstrate with a plot that the limit here has an entry converging to `probability_of_a_one(1)`, and that it agrees with what the matrix exponential implementation in SciPy gives._ The upshot of all of this is that we can model the substitution processes with a probability transition matrix of the form$$P(t) = e^{t Q}.$$As in the previous homework, the 0,0 entry of this matrix is the probability of a site being in state 0 after time $t$ given that it was in state 0 to begin with.The 1,0 entry of this matrix is the probability of the site being in state 1 after time $t$ given that it was in state 0 to begin with.We will describe in class how using a Q matrix allows for flexible modeling of DNA substitution (beyond the simple toy model we used here)._Exercise: We are going to perform likelihood-based branch-length inference using a Q matrix! Say we have a sequence that starts in state `00000` and after some time $t$ evolves to be `00101`. Making the standard phylogenetic assumptions,  estimate $t$ by plotting the likelihood function and eyeballing a peak._<jupyter_code>def likelihood(t):
    FILL_IN_MISSING_HERE

x_values = np.linspace(1., 3.)
likelihoods = np.array([likelihood(x) for x in x_values])
pd.DataFrame({"x": x_values, "likelihood": likelihoods}).plot(x="x", y="likelihood")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07198292588280947, 'fraction_numerical': 0.011253395421032208, 'mean_word_length': 5.427680798004988, 'pattern_counts': {'":': 4, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '542366', 'license': 'MIT', 'n_tokens_mistral': 1475, 'n_tokens_neox': 1384, 'n_words': 684, 'path': 'notebooks/homework_2.ipynb', 'repo_name': 'matsen/teaching', 'size': '5154'}