logit.R-20170901

library(caret) 
data=as.matrix(oadata2)
colnames(data)=c('y','x1','x2','x17',paste("x",3:16,sep=""))

c=5000 #验证c次
set.seed(5000)
r1=matrix(0,c,1)
r2=matrix(0,c,1)
p1=matrix(0,12,c)

for (k in 1:c) {
  index <-createDataPartition(data[,1], time=1, p=0.8, list=F)
  train=data[index, ]
  test=data[-index, ]
  colnames(train)=c('y','x1','x2','x17',paste("x",3:16,sep=""))
  colnames(test)=c('y','x1','x2','x17',paste("x",3:16,sep=""))

library(mlogit)  
#logit
lm <- glm(y~0+x1+x2+x17+as.factor(x3)+as.factor(x4)+as.factor(x5)+as.factor(x6)+as.factor(x7)+as.factor(x8)
         +as.factor(x9)+as.factor(x10)+as.factor(x11)+as.factor(x12)+as.factor(x13)+as.factor(x14)
         +as.factor(x15)+as.factor(x16),data=data.frame(train1))
summary(lm)
predict(lm,newx=test[,2:18])

#logit, Elastic Net, 0<alpha<1
library(MASS)
library(glmnet)
l=cv.glmnet(x=cbind(train[,'x1'],train[,'x2'],train[,'x17'],as.factor(train[,'x3']), 
         as.factor(train[,'x4']),as.factor(train[,'x5']),as.factor(train[,'x6']),
         as.factor(train[,'x7']),as.factor(train[,'x8']),as.factor(train[,'x9']),
         as.factor(train[,'x10']),as.factor(train[,'x11']),as.factor(train[,'x12']),
         as.factor(train[,'x13']),as.factor(train[,'x14']),as.factor(train[,'x15']),
         as.factor(train[,'x16'])),y=as.factor(train[,1]), family="multinomial",
         alpha=0.5, intercept=FALSE,type.measure="mse")
l=cv.glmnet(x=train[,2:18],y=as.factor(train[,1]), family="multinomial",
            alpha=0.5,intercept=FALSE, type.measure="mse")
print(l)
l$lambda.min
l$lambda.1se
coef(l,s=l$lambda.min)

coefficients=coef(l,s=l$lambda.min)
Active.Index<-which(coefficients!=0)     #系数不为0的特征索引
Active.coefficients<-coefficients[Active.Index]   #系数不为0的特征系数值

predict(l,s=l$lambda.min,newx=test[,2:18])
l.ytest=predict(l,s=l$lambda.min,newx=test[,2:18])
mean((l.ytest-test[,1])^2) 
c(0,1,2,3,4)
l.ysum=data.frame(rowSums(l.ytest*c(0,1,2,3,4)))

c(l$lambda.min, l$lambda.1se)
predict(l, newx=test[,2:18], type="response",s="lambda.1se")
predict(l, newx=test[,2:18], type="response",s="lambda.min")

# CV for 11 alpha value,取alpha的值
mse=matrix(0,10,1)
for (i in 0:10) {
  li=cv.glmnet(x=train[,2:18],y=as.factor(train[,1]), family="multinomial", alpha=i/10,
              intercept=FALSE, type.measure="class")
  assign(paste("l", i, sep=""),li)
  yhat=predict(li, s=l0$lambda.1se, newx=test[,2:18])
  assign(paste("yhat", i, sep=""),yhat)
  mse[i]=mean((test[,1] - yhat)^2)
}
print(mse) #取使mse最小的alpha
# Plot Solution Paths
par(mfrow=c(3,1))
plot(l10, main="LASSO")
plot(l0, main="Ridge")
plot(l5, main="Elastic Net")

l.lasso <- glmnet(train[,2:18], train[,1], family="gaussian", alpha=1)
l.ridge <- glmnet(train[,2:18], train[,1], family="gaussian", alpha=0)
l.elnet <- glmnet(train[,2:18], train[,1], family="gaussian", alpha=.5)

# 10-fold Cross validation for each alpha = 0, 0.1, ... , 0.9, 1.0
l.lasso.cv <- cv.glmnet(train[,2:18], train[,1], type.measure="mse", alpha=1, 
                          family="gaussian")
l.ridge.cv <- cv.glmnet(x.train, y.train, type.measure="mse", alpha=0,
                          family="gaussian")
l.elnet.cv <- cv.glmnet(x.train, y.train, type.measure="mse", alpha=.5,
                          family="gaussian")
# Plot solution paths:
par(mfrow=c(3,2))
# For plotting options, type '?plot.glmnet' in R console
plot(l.lasso, xvar="lambda")
plot(l10, main="LASSO")

plot(l.ridge, xvar="lambda")
plot(l0, main="Ridge")

plot(l.elnet, xvar="lambda")
plot(l5, main="Elastic Net")



#guassion 
g=glmnet(x=train[,2:18], y=train[,1], family="gaussian")
print(g)
coef(g,s=min(g$lambda))
coef(g, s=c(g$lambda[16],0.5))
plot(g, xvar="lambda", label=TRUE)
predict(g, newx=test[,2:18], s=c(g$lambda[16],0.5))

library(GPfit)
library(lhs)
g <- GP_fit(train[,2:18], train[,1])

# m<-svm(train[,2:18],train[,1])
# #assign(paste("m",k,sep=""),m)
# #summary(m)
# ytest2=predict(m,test[,2:18])
# #assign(paste("ytest2_",k,sep=""),ytest2)
# resi2=abs(ytest2-test[,1])/test[,1]
# r2[k]=mean(resi2) #误差
# #print(mean(resi2))

mean(r1)
#mean(r2)