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

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

Course Description
This course is a survey of statistical learning methods and will cover major techniques and concepts for both supervised and unsupervised learning. Topics covered include penalized regression and classification, support vector machines, kernel methods, model selection, matrix factorizations, graphical models, clustering, boosting, random forests, and ensemble learning. Students will learn how and when to apply statistical learning techniques, their comparative strengths and weaknesses, and how to critically evaluate the performance of learning algorithms. Students completing this course should be able to

• (i) apply basic statistical learning methods to build predictive models or perform exploratory analysis,
• (ii) properly tune and select statistical learning models,
• (iii) correctly assess model fit and error, and
• (iv) build an ensemble of learning algorithms.
  

本帖隐藏的内容

Genevera Allen Machine Learning Rice University.zip (4.45 MB, 需要: 1 个论坛币)

2016-5-29 22:43:04 上传

需要: 1 个论坛币  [购买]

扫码加我 拉你入群

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

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

分享0 收藏3 回帖

关键词：Statistical statistica statistic Learning earning techniques machines learning concepts include

1. #########################
   
2. #lecture 2 code
   
3. #ozone data
   
4. 
   
5. #load in ozone data (taken from "cosso" R package)
   
6. ozone = read.csv("ozone.csv")
   
7. 
   
8. #take a look at the data
   
9. plot(ozone,pch=16,cex=.5)
   
10. 
    
11. #directly fit least squares
    
12. Y = as.numeric(ozone[,1])
    
13. aX = cbind(rep(1,nrow(ozone)),as.matrix(ozone[,-1]))
    
14. 
    
15. #estimated coefficients
    
16. betahat = solve(t(aX)%*%aX)%*%t(aX)%*%Y
    
17. 
    
18. #fitted data
    
19. Yhat = aX%*%betahat
    
20. 
    
21. #hat matrix
    
22. H = aX%*%solve(t(aX)%*%aX)%*%t(aX)
    
23. evH =  eigen(H)
    
24. evH$values[1:10]
    
25. 
    
26. #fit linear model using "lm"
    
27. fit = lm(ozone~.,data=ozone)
    
28. summary(fit)
    
29. 
    
30. cbind(betahat,t(t(fit$coefficients)))
    
31. 
    
32. 
    
33. #diagnostic plots
    
34. layout(matrix(c(1,2,3,4),2,2))
    
35. plot(fit)
    
36. dev.off()

复制代码

1. ###########################
   
2. #lecture 3 code
   
3. #ozone data
   
4. 
   
5. #load in ozone data (taken from "cosso" R package)
   
6. ozone = read.csv("ozone.csv")
   
7. 
   
8. #take a look at the data
   
9. plot(ozone,pch=16,cex=.5)
   
10. 
    
11. 
    
12. #ridge regression
    
13. Y = as.numeric(ozone[,1]); Y = Y - mean(Y)
    
14. X = as.matrix(ozone[,-1]); X = scale(X,center=T,scale=F)
    
15. 
    
16. #ridge regression solution
    
17. lam = nrow(ozone)
    
18. 
    
19. betals = solve(t(X)%*%X)%*%t(X)%*%Y
    
20. betar = solve(t(X)%*%X + diag(rep(lam,8)))%*%t(X)%*%Y
    
21. cbind(betals,betar)
    
22. 
    
23. 
    
24. #ridge regression coefficient paths
    
25. lambdas = exp(seq(log(.01),log(10*nrow(ozone)),l=100))
    
26. betasr = matrix(0,length(lambdas),8)
    
27. for(i in 1:length(lambdas))
    
28. {
    
29.   betasr[i,] = solve(t(X)%*%X + diag(rep(lambdas[i],8)))%*%t(X)%*%Y
    
30. }
    
31. 
    
32. plot(c(1,length(lambdas)),range(betasr),type="n",ylab="Coefficients",xlab="Lambda Index")
    
33. for(j in 1:8)
    
34. {
    
35.   lines(betasr[length(lambdas):1,j],col=j)
    
36. }
    
37. legend(0,20,legend=names(ozone)[2:9],col=1:9,lty=rep(1,9))
    
38. 
    
39. 
    
40. ################################
    
41. #looking at principal components
    
42. svdx = svd(X)
    
43. 
    
44. #scatterplots of samples PCs
    
45. par(mar=c(1,1,1,1))
    
46. layout(matrix(1:25,5,5))
    
47. mycols = rainbow(length(Y))
    
48. orY = order(Y)
    
49. for(i in 1:5)
    
50. {
    
51.   for(j in 1:5)
    
52.     {
    
53.       plot(svdx$u[,i],svdx$u[,j],type="p",pch=16,col=mycols[orY])
    
54.     }
    
55. }
    
56. 
    
57. #amount of variance explained
    
58. varex = 0; cumvarex = 0;
    
59. for(i in 1:8)
    
60. {
    
61.   varex[i] = svdx$d[i]/sum(svdx$d)
    
62.   cumvarex[i] = sum(varex)
    
63. }
    
64. par(mfrow=c(1,2))
    
65. par(mar=c(5,4,4,2))
    
66. barplot(varex,ylab="Amount of Var Explained",xlab="PCs")
    
67. barplot(cumvarex,ylab="Cummulative Var Explained",xlab="PCs")
    
68. 
    
69. 
    
70. #PC direction weights
    
71. par(mfrow=c(3,2))
    
72. par(mar=c(5,4,3,2))
    
73. for(i in 1:6)
    
74.   {
    
75.     barplot(svdx$v[,i],names.arg=names(ozone)[2:9])
    
76.   }
    
77. 
    
78. 
    
79. #ridge paths again
    
80. dev.off()
    
81. plot(c(1,length(lambdas)),range(betasr),type="n",ylab="Coefficients",xlab="Lambda Index")
    
82. for(j in 1:8)
    
83. {
    
84.   lines(betasr[length(lambdas):1,j],col=j)
    
85. }
    
86. legend(0,20,legend=names(ozone)[2:9],col=1:9,lty=rep(1,9))
    
87. 
    
88. 
    
89. 
    
90. #####################
    
91. #least squares on derived inputs
    
92. 
    
93. #PC regression
    
94. betapcr = diag(svdx$d)%*%t(svdx$u)%*%Y
    
95. 
    
96. 
    
97. #PLS regression
    
98. plsfunc = function(x,y)
    
99. {
    
100.   p = ncol(x); n = nrow(x);
     
101.   M = t(x)%*%y
     
102.   Z = NULL; V = NULL; P = NULL;
     
103.   for(k in 1:p)
     
104.     {
     
105.       svdm = svd(M)
     
106.       z = x%*%svdm$u
     
107.       z = z*as.numeric(1/sqrt(t(z)%*%z))
     
108.       V = cbind(V,svdm$u)
     
109.       p = t(x)%*%z/as.numeric(t(z)%*%z)
     
110.       P = cbind(P,p);
     
111.       Z = cbind(Z,z);
     
112.       M = M - P%*%solve(t(P)%*%P)%*%t(P)%*%M;
     
113.     }
     
114.   return(list(Z=Z,V=V))
     
115. }
     
116. 
     
117. plsx = plsfunc(X,Y)
     
118. 
     
119. 
     
120. #scatterplots of pls components
     
121. par(mar=c(1,1,1,1))
     
122. layout(matrix(1:25,5,5))
     
123. mycols = rainbow(length(Y))
     
124. orY = order(Y)
     
125. for(i in 1:5)
     
126. {
     
127.   for(j in 1:5)
     
128.     {
     
129.       plot(plsx$Z[,i],plsx$Z[,j],type="p",pch=16,col=mycols[orY])
     
130.     }
     
131. }
     
132. 
     
133. betapls = t(plsx$Z)%*%Y
     
134. 
     
135. cbind(betapcr,betapls)

复制代码

1. ###########################
   
2. #lecture 4 code
   
3. #ozone data
   
4. #ridge regression, PCA, PC regression & PLS regression
   
5. 
   
6. 
   
7. #load in ozone data (taken from "cosso" R package)
   
8. ozone = read.csv("ozone.csv")
   
9. 
   
10. #take a look at the data
    
11. plot(ozone,pch=16,cex=.5)
    
12. 
    
13. ##################
    
14. #ridge regression
    
15. 
    
16. Y = as.numeric(ozone[,1]); Y = Y - mean(Y)
    
17. X = as.matrix(ozone[,-1]); X = scale(X,center=T,scale=F)
    
18. 
    
19. #ridge regression solution
    
20. lam = nrow(ozone)*.1
    
21. 
    
22. betals = solve(t(X)%*%X)%*%t(X)%*%Y
    
23. betar = solve(t(X)%*%X + diag(rep(lam,8)))%*%t(X)%*%Y
    
24. cbind(betals,betar)
    
25. 
    
26. 
    
27. #ridge regression coefficient paths
    
28. lambdas = exp(seq(log(.01),log(10*nrow(ozone)),l=100))
    
29. betasr = matrix(0,length(lambdas),8)
    
30. for(i in 1:length(lambdas))
    
31. {
    
32.   betasr[i,] = solve(t(X)%*%X + diag(rep(lambdas[i],8)))%*%t(X)%*%Y
    
33. }
    
34. 
    
35. plot(c(1,length(lambdas)),range(betasr),type="n",ylab="Coefficients",xlab="Lambda Index")
    
36. for(j in 1:8)
    
37. {
    
38.   lines(betasr[length(lambdas):1,j],col=j)
    
39. }
    
40. legend(0,20,legend=names(ozone)[2:9],col=1:9,lty=rep(1,9))
    
41. 
    
42. 
    
43. ################################
    
44. #looking at principal components
    
45. svdx = svd(X)
    
46. 
    
47. #scatterplots of samples PCs
    
48. par(mar=c(1,1,1,1))
    
49. layout(matrix(1:25,5,5))
    
50. mycols = rainbow(length(Y))
    
51. orY = order(Y)
    
52. for(i in 1:5)
    
53. {
    
54.   for(j in 1:5)
    
55.     {
    
56.       plot(svdx$u[,i],svdx$u[,j],type="p",pch=16,col=mycols[orY])
    
57.     }
    
58. }
    
59. 
    
60. #amount of variance explained
    
61. varex = 0; cumvarex = 0;
    
62. for(i in 1:8)
    
63. {
    
64.   varex[i] = svdx$d[i]/sum(svdx$d)
    
65.   cumvarex[i] = sum(varex)
    
66. }
    
67. par(mfrow=c(1,2))
    
68. par(mar=c(5,4,4,2))
    
69. barplot(varex,ylab="Amount of Var Explained",xlab="PCs")
    
70. barplot(cumvarex,ylab="Cummulative Var Explained",xlab="PCs")
    
71. 
    
72. 
    
73. #PC direction weights
    
74. par(mfrow=c(3,2))
    
75. par(mar=c(5,4,3,2))
    
76. for(i in 1:6)
    
77.   {
    
78.     barplot(svdx$v[,i],names.arg=names(ozone)[2:9])
    
79.   }
    
80. 
    
81. 
    
82. #ridge paths again
    
83. dev.off()
    
84. plot(c(1,length(lambdas)),range(betasr),type="n",ylab="Coefficients",xlab="Lambda Index")
    
85. for(j in 1:8)
    
86. {
    
87.   lines(betasr[length(lambdas):1,j],col=j)
    
88. }
    
89. legend(0,20,legend=names(ozone)[2:9],col=1:9,lty=rep(1,9))
    
90. 
    
91. 
    
92. 
    
93. #####################
    
94. #least squares on derived inputs
    
95. 
    
96. #PC regression
    
97. betapcr = diag(svdx$d)%*%t(svdx$u)%*%Y
    
98. 
    
99. 
    
100. #PLS regression
     
101. plsfunc = function(x,y)
     
102. {
     
103.   p = ncol(x); n = nrow(x);
     
104.   M = t(x)%*%y
     
105.   Z = NULL; V = NULL; P = NULL;
     
106.   for(k in 1:p)
     
107.     {
     
108.       svdm = svd(M)
     
109.       z = x%*%svdm$u
     
110.       z = z*as.numeric(1/sqrt(t(z)%*%z))
     
111.       V = cbind(V,svdm$u)
     
112.       p = t(x)%*%z/as.numeric(t(z)%*%z)
     
113.       P = cbind(P,p);
     
114.       Z = cbind(Z,z);
     
115.       M = M - P%*%solve(t(P)%*%P)%*%t(P)%*%M;
     
116.     }
     
117.   return(list(Z=Z,V=V))
     
118. }
     
119. 
     
120. plsx = plsfunc(X,Y)
     
121. 
     
122. 
     
123. #scatterplots of pls components
     
124. par(mar=c(1,1,1,1))
     
125. layout(matrix(1:25,5,5))
     
126. mycols = rainbow(length(Y))
     
127. orY = order(Y)
     
128. for(i in 1:5)
     
129. {
     
130.   for(j in 1:5)
     
131.     {
     
132.       plot(plsx$Z[,i],plsx$Z[,j],type="p",pch=16,col=mycols[orY])
     
133.     }
     
134. }
     
135. 
     
136. betapls = t(plsx$Z)%*%Y
     
137. 
     
138. cbind(betapcr,betapls)
     
139. 
     
140. 
     
141. ########################################
     
142. #feature selection in linear models
     
143. 
     
144. #libraries required
     
145. library(MASS)
     
146. library(leaps)
     
147. library(glmnet)
     
148. 
     
149. #####################
     
150. #Algorithmic
     
151. 
     
152. #best subsets selection
     
153. fitbsub = regsubsets(x=ozone[,2:9],y=ozone[,1])
     
154. summary(fitbsub)
     
155. 
     
156. #forward step-wise  - via BIC
     
157. fit0 = lm(ozone~1,data=ozone)
     
158. fitf = stepAIC(fit0,scope=ozone~wind+temp+invHt+press+hum+vis+milPress+invTemp,direction="forward",data=ozone,k=log(nrow(ozone)))
     
159. summary(fitf)
     
160. 
     
161. #backwards step-wise - via BIC
     
162. fit = lm(ozone~.,data=ozone)
     
163. fitb = stepAIC(fit,direction="backward",data=ozone,k=log(nrow(ozone)))
     
164. summary(fitb)

复制代码

感谢分享

1. ###########################
   
2. #lecture 5 code
   
3. #ozone data
   
4. 
   
5. ozone = read.csv("ozone.csv")
   
6. 
   
7. ##################
   
8. #ridge regression
   
9. 
   
10. Y = as.numeric(ozone[,1]); Y = Y - mean(Y)
    
11. X = as.matrix(ozone[,-1]); X = scale(X,center=T,scale=F)
    
12. 
    
13. #ridge regression coefficient paths
    
14. lambdas = exp(seq(log(.01),log(10*nrow(ozone)),l=100))
    
15. betasr = matrix(0,length(lambdas),8)
    
16. for(i in 1:length(lambdas))
    
17. {
    
18.   betasr[i,] = solve(t(X)%*%X + diag(rep(lambdas[i],8)))%*%t(X)%*%Y
    
19. }
    
20. 
    
21. plot(c(1,length(lambdas)),range(betasr),type="n",ylab="Coefficients",xlab="Lambda Index")
    
22. for(j in 1:8)
    
23. {
    
24.   lines(betasr[length(lambdas):1,j],col=j)
    
25. }
    
26. legend(0,20,legend=names(ozone)[2:9],col=1:9,lty=rep(1,9))
    
27. 
    
28. plot(ozone)
    
29. 
    
30. ########################################
    
31. #feature selection in linear models
    
32. 
    
33. #libraries required
    
34. library(MASS)
    
35. library(leaps)
    
36. library(glmnet)
    
37. 
    
38. #####################
    
39. #Algorithmic
    
40. 
    
41. #best subsets selection
    
42. fitbsub = regsubsets(x=ozone[,2:9],y=ozone[,1])
    
43. summary(fitbsub)
    
44. 
    
45. #forward step-wise  - via BIC
    
46. fit0 = lm(ozone~1,data=ozone)
    
47. fitf = stepAIC(fit0,scope=ozone~wind+temp+invHt+press+hum+vis+milPress+invTemp,direction="forward",data=ozone,k=log(nrow(ozone)))
    
48. summary(fitf)
    
49. 
    
50. #backwards step-wise - via BIC
    
51. fit = lm(ozone~.,data=ozone)
    
52. fitb = stepAIC(fit,direction="backward",data=ozone,k=log(nrow(ozone)))
    
53. summary(fitb)
    
54. 
    
55. #######################
    
56. #L1 Regularization
    
57. 
    
58. fit0 = lm(ozone~.-1,data=ozone)
    
59. Y = as.numeric(ozone[,1])
    
60. X = as.matrix(ozone[,-1])
    
61. 
    
62. lam = 1
    
63. fitl = glmnet(x=X,y=Y,family="gaussian",lambda=lam,alpha=1)
    
64. cbind(fit0$coef,as.matrix(fitl$beta))
    
65. 
    
66. #lasso paths
    
67. fitl = glmnet(x=X,y=Y,family="gaussian",alpha=1)
    
68. plot(fitl,col=1:8)
    
69. legend(0,19,legend=names(ozone)[2:9],col=1:8,lty=rep(1,8),cex=.8)

复制代码

1. ###########################
   
2. #lecture 6 code
   
3. #ozone data
   
4. #comparison: sparse regression methods
   
5. 
   
6. #libraries required
   
7. library(glmnet)
   
8. library(ncvreg)
   
9. 
   
10. #load in ozone data
    
11. ozone = read.csv("ozone.csv")
    
12. Y = as.numeric(ozone[,1]); Y = Y - mean(Y)
    
13. X = as.matrix(ozone[,-1]); X = scale(X,center=T,scale=F)
    
14. 
    
15. #Lasso
    
16. lam = 1
    
17. fit0 = lm(Y~X-1)
    
18. fitl = glmnet(x=X,y=Y,family="gaussian",lambda=lam,alpha=1)
    
19. cbind(fit0$coef,as.matrix(fitl$beta))
    
20. 
    
21. #Lasso paths
    
22. fitl = glmnet(x=X,y=Y,family="gaussian",alpha=1)
    
23. plot(fitl,col=1:8)
    
24. legend(0,19,legend=names(ozone)[2:9],col=1:8,lty=rep(1,8),cex=.8)
    
25. 
    
26. 
    
27. 
    
28. ###############################
    
29. #least squares, lasso, adaptive lasso, SCAD, ridge, elastic net, MC+
    
30. lam = 1
    
31. 
    
32. betals = solve(t(X)%*%X)%*%t(X)%*%Y
    
33. betar = solve(t(X)%*%X + diag(rep(lam/2*nrow(ozone),8)))%*%t(X)%*%Y
    
34. fitl = glmnet(x=X,y=Y,family="gaussian",lambda=lam,alpha=1)
    
35. fital = glmnet(x=X,y=Y,family="gaussian",lambda=lam,alpha=1,penalty.factor=1/abs(betals))
    
36. fitel = glmnet(x=X,y=Y,family="gaussian",lambda=lam,alpha=.5)
    
37. fitscad = ncvreg(X,Y,family="gaussian",penalty="SCAD",lambda=lam)
    
38. fitmcp = ncvreg(X,Y,family="gaussian",penalty="MCP",lambda=lam)
    
39. mat = cbind(betals,betar,as.matrix(fitl$beta),as.matrix(fital$beta),as.matrix(fitel$beta),fitscad$beta[-1],fitmcp$beta[-1])
    
40. colnames(mat) = c("LS","Ridge","Lasso","A-Lasso","EL","SCAD","MC+")
    
41. mat
    
42. 
    
43. 
    
44. #############################
    
45. #compare ridge, lasso, elastic net & SCAD regualrization paths
    
46. 
    
47. par(mfrow=c(2,3))
    
48. par(mar=c(5,4,3,2))
    
49. betals = solve(t(X)%*%X)%*%t(X)%*%Y
    
50. lambdas = exp(seq(log(.01),log(100*nrow(ozone)),l=100))
    
51. betasr = matrix(0,length(lambdas),8)
    
52. for(i in 1:length(lambdas))
    
53. {
    
54.   betasr[i,] = solve(t(X)%*%X + diag(rep(lambdas[i],8)))%*%t(X)%*%Y
    
55. }
    
56. plot(c(1,length(lambdas)),range(betals),type="n",ylab="Coefficients",xlab="Lambda Index",main="Ridge")
    
57. for(j in 1:8)
    
58. {
    
59.   lines(betasr[length(lambdas):1,j],col=j)
    
60. }
    
61. legend(0,20,legend=names(ozone)[2:9],col=1:9,lty=rep(1,9),cex=.75)
    
62. 
    
63. fitl = glmnet(x=X,y=Y,family="gaussian",alpha=1)
    
64. plot(fitl,col=1:8,main="Lasso")
    
65. legend(0,20,legend=names(ozone)[2:9],col=1:8,lty=rep(1,8),cex=.75)
    
66. 
    
67. fitel = glmnet(x=X,y=Y,family="gaussian",alpha=.5)
    
68. plot(fitel,col=1:8,main="EL alpha=.5")
    
69. legend(0,20,legend=names(ozone)[2:9],col=1:8,lty=rep(1,8),cex=.75)
    
70. 
    
71. fitel = glmnet(x=X,y=Y,family="gaussian",alpha=.25)
    
72. plot(fitel,col=1:8,main="EL alpha=.25")
    
73. legend(0,20,legend=names(ozone)[2:9],col=1:8,lty=rep(1,8),cex=.75)
    
74. 
    
75. fitscad = ncvreg(X,Y,family="gaussian",penalty="SCAD")
    
76. plot(fitscad,col=1:8,main="SCAD",shade=F)
    
77. legend(6,30,legend=names(ozone)[2:9],col=1:8,lty=rep(1,8),cex=.75)
    
78. 
    
79. fitmcp = ncvreg(X,Y,family="gaussian",penalty="MCP")
    
80. plot(fitmcp,col=1:8,main="MC+",shade=F)
    
81. legend(6,30,legend=names(ozone)[2:9],col=1:8,lty=rep(1,8),cex=.75)

复制代码

1. %%%%%%%%%%%%%%%%%%%%%%%
   
2. %lecture 8 code
   
3. %LDA
   
4. 
   
5. dat = dlmread('zip.train');
   
6. x = dat(:,2:end);
   
7. [n,p] = size(x);
   
8. 
   
9. imshow(flipud(rot90(reshape(dat(1,2:end),16,16))))
   
10. 
    
11. 
    
12. nclass = length(unique(dat(:,1)));
    
13. Y = zeros(n,nclass);
    
14. nK = 0;
    
15. for i=1:nclass
    
16.     Y(dat(:,1)==(i-1),i) = 1;
    
17.     nK(i) = sum(Y(:,i));
    
18. end
    
19. 
    
20. 
    
21. X = x - ones(n,1)*mean(x);
    
22. mvecs = zeros(p,nclass);
    
23. for i=1:nclass
    
24.     mvecs(:,i) = mean(X(Y(:,i)==1,:));
    
25. end
    
26. 
    
27. SigB = zeros(p,p);
    
28. for i=1:nclass
    
29.     SigB = SigB + nK(i)*mvecs(:,i)*mvecs(:,i)';
    
30. end
    
31. 
    
32. SigT = X'*X;
    
33. SigW = SigT - SigB;
    
34. 
    
35. [V,D] = eig(SigB,SigW);
    
36. 
    
37. 
    
38. scatter(X*V(:,1),X*V(:,2),'SizeData',1)
    
39. axis tight
    
40. text(X*V(:,1),X*V(:,2),num2str(dat(:,1)))
    
41. 
    
42. 
    
43. scatter(X*V(:,3),X*V(:,4),'SizeData',1)
    
44. axis tight
    
45. text(X*V(:,3),X*V(:,4),num2str(dat(:,1)))
    
46. 
    
47. 
    
48. scatter(X*V(:,5),X*V(:,6),'SizeData',1)
    
49. axis tight
    
50. text(X*V(:,5),X*V(:,6),num2str(dat(:,1)))
    
51. 
    
52. 
    
53. [classNB,errNB] = classify(X,X,dat(:,1),'diaglinear');
    
54. [classLDA,errLDA] = classify(X,X,dat(:,1),'linear');
    
55. [classGLDA,errGLDA] = classify(X*V(:,1:10),X*V(:,1:10),dat(:,1),'diaglinear');
    
56. 
    
57. [errNB errLDA  errGLDA]
    
58. 
    
59. 
    
60. tdat = dlmread('zip.test');
    
61. [tn,tp] = size(tdat);
    
62. tx = tdat(:,2:end);
    
63. tX = tx - ones(size(tdat,1),1)*mean(tx);
    
64. 
    
65. [classNB,errNB] = classify(tX,X,dat(:,1),'diaglinear');
    
66. [classLDA,errLDA] = classify(tX,X,dat(:,1),'linear');
    
67. [classGLDA,errGLDA] = classify(tX*V(:,1:10),X*V(:,1:10),dat(:,1),'diaglinear');
    
68. 
    
69. [sum(tdat(:,1)~=classNB)/tn,  sum(tdat(:,1)~=classLDA)/tn,  sum(tdat(:,1)~=classGLDA)/tn]

复制代码

1. ###############
   
2. #lecture 9 code
   
3. 
   
4. #univariate logistic regression - simulated example
   
5. #note: change  coefficients to understand logistic function
   
6. n = 250; p = 1;
   
7. x = matrix(rnorm(n*p),n,p)
   
8. beta0 = 0
   
9. beta = 1
   
10. eps = matrix(rnorm(n),n,1)
    
11. probs = exp(beta0 + x*beta + eps)/(1 + exp(beta0 + x*beta + eps))
    
12. Y = as.numeric(probs>.5)
    
13. 
    
14. #fitting logistic
    
15. fit = glm(Y ~ x,family="binomial")
    
16. summary(fit)
    
17. plot(x,Y)
    
18. xs = seq(min(x),max(x),l=1000)
    
19. pihat = exp(fit$coefficients[1] + xs*fit$coefficients[2])/(1 + exp(fit$coefficients[1] + xs*fit$coefficients[2]))
    
20. lines(xs,pihat)
    
21. lines(c(min(x),max(x)),c(.5,.5),lty=2)
    
22. 
    
23. 
    
24. ###############
    
25. #looking at spam data
    
26. 
    
27. #reading in data
    
28. data = read.csv("spam_dat.csv",header=FALSE)
    
29. ndat = read.delim("spam_vars.txt",header=FALSE)
    
30. 
    
31. #parsing variable names
    
32. nams = NULL
    
33. for(i in 1:nrow(ndat))
    
34. {
    
35.   vec = strsplit(as.character(ndat[i,]),split="_")
    
36.   for(j in 1:length(vec[[1]]))
    
37.     {
    
38.       if(length(grep(":",vec[[1]][j]))>0)
    
39.         {
    
40.           vars = strsplit(vec[[1]][j],split=":")
    
41.           nams = c(nams,vars[[1]][1])
    
42.         }
    
43.     }
    
44. }
    
45. 
    
46. Y = data[,58]
    
47. n = length(Y)
    
48. sum(Y)/n
    
49. X = as.matrix(log(1 + data[,1:57]))
    
50. colnames(X) = nams
    
51. X = scale(X)/sqrt(n-1)
    
52. dat = data.frame(Y,X)
    
53. 
    
54. #taking a subset of data
    
55. sdat = data.frame(Y,dat$george,dat$meeting,dat$total,dat$re,dat$edu,dat$free,dat$your)
    
56. plot(sdat,pch=16,cex=.5)
    
57. 
    
58. #########
    
59. #fit logistic model
    
60. 
    
61. fits = glm(Y~.,data=sdat,family="binomial",maxit=50)
    
62. summary(fits)
    
63. 
    
64. #individual variable total
    
65. fit = glm(Y~dat.total,data=sdat,family="binomial")
    
66. summary(fit)
    
67. plot(sdat$dat.total,sdat$Y)
    
68. xs = seq(min(sdat$dat.total),max(sdat$dat.total),l=1000)
    
69. pihat = exp(fit$coefficients[1] + xs*fit$coefficients[2])/(1 + exp(fit$coefficients[1] + xs*fit$coefficients[2]))
    
70. lines(xs,pihat)
    
71. lines(c(min(xs),max(xs)),c(.5,.5),lty=2)
    
72. 
    
73. 
    
74. #individual variable george
    
75. fit = glm(Y~dat.george,data=sdat,family="binomial")
    
76. summary(fit)
    
77. plot(sdat$dat.george,sdat$Y)
    
78. xs = seq(min(sdat$dat.george),max(sdat$dat.george),l=1000)
    
79. pihat = exp(fit$coefficients[1] + xs*fit$coefficients[2])/(1 + exp(fit$coefficients[1] + xs*fit$coefficients[2]))
    
80. lines(xs,pihat)
    
81. lines(c(min(xs),max(xs)),c(.5,.5),lty=2)
    
82. 
    
83. 
    
84. #individual variable free
    
85. fit = glm(Y~dat.free,data=sdat,family="binomial")
    
86. summary(fit)
    
87. plot(sdat$dat.free,sdat$Y)
    
88. xs = seq(min(sdat$dat.free),max(sdat$dat.free),l=1000)
    
89. pihat = exp(fit$coefficients[1] + xs*fit$coefficients[2])/(1 + exp(fit$coefficients[1] + xs*fit$coefficients[2]))
    
90. lines(xs,pihat)
    
91. lines(c(min(xs),max(xs)),c(.5,.5),lty=2)

复制代码

1. ###############
   
2. #lecture 9 code
   
3. #classifying spam emails - HP data set
   
4. 
   
5. library(glmnet)
   
6. library(ncvreg)
   
7. 
   
8. 
   
9. #reading in data
   
10. data = read.csv("spam_dat.csv",header=FALSE)
    
11. ndat = read.delim("spam_vars.txt",header=FALSE)
    
12. 
    
13. #parsing variable names
    
14. nams = NULL
    
15. for(i in 1:nrow(ndat))
    
16. {
    
17.   vec = strsplit(as.character(ndat[i,]),split="_")
    
18.   for(j in 1:length(vec[[1]]))
    
19.     {
    
20.       if(length(grep(":",vec[[1]][j]))>0)
    
21.         {
    
22.           vars = strsplit(vec[[1]][j],split=":")
    
23.           nams = c(nams,vars[[1]][1])
    
24.         }
    
25.     }
    
26. }
    
27. 
    
28. Y = data[,58]
    
29. n = length(Y)
    
30. sum(Y)/n
    
31. X = as.matrix(log(1 + data[,1:57]))
    
32. colnames(X) = nams
    
33. X = scale(X)/sqrt(n-1)
    
34. dat = data.frame(Y,X)
    
35. 
    
36. #taking a subset of data
    
37. sdat = data.frame(Y,dat$george,dat$meeting,dat$total,dat$re,dat$edu,dat$free,dat$your)
    
38. plot(sdat,pch=16,cex=.5)
    
39. 
    
40. 
    
41. #########
    
42. #fit logistic model
    
43. 
    
44. fits = glm(Y~.,data=sdat,family="binomial",maxit=50)
    
45. summary(fits)
    
46. 
    
47. 
    
48. ###########
    
49. #sparse logistic regression
    
50. 
    
51. #penalized logistic - coefficients
    
52. lam = .1
    
53. 
    
54. Xs = as.matrix(sdat[,2:8])
    
55. sfitl = glmnet(x=Xs,y=Y,family="binomial",lambda=lam,alpha=1)
    
56. sfitr = glmnet(x=Xs,y=Y,family="binomial",lambda=lam,alpha=0)
    
57. sfitel = glmnet(x=Xs,y=Y,family="binomial",lambda=lam,alpha=.5)
    
58. sfitscad = ncvreg(Xs,Y,family="binomial",penalty="SCAD",lambda=lam)
    
59. mat = cbind(fits$coefficients[-1],as.matrix(sfitl$beta),as.matrix(sfitr$beta),as.matrix(sfitel$beta),sfitscad$beta[-1])
    
60. colnames(mat) = c("Logistic","Lasso","Ridge","EL","SCAD")
    
61. mat
    
62. 
    
63. #penalized logistic - regularization paths
    
64. sfitl = glmnet(x=Xs,y=Y,family="binomial",alpha=1)
    
65. sfitr = glmnet(x=Xs,y=Y,family="binomial",alpha=0)
    
66. sfitel = glmnet(x=Xs,y=Y,family="binomial",alpha=.5)
    
67. sfitscad = ncvreg(Xs,Y,family="binomial",,penalty="SCAD")
    
68. 
    
69. par(mfrow=c(2,2))
    
70. plot(sfitl,col=1:7,main="L1")
    
71. legend(0,-100,legend=names(sdat)[2:8],col=1:7,lty=rep(1,7),cex=.75)
    
72. 
    
73. plot(sfitr,col=1:7,main="Ridge")
    
74. legend(0,-10,legend=names(sdat)[2:8],col=1:7,lty=rep(1,7),cex=.75)
    
75. 
    
76. plot(sfitel,col=1:7,main="Elastic Net")
    
77. legend(0,-100,legend=names(sdat)[2:8],col=1:7,lty=rep(1,7),cex=.75)
    
78. 
    
79. plot(sfitscad,col=1:7,main="SCAD",shade=F)
    
80. legend(.22,-100,legend=names(sdat)[2:8],col=1:7,lty=rep(1,7),cex=.75)
    
81. 
    
82. 
    
83. #note: this takes a while to run
    
84. #full data L1 regularizaiton paths
    
85. fit1 = glmnet(x=X,y=Y,family="binomial")
    
86. plot(fit1)

复制代码