In [1]:

    

# set some ipython notebook properties
%matplotlib inline

# set degree of verbosity (adapt to INFO for more verbose output)
import logging
logging.basicConfig(level=logging.WARNING)

# set figure sizes
import pylab
pylab.rcParams['figure.figsize'] = (10.0, 8.0)

# set display width for pandas data frames
import pandas as pd
pd.set_option('display.width', 1000)


    

In [27]:

    

import numpy as np
K = np.array([[1., 2, 3],[2,1,5],[3,5,1]])
S,U = np.linalg.eigh(K)
print S
print U


    

    




[-4.12839474 -0.71069159  7.83908632]
[[ 0.17830987  0.86914323  0.46129777]
 [ 0.65036916 -0.45590689  0.60759268]
 [-0.7383939  -0.19167408  0.64655665]]

In [14]:

    

np.eye(3) * S


    

    
Out[14]:





array([[-4.12839474, -0.        ,  0.        ],
       [-0.        , -0.71069159,  0.        ],
       [-0.        , -0.        ,  7.83908632]])

In [17]:

    

U


    

    
Out[17]:





array([[ 0.17830987,  0.86914323,  0.46129777],
       [ 0.65036916, -0.45590689,  0.60759268],
       [-0.7383939 , -0.19167408,  0.64655665]])

In [16]:

    

np.dot(np.dot(U,np.eye(3) * S),U.T)


    

    
Out[16]:





array([[ 1.,  2.,  3.],
       [ 2.,  1.,  5.],
       [ 3.,  5.,  1.]])

In [21]:

    

print np.linalg.inv(K)
print np.dot(np.dot(U,np.eye(3) / S),U.T)


    

    




[[-1.04347826  0.56521739  0.30434783]
 [ 0.56521739 -0.34782609  0.04347826]
 [ 0.30434783  0.04347826 -0.13043478]]
[[-1.04347826  0.56521739  0.30434783]
 [ 0.56521739 -0.34782609  0.04347826]
 [ 0.30434783  0.04347826 -0.13043478]]

In [26]:

    

V = .2 * K + .8 * np.eye(3)
print V
print np.linalg.inv(V)
print .2 * np.dot(np.dot(U,np.eye(3) / S),U.T) + .8 * np.eye(3)


    

    




[[ 1.   0.4  0.6]
 [ 0.4  1.   1. ]
 [ 0.6  1.   1. ]]
[[  3.77475828e-15  -5.00000000e+00   5.00000000e+00]
 [ -5.00000000e+00  -1.60000000e+01   1.90000000e+01]
 [  5.00000000e+00   1.90000000e+01  -2.10000000e+01]]
[[ 0.59130435  0.11304348  0.06086957]
 [ 0.11304348  0.73043478  0.00869565]
 [ 0.06086957  0.00869565  0.77391304]]

In [33]:

    

print np.linalg.inv(13 * K)
print 1./13 * np.linalg.inv(K)


    

    




[[-0.08026756  0.04347826  0.02341137]
 [ 0.04347826 -0.02675585  0.00334448]
 [ 0.02341137  0.00334448 -0.01003344]]
[[-0.08026756  0.04347826  0.02341137]
 [ 0.04347826 -0.02675585  0.00334448]
 [ 0.02341137  0.00334448 -0.01003344]]

In [38]:

    

Kn = K/factor
print Kn


    

    




[[ 0.33333333  0.66666667  1.        ]
 [ 0.66666667  0.33333333  1.66666667]
 [ 1.          1.66666667  0.33333333]]

In [39]:

    

np.std(Kn)


    

    
Out[39]:





0.49965694678637196

In [7]:

    

s,u = np.linalg.eigh(np.eye(3))
print s
print u


    

    




[ 1.  1.  1.]
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]

In [8]:

    

np.dot(np.dot(u,np.eye(3) * s),u.T)


    

    
Out[8]:





array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

In [15]:

    

u2 = 1
print np.dot(u2,np.eye(3)*s)


    

    




[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]

In [27]:

    

#Why can we assume that h2 and e2 add to one?
#Because y is unit standardized then its s.d. is 1
# and K's diagional is N

covar = np.array([[float(num)] for num in xrange(490)])
np.random.seed(0)
pheno = covar * 2.0 + 100 + np.random.normal(size=len(covar))*10
print np.mean(pheno), np.std(pheno)


    

    




588.713280344 283.075919246

In [32]:

    

import pylab
pylab.plot(covar, pheno,".")
pylab.show()


    

In [ ]: