

In [1]:

    
import randomCorr

c = randomCorr.mkRandCorr(40, 10**-3)

Exemplo de matriz de correlação aleatória

O método utilizado é o descrito no artigo:

http://www.plosone.org/article/info:doi/10.1371/journal.pone.0048902



In [2]:

    
pcolor(c)









    Out[2]:





<matplotlib.collections.PolyCollection at 0x2a84650>



In [4]:

    
plot(np.linalg.eigvals(c))









    Out[4]:





[<matplotlib.lines.Line2D at 0x4172550>]

Método usual para gerar matrizes de correlação randômicas



In [3]:

    
def rand_ortho_matrix(eigen_values):
    size = len(eigen_values)
    A = np.mat(np.random.random((size,size)))
    Q, R = np.linalg.qr(A)
    
    return Q * np.diag(eigen_values) * Q.T

n_traits = 40
rand_eigen = np.array(np.random.random(n_traits))
p = rand_ortho_matrix(rand_eigen)

pcolor(numpy.array(p))









    Out[3]:





<matplotlib.collections.PolyCollection at 0x3df0050>



In [2]:

    
def CalcR2(Matrix):
   tr = Matrix.shape[1]
   x, y = np.asarray(np.invert(np.tri(tr, tr, dtype=bool)), dtype=float).nonzero()
   R2Tot = np.mean(Matrix[x, y] * Matrix[x, y])
   return R2Tot



In [3]:

    
def flexibility(matrix1, num_vectors=1000):
    traits = matrix1.shape[0]
    rand_vec = np.random.multivariate_normal(np.zeros(traits),
                                             np.identity(traits, float),
                                             num_vectors).T
    
    rand_vec = rand_vec/np.sqrt((rand_vec*rand_vec).sum(0))
    
    delta_z1 = np.dot(matrix1, rand_vec)

    ndelta_z1 = delta_z1/np.sqrt((delta_z1*delta_z1).sum(0))

    return np.mean(np.diag(np.dot(ndelta_z1.T, rand_vec)))



In [76]:

    
s = 1000
r2 = np.zeros(s)
flex = np.zeros(s)
isoc = np.zeros(s)

for i in xrange(s):
    c = randomCorr.mkRandCorr(40, 10**-3)  
    r2[i] = CalcR2(c)
    flex[i] = flexibility(c)
    isoc[i] = np.abs(np.dot(eig(c)[1][:,0], iso))

Histograma correlação com vetor isométrico



In [77]:

    
hist(isoc)









    Out[77]:





(array([277, 239, 189, 138,  85,  39,  21,   7,   2,   3]),
 array([  3.80143281e-04,   5.41018694e-02,   1.07823596e-01,
         1.61545322e-01,   2.15267048e-01,   2.68988774e-01,
         3.22710500e-01,   3.76432226e-01,   4.30153953e-01,
         4.83875679e-01,   5.37597405e-01]),
 <a list of 10 Patch objects>)

Histograma R2



In [78]:

    
hist(r2)









    Out[78]:





(array([ 15,  99, 203, 277, 212, 113,  51,  21,   7,   2]),
 array([ 0.15437342,  0.16973832,  0.18510322,  0.20046811,  0.21583301,
        0.23119791,  0.2465628 ,  0.2619277 ,  0.2772926 ,  0.2926575 ,
        0.30802239]),
 <a list of 10 Patch objects>)

Histograma da flexibilidade



In [79]:

    
hist(flex)









    Out[79]:





(array([  3,  20,  54, 126, 210, 267, 174,  97,  42,   7]),
 array([ 0.28946314,  0.29843751,  0.30741188,  0.31638625,  0.32536062,
        0.33433499,  0.34330937,  0.35228374,  0.36125811,  0.37023248,
        0.37920685]),
 <a list of 10 Patch objects>)



In [81]:

    
scatter(r2,flex)









    Out[81]:





<matplotlib.collections.PathCollection at 0x77c0290>



In [43]:

    
corrcoef(r2, flex)









    Out[43]:





array([[ 1.        , -0.94217028],
       [-0.94217028,  1.        ]])



In [65]:

    
p = np.dot(c, eig(c)[1][:,0])



In [66]:

    
np.sqrt(np.dot(p, p))









    Out[66]:





14.85366928000486



In [64]:

    
eigvals(c)









    Out[64]:





array([  1.48536693e+01,   9.00525414e+00,   6.62282829e+00,
         3.48049664e+00,   1.91936855e+00,   1.55731768e+00,
         1.10373711e+00,   5.31852658e-01,   4.31979031e-01,
         1.63227311e-01,   1.39492590e-01,   7.87821502e-02,
         6.57679478e-02,   2.93080476e-02,   7.58424751e-03,
         4.27840324e-03,   1.65061907e-03,   7.47785040e-04,
         5.04088158e-04,   3.14290933e-04,   3.08745768e-04,
         2.71951107e-04,   2.34327079e-04,   1.89527289e-04,
         1.74323924e-04,   1.61352610e-04,   1.20677024e-04,
         1.11145417e-04,   8.05011263e-05,   7.09764560e-05,
         4.97751631e-05,   2.69545672e-05,   1.72473102e-05,
         1.31837392e-05,   5.46194444e-06,   2.14067989e-06,
         5.49471321e-07,   2.51270615e-07,   4.25199418e-08,
         2.57099037e-10])



In [22]:

    
ms =[m1,m2,m3,m4]



In [71]:

    
params = []

for m in ms:
    mb = mkB(m)
    mp = calc_params(mb)
    params.append((mb, mp))



In [79]:

    
plot(params[0][1][:,0])









    Out[79]:





[<matplotlib.lines.Line2D at 0x8852710>]



In [80]:

    
plot(params[1][1][:,0])









    Out[80]:





[<matplotlib.lines.Line2D at 0x87ffbd0>]



In [81]:

    
plot(params[2][1][:,0])









    Out[81]:





[<matplotlib.lines.Line2D at 0x8d080d0>]



In [82]:

    
plot(params[3][1][:,0])









    Out[82]:





[<matplotlib.lines.Line2D at 0x9003590>]



In [102]:

    
diff = (m4-m3)/100
diff









    Out[102]:





array([[  0.00000000e+00,  -6.22541420e-04,   8.67440500e-05, ...,
         -2.62268036e-03,   1.36925238e-03,  -5.29693630e-04],
       [ -6.22541420e-04,   0.00000000e+00,  -7.83092710e-04, ...,
         -3.65399300e-05,   1.39710180e-03,  -1.12024273e-03],
       [  8.67440500e-05,  -7.83092710e-04,   0.00000000e+00, ...,
          9.54287040e-04,   4.27959982e-03,   6.99147560e-04],
       ..., 
       [ -2.62268036e-03,  -3.65399300e-05,   9.54287040e-04, ...,
          0.00000000e+00,  -4.75033010e-04,   2.14057603e-03],
       [  1.36925238e-03,   1.39710180e-03,   4.27959982e-03, ...,
         -4.75033010e-04,   0.00000000e+00,   2.15962038e-03],
       [ -5.29693630e-04,  -1.12024273e-03,   6.99147560e-04, ...,
          2.14057603e-03,   2.15962038e-03,   0.00000000e+00]])



In [103]:

    
m_min = m3
for i in xrange(101):
    pcolor(m_min)
    savefig('{}.png'.format(i))
    m_min += diff



In [56]:

    
import randomCorr as rc



In [68]:

    
m1 = np.loadtxt(open("../matrices/1-gorilla.csv","rb"),delimiter=",")
m2 = np.loadtxt(open("../matrices/2-peramelidae.csv","rb"),delimiter=",")
m3 = np.loadtxt(open("../matrices/3-molossidae.csv","rb"),delimiter=",")
m4 = np.loadtxt(open("../matrices/4-hyllostomidae.csv","rb"),delimiter=",")
#m5 = np.loadtxt(open("../matrices/5.csv","rb"),delimiter=",")
#m6 = np.loadtxt(open("../matrices/6.csv","rb"),delimiter=",")

ms =[m1,m2,m3,m4]



In [70]:

    
bp = []

for i in xrange(4):
    b_n = rc.triang_decomp(ms[i])
    p_n = rc.calc_params(b_n)
    bp.append((b_n, p_n))



In [59]:

    
for i in xrange(6):
    print np.max(ms[i] - np.dot(bp[i][0], bp[i][0].T))









    



4.4408920985e-16
5.55111512313e-16
7.77156117238e-16
6.66133814775e-16
3.46944695195e-18
0.238453007



In [45]:

    
for i in xrange(4):
    bp_i = rc.triang_from_params(bp[i][1])
    print np.max(ms[i] - np.dot(bp_i, bp_i.T))









    



9.11302650008e-09
5.11961184557e-09
7.61498694657e-09
9.09617700606e-09



In [58]:

    
diff = (bp[1][1] - bp[0][1])/100

p = bp[0][1]
m0 = ms[0]

for i in xrange(101):
    new_b = rc.triang_from_params(p)
    m0 = np.dot(new_b, new_b.T)
    pcolor(m0)
    savefig('{}.png'.format(i))
    p += diff



In [72]:

    
diff = (ms[1] - ms[0])/100

m0 = ms[0]

for i in xrange(101):
    pcolor(m0)
    savefig('m-{:03d}.png'.format(i))
    m0 += diff



In [39]:

    
pcolor(m1)









    Out[39]:





<matplotlib.collections.PolyCollection at 0xe930510>



In [23]:

    
pcolor(ms[0])









    Out[23]:





<matplotlib.collections.PolyCollection at 0x59e6e90>



In [14]:

    
tm = rc.triang_from_params(bp[0][1] + 100*diff)
pcolor(np.dot(tm, tm.T))









    Out[14]:





<matplotlib.collections.PolyCollection at 0x50f8490>



In [135]:

    
m0 = ms[0]
s = 100
r2 = np.zeros(s)
flex = np.zeros(s)
isoc1 = np.zeros(s)
isoc2 = np.zeros(s)

iso = np.ones(m0.shape[0])/np.sqrt(m0.shape[0])
diff = (bp[3][1] - bp[0][1])/100
p = bp[0][1]

for i in xrange(100):
    r2[i] = CalcR2(m0)
    flex[i] = flexibility(m0)
    isoc1[i] = np.abs(np.dot(eig(m0)[1][:,0], iso)) 
    isoc2[i] = np.abs(np.dot(eig(m0)[1][:,1], iso)) 

    p += diff
    new_b = rc.triang_from_params(p)
    m0 = np.dot(new_b, new_b.T)



In [156]:

    
scatter(r2[0], flex[0])









    Out[156]:





<matplotlib.collections.PathCollection at 0x78ede90>



In [37]:

    
CalcR2(ms[0])









    Out[37]:





0.071071941176470629



In [14]:

    
eig(ms[2])[1][:,0]









    Out[14]:





array([-0.14639117, -0.18927692, -0.17999167, -0.11383435, -0.13866206,
       -0.1512532 , -0.03624093, -0.20780105, -0.25858205, -0.1263813 ,
       -0.2860151 , -0.23674419, -0.24766999, -0.15076605, -0.26594475,
       -0.19303911, -0.21734474, -0.08107354, -0.07217581, -0.16727413,
       -0.08215141, -0.16877756, -0.15115711, -0.12454392, -0.13828721,
       -0.05198079, -0.22440013, -0.13631722, -0.21921859, -0.16375375,
       -0.0965526 , -0.1552577 , -0.1792374 , -0.05408101, -0.05769688])



In [15]:

    
eig(ms[3])[1][:,0]









    Out[15]:





array([ 0.00548838, -0.12234011, -0.20956893, -0.14645006, -0.23492048,
       -0.20782179, -0.10886239, -0.16620609, -0.25482124, -0.03853063,
       -0.21770677, -0.06025946, -0.25561981, -0.26166642, -0.27268206,
       -0.25697442, -0.24683595, -0.07062001, -0.11089805, -0.19085818,
       -0.09255424, -0.08386025, -0.05861632, -0.11048438, -0.12644895,
       -0.07465616, -0.20866868, -0.02547831, -0.16566946, -0.15009225,
       -0.16838086, -0.05432352, -0.27710192, -0.13242886, -0.0538359 ])



In [44]:

    
def calc_path(i, j, s=100, pcs=[0]):
    m0 = ms[i]
    r2 = np.zeros(s)
    flex = np.zeros(s)
    isoc = []
    
    iso = np.ones(m0.shape[0])/np.sqrt(m0.shape[0])
    diff = (bp[j][1] - bp[i][1])/100
    p = bp[i][1]
    
    for i in xrange(100):
        r2[i] = CalcR2(m0)
        flex[i] = flexibility(m0)
        for j in pcs:
            isoc.append(np.abs(np.dot(eig(m0)[1][:,j], iso)))
        
        p += diff
        new_b = rc.triang_from_params(p)
        m0 = np.dot(new_b, new_b.T)
    
    return r2, flex, isoc



In [71]:

    
rn, fn, ino = calc_path(2,3)



In [73]:

    
scatter(rn, fn, c=range(0,100))
annotate('3', xy=(rn[0], fn[0]),  xycoords='data',
                xytext=(0.05, 0.06), textcoords='axes fraction',
                arrowprops=dict(facecolor='blue', shrink=0.05),
                horizontalalignment='right', verticalalignment='top',
                )

annotate('4', xy=(rn[-1], fn[-1]),  xycoords='data',
                xytext=(0.05, 0.2), textcoords='axes fraction',
                arrowprops=dict(facecolor='red', shrink=0.05),
                horizontalalignment='right', verticalalignment='top',
                )
savefig('3-4-r2_vs_flex.png')



In [118]:

    
y = np.sqrt(np.diagonal(ms[]))
y = (y[:, np.newaxis] * y)
corrP = ms[5] / y



In [119]:

    
ms[4] = corrP



In [125]:

    
plt.ax









    Out[125]:





array([[ 1.        , -0.13356862,  0.2657388 , ...,  0.07197488,
        -0.05065331,  0.17577846],
       [-0.13356862,  1.        ,  0.45599007, ...,  0.08454455,
         0.32414673, -0.06027531],
       [ 0.2657388 ,  0.45599007,  1.        , ...,  0.07615338,
         0.1457147 ,  0.10191853],
       ..., 
       [ 0.07197488,  0.08454455,  0.07615338, ...,  1.        ,
         0.28592874,  0.31562355],
       [-0.05065331,  0.32414673,  0.1457147 , ...,  0.28592874,
         1.        ,  0.22368994],
       [ 0.17577846, -0.06027531,  0.10191853, ...,  0.31562355,
         0.22368994,  1.        ]])



In [109]:

    
t4 = rc.triang_decomp(ms[4])



In [111]:

    
np.max(np.dot(t4, t4.T) - ms[4])









    Out[111]:





7.7715611723760958e-16



In [154]:

    
np.savetxt('np5.csv', ms[4], delimiter=',')



In [47]:

    
#m7 = np.loadtxt(open("../matrices/7.csv","rb"),delimiter=";")
#m8 = np.loadtxt(open("../matrices/8.csv","rb"),delimiter=";")

m7 = rc.rand_corr(35, 10**-3)
m8 = rc.rand_corr(35, 10**-3)

ms =[m7, m8]



In [31]:

    
bp = []

for i in xrange(2):
    b_n = rc.triang_decomp(ms[i])
    p_n = rc.calc_params(b_n)
    bp.append((b_n, p_n))



In [38]:

    
import randomCorr as rc
reload(rc)









    Out[38]:





<module 'randomCorr' from '/home/walrus/corr-param/randomCorr.pyc'>



In [43]:

    
r,f,i = rc.calc_path(ms, bp, 0,1, pcs=[0])



In [44]:

    
plot(r)









    Out[44]:





[<matplotlib.lines.Line2D at 0x52e6990>]



In [41]:

    
plot(f)









    Out[41]:





[<matplotlib.lines.Line2D at 0x4ee0d10>]



In [46]:

    
plot(i)









    Out[46]:





[<matplotlib.lines.Line2D at 0x583bc50>]



In [14]:

    
scatter(r, f, c=range(0,100))
annotate('7', xy=(r[0], f[0]),  xycoords='data',
                xytext=(0.05, 0.06), textcoords='axes fraction',
                arrowprops=dict(facecolor='blue', shrink=0.05),
                horizontalalignment='right', verticalalignment='top',
                )

annotate('8', xy=(r[-1], f[-1]),  xycoords='data',
                xytext=(0.05, 0.2), textcoords='axes fraction',
                arrowprops=dict(facecolor='red', shrink=0.05),
                horizontalalignment='right', verticalalignment='top',
                )
savefig('7-8_r2_vs_flex.png')



In [32]:

    
def calc_path(mi, mj, s=100, pcs=[0]):
    m0 = mi
    r2 = np.zeros(s)
    flex = np.zeros(s)
    isoc = []
    
    iso = np.ones(m0.shape[0])/np.sqrt(m0.shape[0])
    diff = (mj - mi)/100
    
    for i in xrange(100):
        r2[i] = CalcR2(m0)
        flex[i] = flexibility(m0)
        for j in pcs:
            isoc.append(np.abs(np.dot(eig(m0)[1][:,j], iso)))
        
        m0 += diff
    
    return r2, flex, isoc



In [48]:

    
r,f,i=calc_path(ms[0], ms[1])



In [49]:

    
plot(r)









    Out[49]:





[<matplotlib.lines.Line2D at 0x5a76950>]



In [50]:

    
plot(f)









    Out[50]:





[<matplotlib.lines.Line2D at 0x5c8e650>]



In [51]:

    
plot(i)









    Out[51]:





[<matplotlib.lines.Line2D at 0x5cc3350>]



In [53]:

    
scatter(r,f)









    Out[53]:





<matplotlib.collections.PathCollection at 0x61d0c10>



In [54]:

    
m = rc.rand_corr(35, 10**-3)



In [31]:

    
import randomCorr as rc

r2 = []
flex = []
ms = []
for x in xrange(1000):
    m = rc.rand_corr(35, 10**-3)
    r2.append(rc.calc_r2(m))
    flex.append(rc.flexibility(m))
    ms.append(m)



In [61]:

    
lm = map(lambda m: m[3][7], ms)
hist(lm)









    Out[61]:





(array([ 48,  79, 101, 112, 134, 149, 140, 111,  86,  40]),
 array([-0.98154615, -0.78466887, -0.58779159, -0.39091432, -0.19403704,
        0.00284024,  0.19971751,  0.39659479,  0.59347207,  0.79034934,
        0.98722662]),
 <a list of 10 Patch objects>)



In [27]:

    
r_f = zip(r2,flex)
rf = filter(lambda r: r[0] > 0.2 and r[0] < 0.23, r_f)
rs = map(lambda i: i[0], rf)
fs = map(lambda i: i[1], rf)
scatter(rs,fs, c=range(0, len(rs)))
print corrcoef(rs,fs)









    



[[ 1.         -0.71620113]
 [-0.71620113  1.        ]]



In [29]:

    
scatter(r2, flex, c=range(1000))
print corrcoef(r2, flex)









    



[[ 1.         -0.94002555]
 [-0.94002555  1.        ]]



In [ ]:

    
np.mean(



In [ ]: