[Lecture Notes]Genevera Allen:Statistical Learning using R and Matlab

1. %%%%%%%%%%%%%%
   
2. %lecture 11 code
   
3. %matlab script - svm exmaple
   
4. %primal form problems
   
5. 
   
6. cd ~/cvx
   
7. cvx_setup
   
8. cd ~/classes/stat640_2015/
   
9. 
   
10. 
    
11. %%%%%%%%%%%
    
12. %separable case - optimal separating hyperplanes
    
13. 
    
14. clf
    
15. var = 1.25;
    
16. 
    
17. n = 100; mu1 = [3; 8]; mu2 = [8; 3];
    
18. x1 = repmat(mu1',n/2,1) + randn(n/2,2)*sqrt(var);
    
19. x2 = repmat(mu2',n/2,1) + randn(n/2,2)*sqrt(var);
    
20. x = [x1; x2];
    
21. Y = [repmat(1,50,1); repmat(-1,50,1)];
    
22. 
    
23. hold on
    
24. scatter(x1(:,1),x1(:,2),'or','filled')
    
25. scatter(x2(:,1),x2(:,2),'ob','filled')
    
26. axis([min(min(x(:,1))) max(max(x(:,1))) min(min(x(:,2))) max(max(x(:,2)))]);
    
27. hold off
    
28. 
    
29. cvx_begin
    
30. cvx_precision('high')
    
31. variable beta0;
    
32. variable beta(2);
    
33. minimize( .5*square_pos(norm(beta)));
    
34. subject to
    
35. Y.*(x*beta + beta0) >= 1;
    
36. cvx_end
    
37. 
    
38. 
    
39. hold on
    
40. scatter(x1(:,1),x1(:,2),'or','filled')
    
41. scatter(x2(:,1),x2(:,2),'ob','filled')
    
42. axis([min(min(x(:,1))) max(max(x(:,1))) min(min(x(:,2))) max(max(x(:,2)))]);
    
43. xs = linspace(min(min(x)),max(max(x)),1000);
    
44. plot(xs,(-beta0 - xs*beta(1))/beta(2),'k')
    
45. M = 1/norm(beta);
    
46. plot(xs,(-beta0 - xs*beta(1)+1)/beta(2) ,'--k')
    
47. plot(xs,(-beta0 - xs*beta(1)-1)/beta(2),'--k')
    
48. hold off
    
49. 
    
50. 
    
51. %%%%%%%%%%%%%%
    
52. %non-seperable case - linear SVMs
    
53. 
    
54. clf
    
55. var = 2.8;
    
56. 
    
57. n = 100; mu1 = [3; 8]; mu2 = [8; 3];
    
58. x1 = repmat(mu1',n/2,1) + randn(n/2,2)*sqrt(var);
    
59. x2 = repmat(mu2',n/2,1) + randn(n/2,2)*sqrt(var);
    
60. x = [x1; x2];
    
61. Y = [repmat(1,50,1); repmat(-1,50,1)];
    
62. 
    
63. hold on
    
64. scatter(x1(:,1),x1(:,2),'or','filled')
    
65. scatter(x2(:,1),x2(:,2),'ob','filled')
    
66. axis([min(min(x(:,1))) max(max(x(:,1))) min(min(x(:,2))) max(max(x(:,2)))]);
    
67. hold off
    
68. 
    
69. 
    
70. gam = 5;
    
71. 
    
72. cvx_begin
    
73. cvx_precision('high')
    
74. variable beta0;
    
75. variable beta(2);
    
76. variable epsilon(n);
    
77. minimize( .5*square_pos(norm(beta)) + gam*sum(epsilon));
    
78. subject to
    
79. for i=1:n
    
80.     Y(i)*(x(i,:)*beta + beta0) >= 1 - epsilon(i);
    
81. end
    
82. epsilon >= 0;
    
83. cvx_end
    
84. 
    
85. clf
    
86. hold on
    
87. scatter(x1(:,1),x1(:,2),'or','filled')
    
88. scatter(x2(:,1),x2(:,2),'ob','filled')
    
89. axis([min(min(x(:,1))) max(max(x(:,1))) min(min(x(:,2))) max(max(x(:,2)))]);
    
90. xs = linspace(min(min(x)),max(max(x)),1000);
    
91. plot(xs,(-beta0 - xs*beta(1))/beta(2),'k')
    
92. M = 1/norm(beta);
    
93. plot(xs,(-beta0 - xs*beta(1)+1)/beta(2) ,'--k')
    
94. plot(xs,(-beta0 - xs*beta(1)-1)/beta(2),'--k')
    
95. ind = epsilon>1e-6;
    
96. scatter(x(ind,1),x(ind,2),'+k','SizeData',150)
    
97. ind = abs(x*beta + beta0 + 1)<1e-4;
    
98. scatter(x(ind,1),x(ind,2),'Xk','SizeData',150)
    
99. ind = abs(x*beta + beta0 - 1)<1e-4;
    
100. scatter(x(ind,1),x(ind,2),'Xk','SizeData',150)
     
101. hold off
     
102. 
     
103. 
     
104. 
     
105. 
     
106. 
     
107. 
     
108. 
     
109. 
     
110. 
     
111. 
     
112. 
     
113. 
     
114. 
     
115. 
     
116. 
     
117. 
     
118. 
     
119. 
     
120. 
     
121. 
     
122. %%%%%%%%%%%%%%
     
123. %lecture 11 code
     
124. %matlab script - svm exmaple
     
125. %primal form problems
     
126. 
     
127. cd ~/cvx
     
128. cvx_setup
     
129. cd ~/classes/stat640_2014/
     
130. 
     
131. 
     
132. %%%%%%%%%%%%%%
     
133. %non-seperable case
     
134. 
     
135. clf
     
136. var = 2.8;
     
137. 
     
138. n = 100; mu1 = [3; 8]; mu2 = [8; 3];
     
139. x1 = repmat(mu1',n/2,1) + randn(n/2,2)*sqrt(var);
     
140. x2 = repmat(mu2',n/2,1) + randn(n/2,2)*sqrt(var);
     
141. x = [x1; x2];
     
142. Y = [repmat(1,50,1); repmat(-1,50,1)];
     
143. 
     
144. hold on
     
145. scatter(x1(:,1),x1(:,2),'or','filled')
     
146. scatter(x2(:,1),x2(:,2),'ob','filled')
     
147. axis([min(min(x(:,1))) max(max(x(:,1))) min(min(x(:,2))) max(max(x(:,2)))]);
     
148. hold off
     
149. 
     
150. 
     
151. gam = .1;
     
152. 
     
153. cvx_begin
     
154. cvx_precision('high')
     
155. variable beta0;
     
156. variable beta(2);
     
157. variable epsilon(n);
     
158. minimize( .5*square_pos(norm(beta)) + gam*sum(epsilon));
     
159. subject to
     
160. for i=1:n
     
161.     Y(i)*(x(i,:)*beta + beta0) >= 1 - epsilon(i);
     
162. end
     
163. epsilon >= 0;
     
164. cvx_end
     
165. 
     
166. clf
     
167. hold on
     
168. scatter(x1(:,1),x1(:,2),'or','filled')
     
169. scatter(x2(:,1),x2(:,2),'ob','filled')
     
170. axis([min(min(x(:,1))) max(max(x(:,1))) min(min(x(:,2))) max(max(x(:,2)))]);
     
171. xs = linspace(min(min(x)),max(max(x)),1000);
     
172. plot(xs,(-beta0 - xs*beta(1))/beta(2),'k')
     
173. M = 1/norm(beta);
     
174. plot(xs,(-beta0 - xs*beta(1)+1)/beta(2) ,'--k')
     
175. plot(xs,(-beta0 - xs*beta(1)-1)/beta(2),'--k')
     
176. ind = epsilon>1e-6;
     
177. scatter(x(ind,1),x(ind,2),'+k','SizeData',150)
     
178. ind = abs(x*beta + beta0 + 1)<1e-4;
     
179. scatter(x(ind,1),x(ind,2),'Xk','SizeData',150)
     
180. ind = abs(x*beta + beta0 - 1)<1e-4;
     
181. scatter(x(ind,1),x(ind,2),'Xk','SizeData',150)
     
182. hold off
     
183. 
     
184. %%%%%%%%%%%%%%%%%%%%%
     
185. %zip code example
     
186. %pair-wise
     
187. 
     
188. zip1 = 3; zip2 = 8;
     
189. 
     
190. trdat = dlmread('zip.train');
     
191. ind = trdat(:,1)==zip1 | trdat(:,1)==zip2;
     
192. X = trdat(ind,2:end);
     
193. Y = trdat(ind,1);
     
194. [n,p] = size(X);
     
195. tsdat = dlmread('zip.test');
     
196. ind = tsdat(:,1)==zip1 | tsdat(:,1)==zip2;
     
197. Xs = tsdat(ind,2:end);
     
198. Ys = tsdat(ind,1);
     
199. 
     
200. C = .01;
     
201. fitsvm = svmtrain(X,Y,'boxconstraint',C);
     
202. class = svmclassify(fitsvm,X);
     
203. trerr = sum(class~=Y)/length(Y);
     
204. class = svmclassify(fitsvm,Xs);
     
205. tserr = sum(class~=Ys)/length(Ys);
     
206. [trerr tserr]
     
207. fitsvm
     
208. 
     
209. colormap(gray)
     
210. for i=1:9
     
211.     subplot(3,3,i)
     
212.     imagesc(flipud(rot90(reshape(fitsvm.SupportVectors(i,:),16,16))))
     
213.     axis off
     
214. end
     
215. 
     
216. 
     
217. [classNB,errNB] = classify(Xs,X,Y,'diaglinear');
     
218. [classLDA,errLDA] = classify(Xs,X,Y,'linear');
     
219. [errNB errLDA tserr]

复制代码

1. %%%%%%%%%%%%%%%%%%%%%%
   
2. %lecture 12 code
   
3. 
   
4. 
   
5. %zip code data - 3's and 8's
   
6. 
   
7. zip1 = 3; zip2 = 8;
   
8. 
   
9. trdat = dlmread('zip.train');
   
10. ind = trdat(:,1)==zip1 | trdat(:,1)==zip2;
    
11. X = trdat(ind,2:end);
    
12. Y = trdat(ind,1);
    
13. [n,p] = size(X);
    
14. tsdat = dlmread('zip.test');
    
15. ind = tsdat(:,1)==zip1 | tsdat(:,1)==zip2;
    
16. Xs = tsdat(ind,2:end);
    
17. Ys = tsdat(ind,1);
    
18. 
    
19. %C-SVM - larger C gives a smaller margin (inversely propotional to gamma)
    
20. 
    
21. C = 1;
    
22. fitsvm = svmtrain(X,Y,'boxconstraint',C);
    
23. class = svmclassify(fitsvm,X);
    
24. trerr = sum(class~=Y)/length(Y);
    
25. class = svmclassify(fitsvm,Xs);
    
26. tserr = sum(class~=Ys)/length(Ys);
    
27. [trerr tserr]
    
28. fitsvm
    
29. 
    
30. 
    
31. %plot some of the support vectors
    
32. colormap(gray)
    
33. for i=1:9
    
34.     subplot(3,3,i)
    
35.     imagesc(-flipud(rot90(reshape(fitsvm.SupportVectors(i,:),16,16))))
    
36.     axis off
    
37. end
    
38. 
    
39. 
    
40. 
    
41. %compare test error to other methods
    
42. 
    
43. [classNB,errNB] = classify(Xs,X,Y,'diaglinear');
    
44. tserrNB = sum(classNB~=Ys)/length(Ys);
    
45. [classLDA,errLDA] = classify(Xs,X,Y,'linear');
    
46. tserrLDA = sum(classLDA~=Ys)/length(Ys);
    
47. [tserrNB tserrLDA tserr]
    
48. 
    
49. 
    
50. 
    
51. 
    
52. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
53. %non-linear 2-class SVM example
    
54. 
    
55. var = 1.25;
    
56. n = 100; mu1 = [0; 0]; cov1 = [1 .8; 1 .8];
    
57. x1 = repmat(mu1',n/2,1) + randn(n/2,2)*chol(cov1)*sqrt(var);
    
58. xs2 = repmat(mu1',n*100,1) + randn(n*100,2)*sqrt(var*3);
    
59. ind = sum((xs2.^2)')' > 5;
    
60. x2 = xs2(ind,:);
    
61. x2 = x2(1:(n/2),:);
    
62. x = [x1; x2];
    
63. Y = [repmat(1,50,1); repmat(-1,50,1)];
    
64. 
    
65. clf
    
66. hold on
    
67. scatter(x1(:,1),x1(:,2),'or','filled')
    
68. scatter(x2(:,1),x2(:,2),'ob','filled')
    
69. axis([min(min(x(:,1))) max(max(x(:,1))) min(min(x(:,2))) max(max(x(:,2)))]);
    
70. hold off
    
71. 
    
72. 
    
73. %linear kernel
    
74. C = 1;
    
75. svms = svmtrain(x,Y,'kernel_function','linear','showplot','true','method','QP','boxconstraint',C);
    
76. 
    
77. %radial kernel
    
78. sig = .5;
    
79. C = .1;
    
80. svms = svmtrain(x,Y,'kernel_function','rbf','showplot','true','method','QP','rbf_sigma',sig,'boxconstraint',C);
    
81. 
    
82. %polynomial kernel
    
83. d = 2;
    
84. C = 1;
    
85. svms = svmtrain(x,Y,'kernel_function','polynomial','showplot','true','method','QP','polyorder',d,'boxconstraint',C);
    
86. 
    
87. 
    
88. 
    
89. 
    
90. %%%%%%%%%%%%%%%%%%%%%%%
    
91. %non-linear regression example
    
92. 
    
93. n = 200; p = 1;
    
94. x = rand(n,1)*10;
    
95. y = sin(x)*4 + randn(n,1);
    
96. 
    
97. scatter(x,y,'filled')
    
98. 
    
99. %linear
    
100. betals = inv(x'*x)*x'*y;
     
101. hold on
     
102. scatter(x,y,'filled')
     
103. xs = linspace(0,10,500);
     
104. plot(xs,xs*betals,'k')
     
105. hold off
     
106. 
     
107. %polynomial
     
108. d = 2; lam = 5;
     
109. k = kern(x,[],ones(p,1),'poly1',d);
     
110. alpha = inv(k + lam*eye(n))*y;
     
111. clf
     
112. hold on
     
113. scatter(x,y,'filled')
     
114. xs = linspace(0,10,500);
     
115. kt = kern(x,xs',ones(p,1),'poly1',d);
     
116. plot(xs,kt*alpha,'b')
     
117. hold off
     
118. 
     
119. %gaussian
     
120. d = 1; lam = .01;
     
121. k = kern(x,[],ones(p,1),'gaussian',d);
     
122. alpha = inv(k + lam*eye(n))*y;
     
123. clf
     
124. hold on
     
125. scatter(x,y,'filled')
     
126. xs = linspace(0,10,500);
     
127. kt = kern(x,xs',ones(p,1),'gaussian',d);
     
128. plot(xs,kt*alpha,'b')
     
129. hold off

复制代码

1. ##################
   
2. #lecture 14 code
   
3. 
   
4. require(glmnet)
   
5. require(Hmisc)
   
6. 
   
7. #model complexity & prediction error
   
8. 
   
9. n = 100
   
10. p = 50
    
11. Btrue = matrix(0,p,1)
    
12. Btrue[1:10] = rnorm(10)*1;
    
13. 
    
14. Xtr = scale(matrix(rnorm(n*p),n,p))
    
15. Ytr = Xtr%*%Btrue + matrix(rnorm(n),n,1)
    
16. Xts = scale(matrix(rnorm(n*p),n,p))
    
17. Yts = Xts%*%Btrue + matrix(rnorm(n),n,1)
    
18. 
    
19. fit = glmnet(x=Xtr,y=Ytr,family="gaussian",standardize=FALSE,nlambda=50,lambda.min.ratio=.001)
    
20. Yhtr = predict(fit,newx=Xtr)
    
21. MSEtr = apply((Yhtr-Ytr%*%matrix(1,1,length(fit$lambda)))^2,2,mean)
    
22. Yhts = predict(fit,newx=Xts)
    
23. MSEts = apply((Yhts-Yts%*%matrix(1,1,length(fit$lambda)))^2,2,mean)
    
24. 
    
25. plot(1:length(fit$lambda),MSEtr,type="l",col=4,xlab="Sparsity",ylab="Error")
    
26. lines(1:length(fit$lambda),MSEtr,col=4)
    
27. lines(1:length(fit$lambda),MSEts,col=2)
    
28. legend(30,5,legend=c("Training Error","Test Error"),col=c(4,2),lty=c(1,1),cex=.75)
    
29. 
    
30. t(t(fit$beta[,15]))
    
31. 
    
32. 
    
33. #Cross-Validation (by hand)
    
34. #Model Selection
    
35. 
    
36. fold = 5
    
37. sam = sample(1:n,n)
    
38. CVerrs = NULL
    
39. for(i in 1:fold)
    
40. {
    
41.   ind = sam[((i-1)*n/fold + 1):(i*n/fold)]
    
42.   Xin = Xtr[-ind,]; Yin = Ytr[-ind]
    
43.   Xout = Xtr[ind,]; Yout = Ytr[ind]
    
44.   fit = glmnet(x=Xin,y=Yin,family="gaussian",standardize=FALSE,nlambda=50,lambda.min.ratio=.001)
    
45.   Yh = predict(fit,newx=Xout)
    
46.   CVerrs = cbind(CVerrs,apply((Yh-Yout%*%matrix(1,1,length(fit$lambda)))^2,2,mean))
    
47. }
    
48. CVerr = apply(CVerrs,1,mean)
    
49. 
    
50. #minimum CV error rule
    
51. lines(1:length(fit$lambda),CVerr,col=1)
    
52. optlam = fit$lambda[which.min(CVerr)]
    
53. 
    
54. #one SE rule
    
55. SE = sqrt(apply(CVerrs,1,var)/fold)
    
56. errbar(1:length(fit$lambda),CVerr,CVerr+SE,CVerr-SE,add=TRUE)
    
57. minSE = CVerr[which.min(CVerr)]+SE[which.min(CVerr)]
    
58. minSEx = which(CVerr<minSE)[1]
    
59. optlam = fit$lambda[minSEx]
    
60. 
    
61. 
    
62. #refit to training data at optimal lambda
    
63. fit = glmnet(x=Xtr,y=Ytr,family="gaussian",standardize=FALSE,lambda=optlam)
    
64. Yhtr = predict(fit,newx=Xtr)
    
65. TRerr = mean( (Yhtr - Ytr)^2 )
    
66. Yhts = predict(fit,newx=Xts)
    
67. TSerr = mean( (Yhts - Yts)^2 )
    
68. 
    
69. optlam
    
70. sum(fit$beta!=0)
    
71. #why doesn't CV find the "true" model?
    
72. 
    
73. TRerr
    
74. TSerr

复制代码

1. #############
   
2. #lecture 16 - PCA
   
3. 
   
4. ncidat = read.table("ncidata.txt")
   
5. labs = read.table("ncilabels.txt")
   
6. colnames(ncidat) = t(labs)
   
7. 
   
8. dim(ncidat)
   
9. unique(colnames(ncidat))
   
10. 
    
11. #PCA - take SVD to get solution
    
12. #center genes, but don't scale
    
13. X = t(scale(t(ncidat),center=TRUE,scale=FALSE))
    
14. sv = svd(t(X));
    
15. U = sv$u
    
16. V = sv$v
    
17. D = sv$d
    
18. Z = t(X)%*%V;
    
19. 
    
20. #PC scatterplots
    
21. cols = as.numeric(as.factor(colnames(ncidat)))
    
22. K = 3
    
23. pclabs = c("PC1","PC2","PC3","PC4")
    
24. par(mfrow=c(1,K))
    
25. for(i in 1:K){
    
26.   j = i+1
    
27.   plot(U[,i],U[,j],type="n",xlab=pclabs[i],ylab=pclabs[j])
    
28.   text(U[,i],U[,j],colnames(X),col=cols)
    
29. }
    
30. 
    
31. #PC loadings - visualize data by limiting to top genes in magnitude in the PC loadings
    
32. aa = grep("grey",colors())
    
33. bb = grep("green",colors())
    
34. cc = grep("red",colors())
    
35. gcol2 = colors()[c(aa[1:30],bb[1:20],rep(cc,2))]
    
36. 
    
37. j = 2
    
38. ord = order(abs(V[,j]),decreasing=TRUE)
    
39. x = as.matrix(X[ord[1:250],])
    
40. heatmap(x,col=gcol2)
    
41. 
    
42. #Variance Explained
    
43. varex = 0
    
44. cumvar = 0
    
45. denom = sum(D^2)
    
46. for(i in 1:64){
    
47.   varex[i] = D[i]^2/denom
    
48.   cumvar[i] = sum(D[1:i]^2)/denom
    
49. }
    
50. 
    
51. #screeplot
    
52. par(mfrow=c(1,2))
    
53. plot(1:64,varex,type="l",lwd=2,xlab="PC",ylab="% Variance Explained")
    
54. plot(1:64,cumvar,type="l",lwd=2,xlab="PC",ylab="Cummulative Variance Explained")
    
55. 
    
56. #######
    
57. #Sparse PCA
    
58. require("PMA")
    
59. 
    
60. spc = SPC(t(X),sumabsv=10,K=4)
    
61. 
    
62. #how many genes selected?
    
63. apply(spc$v!=0,2,sum)
    
64. 
    
65. #PC scatterplots
    
66. cols = as.numeric(as.factor(colnames(ncidat)))
    
67. K = 3
    
68. pclabs = c("SPC1","SPC2","SPC3","SPC4")
    
69. par(mfrow=c(1,K))
    
70. for(i in 1:K){
    
71.   j = i+1
    
72.   plot(spc$u[,i],spc$u[,j],type="n",xlab=pclabs[i],ylab=pclabs[j])
    
73.   text(spc$u[,i],spc$u[,j],colnames(X),col=cols)
    
74. }
    
75. 
    
76. #SPC loadings - visualize data by limiting to gene selected by the sparse PC loadings
    
77. aa = grep("grey",colors())
    
78. bb = grep("green",colors())
    
79. cc = grep("red",colors())
    
80. gcol2 = colors()[c(aa[1:30],bb[1:20],rep(cc,2))]
    
81. 
    
82. j = 1
    
83. ind = which(spc$v[,j]!=0)
    
84. x = as.matrix(X[ind,])
    
85. heatmap(x,col=gcol2)
    
86. 
    
87. #variance explained
    
88. spc$prop.var.explained

复制代码

1. #######################
   
2. #lecutre 17 code
   
3. #comparing dimension reduction methods
   
4. 
   
5. 
   
6. load("UnsupL.Rdata")
   
7. imagedigit = function(vec){
   
8.   mat = matrix(-vec[256:1],16,16)
   
9.   mat = apply(mat,2,rev)
   
10.   image(mat,col=gray((0:64)/64),xaxt="n",yaxt="n",ann=FALSE)
    
11. }
    
12. 
    
13. 
    
14. #pull out 3's
    
15. dat3 = digits[which(rownames(digits)==3),]
    
16. 
    
17. #visulaize
    
18. par(mfrow=c(3,4))
    
19. for(i in 1:12){
    
20.   imagedigit(dat3[i,])
    
21. }
    
22. 
    
23. #PCA - take SVD to get solution
    
24. #don't center and scale to retain interpretation as images
    
25. svd3 = svd(dat3)
    
26. U = svd3$u
    
27. V = svd3$v #PC loadings
    
28. D = svd3$d
    
29. Z = dat3%*%V #PCs
    
30. 
    
31. #PC scatterplot
    
32. par(mfrow=c(1,1))
    
33. plot(Z[,2],Z[,3],pch=16)
    
34. 
    
35. #PC loadings
    
36. par(mfrow=c(1,4))
    
37. for(i in 1:4){
    
38.   imagedigit(V[,i])
    
39. }
    
40. 
    
41. #Variance Explained
    
42. varex = 0
    
43. cumvar = 0
    
44. denom = sum(D^2)
    
45. for(i in 1:256){
    
46.   varex[i] = D[i]^2/denom
    
47.   cumvar[i] = sum(D[1:i]^2)/denom
    
48. }
    
49. 
    
50. #screeplot
    
51. par(mfrow=c(1,2))
    
52. plot(1:256,varex,type="l",lwd=2,xlab="PC",ylab="% Variance Explained")
    
53. plot(1:256,cumvar,type="l",lwd=2,xlab="PC",ylab="Cummulative Variance Explained")
    
54. 
    
55. 
    
56. cumvar[25] #first 25 PCs explain over 90% of variance
    
57. pdat3 = dat3%*%V[,1:25] #projected data - a tenth of the original size
    
58. 
    
59. 
    
60. #######
    
61. #now NMF
    
62. 
    
63. require("NMF")
    
64. 
    
65. K = 15
    
66. nmffit = nmf(dat3+1,rank=K)
    
67. W = basis(nmffit)
    
68. H = coef(nmffit)
    
69. 
    
70. #plot archetypes - try changing K
    
71. par(mfrow=c(3,5))
    
72. for(i in 1:K){
    
73.   imagedigit(H[i,])
    
74. }
    
75. 
    
76. 
    
77. ###########
    
78. #now ICA
    
79. 
    
80. require("fastICA")
    
81. 
    
82. K = 15
    
83. icafit = fastICA(t(dat3),n.comp=K)
    
84. 
    
85. #plot independent source signals - try changing K
    
86. par(mfrow=c(3,5))
    
87. for(i in 1:K){
    
88.   imagedigit(icafit$S[,i])
    
89. }

复制代码

1. ##################
   
2. #lecture 18 - K-means clustering
   
3. 
   
4. 
   
5. #simulate data
   
6. n = 300
   
7. mu1 = c(3,3); mu2 = c(7,4); mu3 = c(6,5.5);
   
8. Sig = matrix(c(1,.5,.5,1),2,2)
   
9. x1 = t(matrix(mu1,2,n/3)) + matrix(rnorm(n),n/3,2)
   
10. xx = matrix(rnorm(n*2/3),n/3,2)
    
11. x2 = t(matrix(mu2,2,n/3)) + xx%*%chol(Sig)
    
12. xx = matrix(rnorm(n*2/3),n/3,2)
    
13. x3 = t(matrix(mu3,2,n/3)) + xx%*%chol(Sig)
    
14. x = rbind(x1,x2,x3)
    
15. Y = c(rep(1,n/3),rep(2,n/3),rep(3,n/3))
    
16. y = factor(Y)
    
17. 
    
18. plot(x[,1],x[,2],col=as.numeric(y),pch=16)
    
19. 
    
20. 
    
21. #try changing k - which clustering looks best?
    
22. k = 3
    
23. km = kmeans(x,centers=k)
    
24. plot(x[,1],x[,2],col=km$cluster,pch=16)
    
25. cens = km$centers
    
26. points(cens[,1],cens[,2],col=1:k,pch=16,cex=3)
    
27. 
    
28. 
    
29. #code to understand K-means algorithm
    
30. require("animation")
    
31. mv.kmeans = function(x,k,cens=NULL){
    
32.   n = nrow(x)
    
33.   if(is.null(cens)){
    
34.       cens = x[sample(1:n,k),]
    
35.     }
    
36.   plot(x[,1],x[,2],pch=16)
    
37.   points(cens[,1],cens[,2],col=1:k,pch=16,cex=3)
    
38.   thr = 1e-6; ind = 1; iter = 1;
    
39.   while( ind>thr)
    
40.     {
    
41.       oldcen = cens
    
42.       km = kmeans(x,centers=cens,iter.max=1,nstart=1,algorithm="MacQueen")
    
43.       plot(x[,1],x[,2],col=km$cluster,pch=16)
    
44.       points(cens[,1],cens[,2],col=1:k,pch=16,cex=3)
    
45.       cens = km$centers
    
46.       #print(cens)
    
47.       plot(x[,1],x[,2],col=km$cluster,pch=16)
    
48.       points(cens[,1],cens[,2],col=1:k,pch=16,cex=3)
    
49.       ind = sum(diag((oldcen-cens)%*%t(oldcen-cens)))
    
50.       #print(ind)
    
51.     }
    
52. }
    
53. 
    
54. #watch K-means algorithm movie
    
55. #start from random starting points
    
56. saveHTML(mv.kmeans(x,3,cens=NULL),img.name="km2015")
    
57. 
    
58. 
    
59. 
    
60. 
    
61. 
    
62. ############
    
63. #K-means for high-dimensional data
    
64. 
    
65. require("ISLR")
    
66. 
    
67. ncidat = t(NCI60$data)
    
68. colnames(ncidat) = NCI60$labs
    
69. 
    
70. dim(ncidat)
    
71. unique(colnames(ncidat))
    
72. 
    
73. 
    
74. K = 9
    
75. km = kmeans(t(ncidat),centers=K)
    
76. 
    
77. #how do we visualize K-means results?
    
78. 
    
79. #PCA - take SVD to get solution
    
80. #center genes, but don't scale
    
81. X = t(scale(t(ncidat),center=TRUE,scale=FALSE))
    
82. sv = svd(t(X));
    
83. U = sv$u
    
84. V = sv$v
    
85. D = sv$d
    
86. Z = t(X)%*%V;
    
87. 
    
88. plot(Z[,1],Z[,2],col=km$cluster,type="n")
    
89. text(Z[,1],Z[,2],colnames(ncidat),cex=.75,col=km$cluster)
    
90. cens = km$centers
    
91. points(cens%*%V[,1],cens%*%V[,2],col=1:K,pch=16,cex=3)
    
92. 
    
93. #Re-run and see if solution changes
    
94. K = 9
    
95. km = kmeans(t(ncidat),centers=K)
    
96. plot(Z[,1],Z[,2],col=km$cluster,type="n")
    
97. text(Z[,1],Z[,2],colnames(ncidat),cex=.75,col=km$cluster)
    
98. cens = km$centers
    
99. points(cens%*%V[,1],cens%*%V[,2],col=1:K,pch=16,cex=3)
    
100. 
     
101. #try different K
     
102. K = 5
     
103. km = kmeans(t(ncidat),centers=K)
     
104. plot(Z[,1],Z[,2],col=km$cluster,type="n")
     
105. text(Z[,1],Z[,2],colnames(ncidat),cex=.75,col=km$cluster)
     
106. cens = km$centers
     
107. points(cens%*%V[,1],cens%*%V[,2],col=1:K,pch=16,cex=3)

复制代码

1. #################
   
2. #lecture 19 code
   
3. #hierarchical clustering
   
4. 
   
5. require("ISLR")
   
6. 
   
7. ncidat = t(NCI60$data)
   
8. colnames(ncidat) = NCI60$labs
   
9. 
   
10. dim(ncidat)
    
11. unique(colnames(ncidat))
    
12. 
    
13. #complete linakge - Euclidean distance
    
14. cols = as.numeric(as.factor(colnames(ncidat)))
    
15. Dmat = dist(t(ncidat))
    
16. com.hclust = hclust(Dmat,method="complete")
    
17. plot(com.hclust,cex=.7,main="Complete Linkage")
    
18. 
    
19. #single linakge
    
20. dev.new()
    
21. sing.hclust = hclust(Dmat,method="single")
    
22. plot(sing.hclust,cex=.7,main="Single Linkage")
    
23. 
    
24. #average linakge
    
25. dev.new()
    
26. ave.hclust = hclust(Dmat,method="average")
    
27. plot(ave.hclust,cex=.7,main="Average Linkage")
    
28. 
    
29. #Ward's linakge
    
30. dev.new()
    
31. ward.hclust = hclust(Dmat,method="ward.D")
    
32. plot(ward.hclust,cex=.7,main="Ward's Linkage")
    
33. 
    
34. #complete linkage with different distances
    
35. dev.new()
    
36. Dmat = dist(t(ncidat),method="manhattan") #L1 distance
    
37. com.hclust = hclust(Dmat,method="complete")
    
38. plot(com.hclust,cex=.7,main="Complete Linkage - L1 Dist")
    
39. 
    
40. 
    
41. ##########
    
42. #Biclustering - Cluster Heatmap
    
43. 
    
44. require("ISLR")
    
45. ncidat = t(NCI60$data)
    
46. colnames(ncidat) = NCI60$labs
    
47. 
    
48. #filter genes using PCA
    
49. X = t(scale(t(ncidat),center=TRUE,scale=FALSE))
    
50. sv = svd(t(X));
    
51. V = sv$v
    
52. 
    
53. #PC loadings - visualize data by limiting to top genes in magnitude in the PC loadings
    
54. aa = grep("grey",colors())
    
55. bb = grep("green",colors())
    
56. cc = grep("red",colors())
    
57. gcol2 = colors()[c(aa[1:30],bb[1:20],rep(cc,2))]
    
58. 
    
59. j = 2
    
60. ord = order(abs(V[,j]),decreasing=TRUE)
    
61. x = as.matrix(X[ord[1:250],])
    
62. 
    
63. #cluster heatmap - uses Ward's linkage (complete is default)
    
64. heatmap(x,col=gcol2,hclustfun=function(x)hclust(x,method="ward.D"))
    
65. 
    
66. ######################################################

复制代码

1. ##################
   
2. #lecture 20 code  - graphical models
   
3. 
   
4. require("diagram")
   
5. require("glasso")
   
6. 
   
7. data = as.matrix(read.table("sachs_data.txt"))
   
8. labs = c("praf",
   
9. "pmek",
   
10. "plcg",
    
11. "PIP2",
    
12. "PIP3",
    
13. "p44/42",
    
14. "pakts473",
    
15. "PKA",
    
16. "PKC",
    
17. "P38",
    
18. "pjnk")
    
19. xc = scale(log(data-min(data)+1),center=TRUE,scale=FALSE)
    
20. n = nrow(data)
    
21. 
    
22. 
    
23. par(mfrow=c(2,3))
    
24. lams = c(.0001, .00025, .0005, .00075, .001, .0025)
    
25. for(i in 1:6){
    
26.   gl = glasso(t(xc)%*%xc/n,rho=lams[i])
    
27.   plotweb(gl$wi,legend=FALSE,length=0,names=labs,main=lams[i])
    
28. }
    
29.   

复制代码

1. ###############
   
2. #lecture 22 - CART
   
3. #classification tree example
   
4. 
   
5. require(rpart)
   
6. 
   
7. #reading in data
   
8. data = read.csv("spam_dat.csv",header=FALSE)
   
9. ndat = read.delim("spam_vars.txt",header=FALSE)
   
10. 
    
11. #parsing variable names
    
12. nams = NULL
    
13. for(i in 1:nrow(ndat))
    
14. {
    
15.   vec = strsplit(as.character(ndat[i,]),split="_")
    
16.   for(j in 1:length(vec[[1]]))
    
17.     {
    
18.       if(length(grep(":",vec[[1]][j]))>0)
    
19.         {
    
20.           vars = strsplit(vec[[1]][j],split=":")
    
21.           nams = c(nams,vars[[1]][1])
    
22.         }
    
23.     }
    
24. }
    
25. 
    
26. Y = data[,58]
    
27. n = length(Y)
    
28. sum(Y)/n
    
29. X = as.matrix(log(1 + data[,1:57]))
    
30. colnames(X) = nams
    
31. X = scale(X)/sqrt(n-1)
    
32. Yf = as.factor(Y)
    
33. levels(Yf) = c("not","spam")
    
34. dat = data.frame(Yf,X)
    
35. 
    
36. 
    
37. #high CP
    
38. ans = rpart(Yf~.,data=dat,method="class",xval=5,cp=.1)
    
39. plot(ans,margin=.05)
    
40. text(ans,use.n=TRUE)
    
41. 
    
42. summary(ans)
    
43. 
    
44. plotcp(ans)
    
45. 
    
46. nams[49:54]
    
47. 
    
48. 
    
49. #medium cp
    
50. ans = rpart(Yf~.,data=dat,method="class",xval=5,cp=.01)
    
51. plot(ans,margin=.05)
    
52. text(ans,use.n=TRUE)
    
53. 
    
54. summary(ans)
    
55. 
    
56. plotcp(ans)
    
57. 
    
58. 
    
59. #low cp
    
60. ans = rpart(Yf~.,data=dat,method="class",xval=5,cp=.005)
    
61. plot(ans,margin=.05)
    
62. text(ans,use.n=TRUE,cex=.75)
    
63. 
    
64. summary(ans)
    
65. 
    
66. plotcp(ans)
    
67. 
    
68. 
    
69. ###############
    
70. #splitting data into training and testing
    
71. 
    
72. 
    
73. trind = sample(1:n,floor(n/2))
    
74. trdat = data.frame(Yf[trind],X[trind,])
    
75. names(trdat)[1] = "Yf"
    
76. tsdat = data.frame(Yf[-trind],X[-trind,])
    
77. names(tsdat)[1] = "Yf"
    
78. 
    
79. 
    
80. ans = rpart(Yf~.,data=tsdat,method="class",xval=5,cp=.005)
    
81. 
    
82. #trees
    
83. 
    
84. ans = rpart(Yf~.,data=trdat,method="class",xval=5,cp=.005)
    
85. trerr = sum(predict(ans,newdata=trdat,type="class")!=trdat$Yf)/nrow(trdat)
    
86. tserr = sum(predict(ans,newdata=tsdat,type="class")!=tsdat$Yf)/nrow(tsdat)
    
87. c(trerr, tserr)

复制代码

1. ####################
   
2. #lecture 23 - boosting
   
3. 
   
4. require(rpart)
   
5. 
   
6. #reading in data
   
7. data = read.csv("spam_dat.csv",header=FALSE)
   
8. ndat = read.delim("spam_vars.txt",header=FALSE)
   
9. 
   
10. #parsing variable names
    
11. nams = NULL
    
12. for(i in 1:nrow(ndat))
    
13. {
    
14.   vec = strsplit(as.character(ndat[i,]),split="_")
    
15.   for(j in 1:length(vec[[1]]))
    
16.     {
    
17.       if(length(grep(":",vec[[1]][j]))>0)
    
18.         {
    
19.           vars = strsplit(vec[[1]][j],split=":")
    
20.           nams = c(nams,vars[[1]][1])
    
21.         }
    
22.     }
    
23. }
    
24. 
    
25. Y = data[,58]
    
26. n = length(Y)
    
27. sum(Y)/n
    
28. X = as.matrix(log(1 + data[,1:57]))
    
29. colnames(X) = nams
    
30. X = scale(X)/sqrt(n-1)
    
31. ind = sample(1:nrow(X),floor(nrow(X)*.3))
    
32. xts = X[ind,]; xtr = X[-ind,]
    
33. yts = Y[ind]; ytr = Y[-ind]
    
34. trdat = data.frame(ytr,xtr)
    
35. tsdat = data.frame(yts,xts)
    
36. 
    
37. 
    
38. #trees
    
39. fit1 = rpart(ytr~.,data=trdat,method="class",xval=5,cp=.01)
    
40. trpred = apply(predict(fit1),1,which.max) - 1
    
41. tree.trerr = sum(abs(ytr - trpred))/length(ytr)
    
42. tspred = apply(predict(fit1,newdata=data.frame(xts)),1,which.max) - 1
    
43. tree.tserr = sum(abs(yts - tspred))/length(yts)
    
44. tree.err = c(tree.trerr, tree.tserr)
    
45. tree.err
    
46. 
    
47. 
    
48. #logistic regression
    
49. fit2 = glm(ytr~.,data=trdat,family="binomial")
    
50. trpred = (sign(predict(fit2)) + 1)/2
    
51. logit.trerr = sum(abs(ytr - trpred))/length(ytr)
    
52. tspred = (sign(predict(fit2,newdata=data.frame(xts))) + 1)/2
    
53. logit.tserr = sum(abs(yts - tspred))/length(yts)
    
54. logit.err = c(logit.trerr, logit.tserr)
    
55. logit.err
    
56. 
    
57. 
    
58. #######boosting
    
59. require(ada)
    
60. 
    
61. #adaboost
    
62. fit3 = ada(x=xtr,y=ytr,loss="exponential",type="discrete",iter=150)
    
63. ab.trerr = (fit3$confusion[2,1] + fit3$confusion[1,2])/length(ytr)
    
64. plot(fit3)
    
65. tspred = round(predict(fit3,newdata=data.frame(xts),type="probs"))
    
66. ab.tserr = sum(abs(yts - (apply(tspred,1,which.max)-1)))/length(yts)
    
67. ab.err = c(ab.trerr, ab.tserr)
    
68. ab.err
    
69. 
    
70. #logitboost
    
71. fit4 = ada(x=xtr,y=ytr,loss="logistic",type="real",iter=150)
    
72. lb.trerr = (fit4$confusion[2,1] + fit4$confusion[1,2])/length(ytr)
    
73. plot(fit4)
    
74. tspred = round(predict(fit4,newdata=data.frame(xts),type="probs"))
    
75. lb.tserr = sum(abs(yts - (apply(tspred,1,which.max)-1)))/length(yts)
    
76. lb.err = c(lb.trerr, lb.tserr)
    
77. lb.err
    
78. 
    
79. errmat = cbind(tree.err,logit.err,ab.err,lb.err)
    
80. errmat

复制代码