## 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 npimport matha = 0.00333333Q = 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 = evecsL = 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.0030.003 -0.01 0.003 0.0030.003 0.003 -0.01 0.0030.003 0.003 0.003 -0.01evals-0.01333332-1.30104260698e-18-0.01333332-0.01333332evecs-0.866 0.5 0.269 0.0240.289 0.5 0.624 0.3140.289 0.5 -0.182 -0.8180.289 0.5 -0.711 0.481U-0.866 0.5 0.269 0.0240.289 0.5 0.624 0.3140.289 0.5 -0.182 -0.8180.289 0.5 -0.711 0.481L-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.013Ui-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.489P-0.01 0.003 0.003 0.0030.003 -0.01 0.003 0.0030.003 0.003 -0.01 0.0030.003 0.003 0.003 -0.01`