[GitHub]Machine Learning for the Web

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

1. Machine Learning For The Web
   
2. 
   
3. Each folder contains the codes discussed in each chapter of the book

复制代码

本帖隐藏的内容

Machine Learning for the Web.zip (99.55 MB)

2017-5-20 19:30:30 上传

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：Learning machine earning GitHub Learn

1. import numpy as np
   
2. from sklearn import mixture
   
3. from scipy.cluster.hierarchy import dendrogram, linkage
   
4. from scipy.cluster.hierarchy import fcluster
   
5. from sklearn.cluster import KMeans
   
6. from sklearn.cluster import MeanShift
   
7. from matplotlib import pyplot as plt
   
8. # generate two clusters: a with 100 points, b with 50:
   
9. np.random.seed(4711)  # for repeatability
   
10. c1 = np.random.multivariate_normal([10, 0], [[3, 1], [1, 4]], size=[100,])
    
11. l1 = np.zeros(100)
    
12. l2 = np.ones(100)
    
13. c2 = np.random.multivariate_normal([0, 10], [[3, 1], [1, 4]], size=[100,])
    
14. 
    
15. print c1.shape
    
16. #add noise:
    
17. np.random.seed(1)  # for repeatability
    
18. noise1x = np.random.normal(0,2,100)
    
19. noise1y = np.random.normal(0,8,100)
    
20. noise2 = np.random.normal(0,8,100)
    
21. c1[:,0] += noise1x
    
22. c1[:,1] += noise1y
    
23. c2[:,1] += noise2
    
24. 
    
25. #
    
26. fig = plt.figure(figsize=(20,15))
    
27. ax = fig.add_subplot(111)
    
28. ax.set_xlabel('x',fontsize=30)
    
29. ax.set_ylabel('y',fontsize=30)
    
30. fig.suptitle('classes',fontsize=30)
    
31. labels = np.concatenate((l1,l2),)
    
32. X = np.concatenate((c1, c2),)
    
33. pp1= ax.scatter(c1[:,0], c1[:,1],cmap='prism',s=50,color='r')
    
34. pp2= ax.scatter(c2[:,0], c2[:,1],cmap='prism',s=50,color='g')
    
35. ax.legend((pp1,pp2),('class 1', 'class2'),fontsize=35)
    
36. fig.savefig('classes.png')
    
37. 
    
38. 
    
39. #start figure
    
40. fig.clf()#reset plt
    
41. fig, ((axis1, axis2), (axis3, axis4)) = plt.subplots(2, 2, sharex='col', sharey='row')
    
42. 
    
43. #k-means
    
44. kmeans = KMeans(n_clusters=2)
    
45. kmeans.fit(X)
    
46. pred_kmeans = kmeans.labels_
    
47. #axis1 = fig.add_subplot(211)
    
48. print 'kmeans:',np.unique(kmeans.labels_)
    
49. print 'kmeans:',homogeneity_completeness_v_measure(labels,pred_kmeans)
    
50. plt.scatter(X[:,0], X[:,1], c=kmeans.labels_, cmap='prism')  # plot points with cluster dependent colors
    
51. axis1.scatter(X[:,0], X[:,1], c=kmeans.labels_, cmap='prism')
    
52. #axis1.set_xlabel('x',fontsize=40)
    
53. axis1.set_ylabel('y',fontsize=40)
    
54. axis1.set_title('k-means',fontsize=20)
    
55. #plt.show()
    
56. 
    
57. 
    
58. #mean-shift
    
59. ms = MeanShift(bandwidth=7)
    
60. ms.fit(X)
    
61. pred_ms = ms.labels_
    
62. axis2.scatter(X[:,0], X[:,1], c=pred_ms, cmap='prism')
    
63. axis2.set_title('mean-shift',fontsize=20)
    
64. 
    
65. print 'ms:',homogeneity_completeness_v_measure(labels,pred_ms)
    
66. print 'ms:',np.unique(ms.labels_)
    
67. 
    
68. #gaussian mixture
    
69. g = mixture.GMM(n_components=2)
    
70. g.fit(X)
    
71. print g.means_
    
72. pred_gmm = g.predict(X)
    
73. print 'gmm:',homogeneity_completeness_v_measure(labels,pred_gmm)
    
74. axis3.scatter(X[:,0], X[:,1], c=pred_gmm, cmap='prism')
    
75. axis3.set_xlabel('x',fontsize=40)
    
76. axis3.set_ylabel('y',fontsize=40)
    
77. axis3.set_title('gaussian mixture',fontsize=20)
    
78. 
    
79. 
    
80. 
    
81. #hierarchical
    
82. # generate the linkage matrix
    
83. Z = linkage(X, 'ward')
    
84. max_d = 20
    
85. pred_h = fcluster(Z, max_d, criterion='distance')
    
86. print 'clusters:',np.unique(pred_h)
    
87. k=2
    
88. fcluster(Z, k, criterion='maxclust')
    
89. print 'h:',homogeneity_completeness_v_measure(labels,pred_h)
    
90. axis4.scatter(X[:,0], X[:,1], c=pred_h, cmap='prism')
    
91. axis4.set_xlabel('x',fontsize=40)
    
92. #axis4.set_ylabel('y',fontsize=40)
    
93. axis4.set_title('hierarchical ward',fontsize=20)
    
94. 
    
95. fig.set_size_inches(18.5,10.5)
    
96. fig.savefig('comp_clustering.png', dpi=100)
    
97. 
    
98. fig.clf()#reset plt
    
99. fig = plt.figure(figsize=(20,15))
    
100. plt.title('Hierarchical Clustering Dendrogram',fontsize=30)
     
101. plt.xlabel('data point index (or cluster index)',fontsize=30)
     
102. plt.ylabel('distance (ward)',fontsize=30)
     
103. dendrogram(
     
104.     Z,
     
105.     truncate_mode='lastp',  # show only the last p merged clusters
     
106.     p=12,  # show only the last p merged clusters
     
107.     leaf_rotation=90.,
     
108.     leaf_font_size=12.,
     
109.     show_contracted=True,  # to get a distribution impression in truncated branches
     
110. )
     
111. fig.savefig('dendrogram.png')
     
112. #plt.show()
     
113. 
     
114. 
     
115. #measures:
     
116. from sklearn.metrics import homogeneity_completeness_v_measure
     
117. from sklearn.metrics import silhouette_score
     
118. 
     
119. res = homogeneity_completeness_v_measure(labels,pred_kmeans)
     
120. print 'kmeans measures, homogeneity:',res[0],' completeness:',res[1],' v-measure:',res[2],' silhouette score:',silhouette_score(X,pred_kmeans)
     
121. res = homogeneity_completeness_v_measure(labels,pred_ms)
     
122. print 'mean-shift measures, homogeneity:',res[0],' completeness:',res[1],' v-measure:',res[2],' silhouette score:',silhouette_score(X,pred_ms)
     
123. res = homogeneity_completeness_v_measure(labels,pred_gmm)
     
124. print 'gaussian mixture model measures, homogeneity:',res[0],' completeness:',res[1],' v-measure:',res[2],' silhouette score:',silhouette_score(X,pred_gmm)
     
125. res = homogeneity_completeness_v_measure(labels,pred_h)
     
126. print 'hierarchical (ward) measures, homogeneity:',res[0],' completeness:',res[1],' v-measure:',res[2],' silhouette score:',silhouette_score(X,pred_h)

复制代码

1. import numpy as np
   
2. from matplotlib import pyplot as plt
   
3. 
   
4. #line y = 2*x
   
5. x = np.arange(1,101,1).astype(float)
   
6. y = 5*np.arange(1,101,1).astype(float)
   
7. #add noise
   
8. noise = np.random.normal(0, 10, 100)
   
9. y += noise
   
10. 
    
11. fig = plt.figure(figsize=(10,10))
    
12. #plot
    
13. plt.plot(x,y,'ro')
    
14. plt.axis([0,102, -20,220])
    
15. plt.quiver(60, 100,10-0, 20-0, scale_units='xy', scale=1)
    
16. plt.arrow(60, 100,10-0, 20-0,head_width=2.5, head_length=2.5, fc='k', ec='k')
    
17. plt.text(70, 110, r'$v^1$', fontsize=20)
    
18. #plt.show()
    
19. 
    
20. #save
    
21. 
    
22. ax = fig.add_subplot(111)
    
23. ax.axis([0,102, -20,220])
    
24. ax.set_xlabel('x',fontsize=40)
    
25. ax.set_ylabel('y',fontsize=40)
    
26. fig.suptitle('2 dimensional dataset',fontsize=40)
    
27. fig.savefig('pca_data.png')
    
28. 
    
29. 
    
30. #calc PCA
    
31. mean_x = np.mean(x)
    
32. mean_y = np.mean(y)
    
33. mean_vector = np.array([[mean_x],[mean_y]])
    
34. u_x = (x- mean_x)/np.std(x)
    
35. u_y = (y-mean_y)/np.std(y)
    
36. sigma = np.cov([u_x,u_y])
    
37. print sigma
    
38. eig_vals, eig_vecs = np.linalg.eig(sigma)
    
39. 
    
40. eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i])
    
41.              for i in range(len(eig_vals))]
    
42.             
    
43. eig_pairs.sort()
    
44. eig_pairs.reverse()
    
45. print eig_pairs
    
46. v1 = eig_pairs[0][1]
    
47. #leading eigenvector:
    
48. x_v1 = v1[0]*np.std(x)+mean_x
    
49. y_v1 = v1[1]*np.std(y)+mean_y
    
50. print x_v1,'-',y_v1,'slope:',(y_v1)/(x_v1)
    
51. 
    
52. from sklearn.decomposition import PCA
    
53. #X = np.array([x,y])
    
54. X = np.array([u_x,u_y])
    
55. X = X.T
    
56. #print X
    
57. pca = PCA(n_components=1)
    
58. pca.fit(X)
    
59. V = pca.components_
    
60. print V,'-',V[0][1]/V[0][0]
    
61. #transform in reduced space
    
62. X_red_sklearn = pca.fit_transform(X)
    
63. print X_red_sklearn.shape
    
64. 
    
65. W = np.array(v1.reshape(2,1))
    
66. X_red = W.T.dot(X.T)
    
67. #check the reduced matrices are equal
    
68. assert X_red.T.all() == X_red_sklearn.all(), 'problem with the pca algorithm'
    
69. print X_red.T[0],'-',X_red_sklearn[0]

复制代码

1. import numpy as np
   
2. from copy import copy
   
3. 
   
4. class HMM:
   
5.     def __init__(self,pi,A,B):
   
6.         self.pi = pi
   
7.         self.A = A
   
8.         self.B = B
   
9.         
   
10.     def MostLikelyStateSequence(self,observations):
    
11.         
    
12.         #calc combinations:
    
13.         N = self.A.shape[0]
    
14.         T = len(observations)
    
15.         sequences = [str(i) for i in range(N)]
    
16.         probs = np.array([self.pi[i]*self.B[i,observations[0]] for i in range(N)])
    
17.         print probs
    
18.         for i in range(1,T):
    
19.             newsequences = []
    
20.             newprobs = np.array([])
    
21.             for s in range(len(sequences)):
    
22.                 for j in range(N):
    
23.                     newsequences.append(sequences[s]+str(j))
    
24.                     bef = int(sequences[s][-1])
    
25.                     tTpprob = probs[s]*self.A[bef,j]*self.B[j,observations[i]]
    
26.                     newprobs = np.append(newprobs,[tTpprob])
    
27.                     print sequences[s]+str(j),'-',tTpprob
    
28.             sequences = newsequences
    
29.             probs = newprobs
    
30.         return max((probs[i],sequences[i]) for i in range(len(sequences)))
    
31.             
    
32.     def ViterbiSequence(self,observations):
    
33.         deltas = [{}]
    
34.         seq = {}
    
35.         N = self.A.shape[0]
    
36.         states = [i for i in range(N)]
    
37.         T = len(observations)
    
38.         #initialization
    
39.         for s in states:
    
40.             deltas[0][s] = self.pi[s]*self.B[s,observations[0]]
    
41.             seq[s] = [s]
    
42.         #compute Viterbi
    
43.         for t in range(1,T):
    
44.             deltas.append({})
    
45.             newseq = {}
    
46.             for s in states:
    
47.                 (delta,state) = max((deltas[t-1][s0]*self.A[s0,s]*self.B[s,observations[t]],s0) for s0 in states)
    
48.                 deltas[t][s] = delta
    
49.                 newseq[s] = seq[state] + [s]
    
50.             seq = newseq
    
51.             
    
52.         (delta,state) = max((deltas[T-1][s],s) for s in states)
    
53.         return  delta,' sequence: ', seq[state]
    
54.         
    
55.     def maxProbSequence(self,observations):
    
56.         N = self.A.shape[0]
    
57.         states = [i for i in range(N)]
    
58.         T = len(observations)
    
59.         M = self.B.shape[1]
    
60.         # alpha_t(i) = P(O_1 O_2 ... O_t, q_t = S_i | hmm)
    
61.         # Initialize alpha
    
62.         alpha = np.zeros((N,T))
    
63.         c = np.zeros(T) #scale factors
    
64.         alpha[:,0] = pi.T * self.B[:,observations[0]]
    
65.         c[0] = 1.0/np.sum(alpha[:,0])
    
66.         alpha[:,0] = c[0] * alpha[:,0]
    
67.         # Update alpha for each observation step
    
68.         for t in range(1,T):
    
69.             alpha[:,t] = np.dot(alpha[:,t-1].T, self.A).T * self.B[:,observations[t]]
    
70.             c[t] = 1.0/np.sum(alpha[:,t])
    
71.             alpha[:,t] = c[t] * alpha[:,t]
    
72. 
    
73.         # beta_t(i) = P(O_t+1 O_t+2 ... O_T | q_t = S_i , hmm)
    
74.         # Initialize beta
    
75.         beta = np.zeros((N,T))
    
76.         beta[:,T-1] = 1
    
77.         beta[:,T-1] = c[T-1] * beta[:,T-1]
    
78.         # Update beta backwards froT end of sequence
    
79.         for t in range(len(observations)-1,0,-1):
    
80.             beta[:,t-1] = np.dot(self.A, (self.B[:,observations[t]] * beta[:,t]))
    
81.             beta[:,t-1] = c[t-1] * beta[:,t-1]
    
82.             
    
83.         norm = np.sum(alpha[:,T-1])
    
84.         seq = ''
    
85.         for t in range(T):
    
86.             g,state = max(((beta[i,t]*alpha[i,t])/norm,i) for i in states)
    
87.             seq +=str(state)
    
88.             
    
89.         return seq
    
90.         
    
91.     def simulate(self,time):
    
92. 
    
93.         def drawFromNormal(probs):
    
94.             return np.where(np.random.multinomial(1,probs) == 1)[0][0]
    
95. 
    
96.         observations = np.zeros(time)
    
97.         states = np.zeros(time)
    
98.         states[0] = drawFromNormal(self.pi)
    
99.         observations[0] = drawFromNormal(self.B[states[0],:])
    
100.         for t in range(1,time):
     
101.             states[t] = drawFromNormal(self.A[states[t-1],:])
     
102.             observations[t] = drawFromNormal(self.B[states[t],:])
     
103.         return observations,states
     
104. 
     
105. 
     
106.     def train(self,observations,criterion):
     
107. 
     
108.         N = self.A.shape[0]
     
109.         T = len(observations)
     
110.         M = self.B.shape[1]
     
111. 
     
112.         A = self.A
     
113.         B = self.B
     
114.         pi = copy(self.pi)
     
115.         
     
116.         convergence = False
     
117.         while not convergence:
     
118. 
     
119.             # alpha_t(i) = P(O_1 O_2 ... O_t, q_t = S_i | hmm)
     
120.             # Initialize alpha
     
121.             alpha = np.zeros((N,T))
     
122.             c = np.zeros(T) #scale factors
     
123.             alpha[:,0] = pi.T * self.B[:,observations[0]]
     
124.             c[0] = 1.0/np.sum(alpha[:,0])
     
125.             alpha[:,0] = c[0] * alpha[:,0]
     
126.             # Update alpha for each observation step
     
127.             for t in range(1,T):
     
128.                 alpha[:,t] = np.dot(alpha[:,t-1].T, self.A).T * self.B[:,observations[t]]
     
129.                 c[t] = 1.0/np.sum(alpha[:,t])
     
130.                 alpha[:,t] = c[t] * alpha[:,t]
     
131. 
     
132.             #P(O=O_0,O_1,...,O_T-1 | hmm)
     
133.             P_O = np.sum(alpha[:,T-1])
     
134.             # beta_t(i) = P(O_t+1 O_t+2 ... O_T | q_t = S_i , hmm)
     
135.             # Initialize beta
     
136.             beta = np.zeros((N,T))
     
137.             beta[:,T-1] = 1
     
138.             beta[:,T-1] = c[T-1] * beta[:,T-1]
     
139.             # Update beta backwards froT end of sequence
     
140.             for t in range(len(observations)-1,0,-1):
     
141.                 beta[:,t-1] = np.dot(self.A, (self.B[:,observations[t]] * beta[:,t]))
     
142.                 beta[:,t-1] = c[t-1] * beta[:,t-1]
     
143. 
     
144.             gi = np.zeros((N,N,T-1));
     
145. 
     
146.             for t in range(T-1):
     
147.                 for i in range(N):
     
148.                     
     
149.                     gamma_num = alpha[i,t] * self.A[i,:] * self.B[:,observations[t+1]].T * \
     
150.                             beta[:,t+1].T
     
151.                     gi[i,:,t] = gamma_num / P_O
     
152.   
     
153.             # gamma_t(i) = P(q_t = S_i | O, hmm)
     
154.             gamma = np.squeeze(np.sum(gi,axis=1))
     
155.             # Need final gamma element for new B
     
156.             prod =  (alpha[:,T-1] * beta[:,T-1]).reshape((-1,1))
     
157.             gamma_T = prod/P_O
     
158.             gamma = np.hstack((gamma,  gamma_T)) #append one Tore to gamma!!!
     
159. 
     
160.             newpi = gamma[:,0]
     
161.             newA = np.sum(gi,2) / np.sum(gamma[:,:-1],axis=1).reshape((-1,1))
     
162.             newB = copy(B)
     
163.             
     
164.             sumgamma = np.sum(gamma,axis=1)
     
165.             for ob_k in range(M):
     
166.                 list_k = observations == ob_k
     
167.                 newB[:,ob_k] = np.sum(gamma[:,list_k],axis=1) / sumgamma
     
168. 
     
169.             if np.max(abs(pi - newpi)) < criterion and \
     
170.                    np.max(abs(A - newA)) < criterion and \
     
171.                    np.max(abs(B - newB)) < criterion:
     
172.                 convergence = True;
     
173.   
     
174.             A[:],B[:],pi[:] = newA,newB,newpi
     
175. 
     
176.         self.A[:] = newA
     
177.         self.B[:] = newB
     
178.         self.pi[:] = newpi
     
179.         self.gamma = gamma
     
180.         
     
181. 
     
182. if __name__ == '__main__':
     
183.       
     
184.     pi = np.array([0.6, 0.4])
     
185.     A = np.array([[0.7, 0.3],
     
186.                            [0.6, 0.4]])
     
187.     B = np.array([[0.7, 0.1, 0.2],
     
188.                            [0.1, 0.6, 0.3]])
     
189.     hmmguess = HMM(pi,A,B)
     
190.     print 'Viterbi sequence:',hmmguess.ViterbiSequence(np.array([0,1,0,2]))
     
191.     print 'max prob sequence:',hmmguess.maxProbSequence(np.array([0,1,0,2]))   
     
192.     #obs,states = hmmguess.simulate(4)
     
193.    
     
194.     hmmguess.train(np.array([0,1,0,2]),0.000001)
     
195. 
     
196.     print 'Estimated initial probabilities:',hmmguess.pi
     
197. 
     
198.     print 'Estimated state transition probabililities:',hmmguess.A
     
199. 
     
200.     print 'Estimated observation probabililities:',hmmguess.B

复制代码

1. import pandas as pd
   
2. import numpy as np
   
3. from sklearn import cross_validation
   
4. from sklearn import svm
   
5. from sklearn.tree import DecisionTreeClassifier
   
6. from sklearn.ensemble import RandomForestClassifier
   
7. from sklearn.naive_bayes import MultinomialNB
   
8. from sklearn.linear_model import LogisticRegression
   
9. from sklearn.neighbors import KNeighborsClassifier
   
10. from sklearn.metrics import f1_score
    
11. from sklearn.metrics import precision_score
    
12. from sklearn.metrics import recall_score
    
13. 
    
14. #read data in
    
15. df = pd.read_csv('data_cars.csv',header=None)
    
16. for i in range(len(df.columns)):
    
17.     df[i] = df[i].astype('category')
    
18. df.head()
    
19. 
    
20. #map catgories to values
    
21. map0 = dict( zip( df[0].cat.categories, range( len(df[0].cat.categories ))))
    
22. #print map0
    
23. map1 = dict( zip( df[1].cat.categories, range( len(df[1].cat.categories ))))
    
24. map2 = dict( zip( df[2].cat.categories, range( len(df[2].cat.categories ))))
    
25. map3 = dict( zip( df[3].cat.categories, range( len(df[3].cat.categories ))))
    
26. map4 = dict( zip( df[4].cat.categories, range( len(df[4].cat.categories ))))
    
27. map5 = dict( zip( df[5].cat.categories, range( len(df[5].cat.categories ))))
    
28. map6 = dict( zip( df[6].cat.categories, range( len(df[6].cat.categories ))))
    
29. 
    
30. cat_cols = df.select_dtypes(['category']).columns
    
31. df[cat_cols] = df[cat_cols].apply(lambda x: x.cat.codes)
    
32. 
    
33. df = df.iloc[np.random.permutation(len(df))]
    
34. print df.head()
    
35. 
    
36. cv = 10
    
37. method = 'linear support vector machine'
    
38. clf = svm.SVC(kernel='linear',C=50)
    
39. y_pred = cross_validation.cross_val_predict(clf, X,Y, cv=cv)
    
40. CalcMeasures(method,y_pred,Y)
    
41. 
    
42. method = 'rbf support vector machine'
    
43. clf = svm.SVC(kernel='rbf',C=50)
    
44. y_pred = cross_validation.cross_val_predict(clf, X,Y, cv=cv)
    
45. CalcMeasures(method,y_pred,Y)
    
46. 
    
47. method = 'poly support vector machine'
    
48. clf = svm.SVC(kernel='poly',C=50)
    
49. y_pred = cross_validation.cross_val_predict(clf, X,Y, cv=cv)
    
50. CalcMeasures(method,y_pred,Y)
    
51. 
    
52. method = 'decision tree'
    
53. clf = DecisionTreeClassifier(random_state=0)
    
54. y_pred = cross_validation.cross_val_predict(clf, X,Y, cv=cv)
    
55. CalcMeasures(method,y_pred,Y)
    
56. 
    
57. method = 'random forest'
    
58. clf = RandomForestClassifier(n_estimators=50,random_state=0,max_features=None)
    
59. y_pred = cross_validation.cross_val_predict(clf, X,Y, cv=cv)
    
60. CalcMeasures(method,y_pred,Y)
    
61. 
    
62. method = 'naive bayes'
    
63. clf = MultinomialNB()
    
64. y_pred = cross_validation.cross_val_predict(clf, X,Y, cv=cv)
    
65. CalcMeasures(method,y_pred,Y)
    
66. 
    
67. method = 'logistic regression'
    
68. clf = LogisticRegression()
    
69. y_pred = cross_validation.cross_val_predict(clf, X,Y, cv=cv)
    
70. CalcMeasures(method,y_pred,Y)
    
71. 
    
72. method = 'k nearest neighbours'
    
73. clf = KNeighborsClassifier(weights='distance',n_neighbors=5)
    
74. y_pred = cross_validation.cross_val_predict(clf, X,Y, cv=cv)
    
75. CalcMeasures(method,y_pred,Y)

复制代码

1. In [1]:
   
2. import numpy as np
   
3. from matplotlib import pyplot as plt
   
4. %matplotlib inline
   
5. In [2]:
   
6. #line y = 2*x
   
7. x = np.arange(1,101,1).astype(float)
   
8. y = 2*np.arange(1,101,1).astype(float)
   
9. #add noise
   
10. noise = np.random.normal(0, 10, 100)
    
11. y += noise
    
12. 
    
13. fig = plt.figure(figsize=(10,10))
    
14. #plot
    
15. plt.plot(x,y,'ro')
    
16. plt.axis([0,102, -20,220])
    
17. plt.quiver(60, 100,10-0, 20-0, scale_units='xy', scale=1)
    
18. plt.arrow(60, 100,10-0, 20-0,head_width=2.5, head_length=2.5, fc='k', ec='k')
    
19. plt.text(70, 110, r'$v^1$', fontsize=20)
    
20. 
    
21. #save
    
22. ax = fig.add_subplot(111)
    
23. ax.axis([0,102, -20,220])
    
24. ax.set_xlabel('x',fontsize=40)
    
25. ax.set_ylabel('y',fontsize=40)
    
26. fig.suptitle('2 dimensional dataset',fontsize=40)
    
27. fig.savefig('pca_data.png')
    
28. 
    
29. In [3]:
    
30. #calc PCA
    
31. mean_x = np.mean(x)
    
32. mean_y = np.mean(y)
    
33. mean_vector = np.array([[mean_x],[mean_y]])
    
34. u_x = (x- mean_x)/np.std(x)
    
35. u_y = (y-mean_y)/np.std(y)
    
36. sigma = np.cov([u_x,u_y])
    
37. print sigma
    
38. eig_vals, eig_vecs = np.linalg.eig(sigma)
    
39. 
    
40. eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i])
    
41.              for i in range(len(eig_vals))]
    
42.             
    
43. eig_pairs.sort()
    
44. eig_pairs.reverse()
    
45. print eig_pairs
    
46. v1 = eig_pairs[0][1]
    
47. #leading eigenvector:
    
48. x_v1 = v1[0]*np.std(x)+mean_x
    
49. y_v1 = v1[1]*np.std(y)+mean_y
    
50. print 'slope:',(y_v1)/(x_v1)
    
51. 
    
52. from sklearn.decomposition import PCA
    
53. #X = np.array([x,y])
    
54. X = np.array([u_x,u_y])
    
55. X = X.T
    
56. #print X
    
57. pca = PCA(n_components=1)
    
58. pca.fit(X)
    
59. V = pca.components_
    
60. print V,'-',V[0][1]/V[0][0]
    
61. #transform in reduced space
    
62. X_red_sklearn = pca.fit_transform(X)
    
63. print X_red_sklearn.shape
    
64. 
    
65. W = np.array(v1.reshape(2,1))
    
66. X_red = W.T.dot(X.T)
    
67. #check the reduced matrices are equal
    
68. assert X_red.T.all() == X_red_sklearn.all(), 'problem with the pca algorithm'
    
69. print X_red.T[0],'-',X_red_sklearn[0]

复制代码

感谢楼主分享

谢谢分享