Sunday, February 28, 2010

Jukes-Cantor (7)

I wrote a little script that derives the eigenvalues and eigenvectors for a simple Q matrix, and then tests whether U Λ U-1 = Q, as we said last time.

It looks good, except that I haven't been able to actually do the exponentiation part yet.


import numpy as np
import math
a = 0.00333333
Q = np.array(
[[ -3*a, a, a, a],
[ a, -3*a, a, a],
[ a, a, -3*a, a],
[ a, a, a, -3*a]])

evals, evecs = np.linalg.eig(Q)
U = evecs
L = np.diag(evals)
Ui = np.linalg.inv(U)
temp = np.dot(U,L)
P = np.dot(temp,Ui)

names = ['Q','evals','evecs','U','L','Ui','P']
for i,v in enumerate([Q,evals,evecs,U,L,Ui,P]):
print names[i]
for row in v:
try:
temp = [str(round(e,3)).rjust(5) for e in row]
print ' '.join(temp)
except:
print row
print



Q
-0.01 0.003 0.003 0.003
0.003 -0.01 0.003 0.003
0.003 0.003 -0.01 0.003
0.003 0.003 0.003 -0.01

evals
-0.01333332
-1.30104260698e-18
-0.01333332
-0.01333332

evecs
-0.866 0.5 0.269 0.024
0.289 0.5 0.624 0.314
0.289 0.5 -0.182 -0.818
0.289 0.5 -0.711 0.481

U
-0.866 0.5 0.269 0.024
0.289 0.5 0.624 0.314
0.289 0.5 -0.182 -0.818
0.289 0.5 -0.711 0.481

L
-0.013 0.0 0.0 0.0
0.0 -0.0 0.0 0.0
0.0 0.0 -0.013 0.0
0.0 0.0 0.0 -0.013

Ui
-0.866 0.542 0.235 0.089
0.5 0.5 0.5 0.5
0.0 0.789 -0.102 -0.688
0.0 0.322 -0.811 0.489

P
-0.01 0.003 0.003 0.003
0.003 -0.01 0.003 0.003
0.003 0.003 -0.01 0.003
0.003 0.003 0.003 -0.01

No comments: