python求落后交叉相关系数

huyy

1. import numpy as np
   
2. from math import sqrt
   
3. 
   
4. import matplotlib.pyplot as plt
   
5. 
   
6. #求落后交叉相关系数
   
7. 
   
8. zry=[]
   
9. 
   
10. zrw=[]
    
11. 
    
12. for j in range(1,10):
    
13. 
    
14. def s(a,b):#协方差
    
15. 
    
16. avofa=0.0
    
17. 
    
18. avofb=0.0
    
19. 
    
20. sumofab=0.0
    
21. 
    
22. for i in range(len(a)):
    
23. 
    
24. avofa+=a[i]/20.
    
25. 
    
26. for i in range(len(b)):
    
27. 
    
28. avofb+=b[i]/20.
    
29. 
    
30. for i in range(len(a)-j):
    
31. 
    
32. temp=(a[i]-avofa)*(b[i+j]-avofb)
    
33. 
    
34. sumofab+=temp/(len(a)-j)
    
35. 
    
36. return sumofab
    
37. 
    
38. def variance(a):#平方的期望-期望的平方
    
39. 
    
40. s1=0;
    
41. 
    
42. s2=0;
    
43. 
    
44. for i in a:
    
45. 
    
46. s1+=i**2;
    
47. 
    
48. s2+=i;
    
49. 
    
50. return float(s1)/len(a)-(float(s2)/len(a))**2
    
51. 
    
52. def r(a,b):#落后交叉相关系数
    
53. 
    
54. return s(a,b)/sqrt((variance(a)*variance(b)))
    
55. 
    
56. 
    
57. 
    
58. f1=open('F:\\shixi3\\shixi312.txt')
    
59. 
    
60. data1=f1.readlines()
    
61. 
    
62. f1.close()
    
63. 
    
64. f2=open('F:\\shixi3\\shixi31.txt')
    
65. 
    
66. data2=f2.readlines()
    
67. 
    
68. f2.close()
    
69. 
    
70. f3=open('F:\\shixi3\\shixi32.txt')
    
71. 
    
72. data3=f3.readlines()
    
73. 
    
74. f3.close()
    
75. 
    
76. 
    
77. 
    
78. ave12=[]
    
79. 
    
80. ave1=[]
    
81. 
    
82. ave2=[]
    
83. 
    
84. for i in range(0,20):
    
85. 
    
86. temp1=data1[i]
    
87. 
    
88. temp2=data2[i]
    
89. 
    
90. temp3=data3[i]
    
91. 
    
92. ave12.append(float(temp1.split()[0]))
    
93. 
    
94. ave1.append(float(temp2.split()[0]))
    
95. 
    
96. ave2.append(float(temp3.split()[0]))
    
97. 
    
98. # print 'j=',j
    
99. 
    
100. # print r(yearave),',',r(winterave)
     
101. 
     
102. zry.append(float(r(ave12,ave1)))
     
103. 
     
104. zrw.append(float(r(ave12,ave2)))
     
105. 
     
106. #print zry,zrw
     
107. 
     
108. x=np.linspace(1,9,9)
     
109. 
     
110. y=np.linspace(-0.4,0.4,9)
     
111. 
     
112. cp=plt.plot(x,zry,marker='.',c='black',markersize=10,label='Month12 and month1')
     
113. 
     
114. cp=plt.plot(x,zrw,marker=6,c='gray',markersize=10,label='Month12 and month2')
     
115. 
     
116. 
     
117. 
     
118. for i in range(0,len(zry)):
     
119. 
     
120. plt.text(x[i],zry[i]+0.005,'%.02f' %zry[i],ha='center',va='bottom',color='black',fontsize=9)
     
121. 
     
122. for i in range(0,len(zrw)):
     
123. 
     
124. plt.text(x[i],zrw[i]+0.005,'%.02f' %zrw[i],ha='center',va='bottom',color='gray',fontsize=9)
     
125. 
     
126. plt.title('Backwardness cross correlation coefficient with different lag length',fontsize=10)
     
127. 
     
128. plt.xlabel('Lag Length',fontsize=10)
     
129. 
     
130. plt.ylabel('Backwardness cross correlation coefficient',fontsize=10)
     
131. 
     
132. 
     
133. 
     
134. plt.legend()
     
135. 
     
136. plt.show()
     

复制代码

又是那隻貓

不知道 python要对齐吗

huyy

又是那隻貓 发表于 2018-5-5 19:03
不知道 python要对齐吗

不要对齐，我好像复制的时候出现了一些问题。。