• A script to simulate sequence evolution (post)

• A peek at the answers for P

_{XX}(t) and P

_{XY}(t) under Jukes-Cantor (post)

• A simple derivation from Higgs & Attwood (post)

• Use to get a relationship between observed changes and distance (post)

• A script that sets up a matrix form of P(t) and cranks it a bit showing the convergence to the stationary distribution (post)

• Going from the solution to the differential equations of JC69 (post)

• A first mention of the equation P(t) = e

^{Q(t)}from Yang's book. (post)

• A test of matrix decomposition into eigenvalues and eigenvectors (post)

I think I've made some progress in understanding models of sequence evolution (though how well I've conveyed this, you'll be the judge). But there is still a way to go yet. I need to understand: (i) why this exponentiation thing seems to be important, (ii) use of the model for likelihood calculations (L ≈ P(data | model), (iii) extension to more realistic models like Kimura's 2-parameter model that treats transitions and transversions differently, and (iv) how to do practical calculations under the models.

Today, I want to do just a bit more with matrix exponentiation. Before, we talked about the idea that you can decompose a matrix (at least a square well-behaved matrix like ours) into something like this

`A = U Λ U`^{-1}

where, critically, the `Λ`

matrix is a diagonal matrix so that for example, its square is:This is supposed to be extensible to exponentiation (or any function), and I was having trouble wrapping my head around that. But then I got reminded that e

^{x}can be written as:

So, if we can do

`A`^{2}

, we can do `e`^{A}

. I'm at my office computer this AM, and to do what follows I needed SciPy installed. So I built and installed Python 2.6, used an available install of the gfortran compiler and FFTW, ignored the sparse matrix stuff, built and installed Numpy and SciPy. All as described here, here and here.

The script (full listing at the end of the post), has two parts. First, we do the decomposition, following the same line as previously, checking that A

^{n}by the eigenvalue method matches what we get for repeated matrix multiplication.

In the second part, I do exponentiation using the eigenvalue method, and check it using three methods from SciPy (reference).

It seems to work fine!

Output:

Listing:

## No comments:

Post a Comment