HIT机器学习实验二逻辑回归R语言参考代码

##一键清空
rm(list=ls())
##Sigmoid函数
Sigmoid<-function(z){
  return (1/(1+exp(-z)))
}
##计算p(y=1|x)
P1<-function(x,w){
  return (1-Sigmoid(t(w)%*%x))
}
##计算梯度函数
l1<-function(x,w,d,n,y){
  i<-0
  h<-cbind(rep(0,n+1))
  repeat{
    i<-i+1
    h<-h-x[,i]%*%(y[i,]-P1(x[,i],w))/d
    if(i==d){
      break
    }
  }
  return (h)
}
##计算梯度函数（带惩罚项）
L1<-function(x,w,d,n,y,k){
  i<-0
  h<-cbind(rep(0,n+1))
  repeat{
    i<-i+1
    h<-h-x[,i]%*%(y[i,]-P1(x[,i],w))/d
    if(i==d){
      break
    }
  }
  return ((h+k*w))
}
##生成海森矩阵函数
l2<-function(x,w,d,n){
  i<-0
  h<-matrix(0,nrow=n+1,ncol=n+1)
  repeat{
    i<-i+1
    h<-h+x[,i]%*%t(x[,i])%*%(P1(x,w)[1,1]*(1-P1(x,w))[1,1]*diag(n+1))/d
    if(i==d){
      break
    }
  }
  return (h)
}
##生成海森矩阵函数（带惩罚项）
L2<-function(x,w,d,n,k){
  i<-0
  h<-matrix(0,nrow=n+1,ncol=n+1)
  repeat{
    i<-i+1
    h<-h+x[,i]%*%t(x[,i])%*%(P1(x,w)[1,1]*(1-P1(x,w))[1,1]*diag(n+1))/d
    if(i==d){
      break
    }
  }
  return (h+k*diag(n+1))
}
##牛顿法优化函数
Newt<-function(x,w,d,n,y,esp){
  i<-0
  repeat{
    i<-i+1
    if(qr(w)$rank!=n+1&w[1,1]!=0){
      break
    }
    w<-w-solve(l2(x,w,d,n))%*%l1(x,w,d,n,y)
    print(w)
    if(i==10000){
      break
    }
  }
return (w)
}
##牛顿法优化函数（带正则项）
Newt1<-function(x,w,d,n,y,k){
  i<-0
  repeat{
    i<-i+1
    if(qr(w)$rank!=n+1&w[1,1]!=0)
      break
    w<-w-solve(L2(x,w,d,n,k))%*%L1(x,w,d,n,y,k)
    if(i==1000){
      break
    }
  }
  return (w)
}
##梯度下降法
down<-function(x,w,d,n,y,ste){
  i<-0
  s<-0
  repeat{
    i<-i+1
    if(norm(l1(x,w,d,n,y)*ste)<s){
      ste<-ste/2
    }
    else{
      w<-w+l1(x,w,d,n,y)*ste
      s<-norm(l1(x,w,d,n,y)*ste)
    }
    if(i==1000){
      return (w)
    }
  }
}
#分别生成两个维度的高斯数据
yr<-rnorm(100,mean=0,sd=0.5)
##yq<-rnorm(100,mean=2,sd=0.3)
yr1<-rnorm(100,mean=2,sd=0.3)
yk<-rnorm(100,mean=3,sd=0.4)
yk1<-rnorm(100,mean=1,sd=0.2)
#定义样本数与维度数
d<-200
n<-2
#初始化目标优化对象w
w<-cbind(rep(0,n+1))
#生成数据集
x<-matrix(c(rep(1,200),yr,yk,yr1,yk1),nrow=n+1,ncol=d,byrow=TRUE)
y<-cbind(c(rep(0,100),rep(1,100)))
#牛顿法优化
w<-cbind(rep(0,n+1))
W<-Newt(x,w,d,n,y,10^(-3))
plot(x[2,],x[3,],main='1.2.1',
     xlab='x1', ylab='x2',col=ifelse(y==1, 'blue', 'red'))
xx <- seq(min(x), max(x), length=200)
yy<- (W[1,]+W[2,]*xx)/(-W[3,])
lines(xx,yy)
#牛顿法优化（带正则项）
w<-cbind(rep(0,n+1))
W1<-Newt1(x,w,d,n,y,10^(-3))
plot(x[2,],x[3,],main='1.2.2',
     xlab='x1', ylab='x2',col=ifelse(y==1, 'blue', 'red'))
xx1<- seq(min(x), max(x), length=200)
yy1<- (W1[1,]+W1[2,]*xx1)/(-W1[3,])
lines(xx1,yy1)
#梯度下降法
plot(x[2,],x[3,],main='1.2.3',
     xlab='x1', ylab='x2',col=ifelse(y==1, 'blue', 'red'))
xx2<- seq(min(x), max(x), length=200)
W2<-down(x,w,d,n,y,10^(-1))
yy2<-(W2[1,]+W2[2,]*xx2)/(-W2[3,])
lines(xx2, yy2, lwd=2, col="black")

yr1<-rnorm(100,mean=0,sd=0.5)
yr11<-rnorm(100,mean=0,sd=0.5)
yk1<-rnorm(100,mean=2,sd=1)
yk11<-rnorm(100,mean=2,sd=1)
x1<-matrix(c(rep(1,200),yr1,yk1,yr11,yk11),nrow=n+1,ncol=d,byrow=TRUE)
y1<-cbind(c(rep(0,100),rep(1,100)))
ww<-cbind(rep(0,n+1))
#牛顿法优化
W3<-Newt(x1,ww,d,n,y1)
plot(x1[2,],x1[3,],main='1.2.4',
     xlab='x1', ylab='x2',col=ifelse(y1==1, 'blue', 'red'))
xx3 <- seq(min(x1), max(x1), length=200)
yy3<- (W3[1,]+W3[2,]*xx3)/(-W3[3,])
lines(xx3,yy3)
#牛顿法优化（带正则项）
ww<-cbind(rep(0,n+1))
W4<-Newt1(x1,ww,d,n,y1,10^(-3))
plot(x1[2,],x1[3,],main='1.2.5',
     xlab='x1', ylab='x2',col=ifelse(y1==1, 'blue', 'red'))
xx4<- seq(min(x1), max(x1), length=200)
yy4<- (W4[1,]+W4[2,]*xx4)/(-W4[3,])
lines(xx4,yy4)
#梯度下降法
ww<-cbind(rep(0,n+1))
plot(x1[2,],x1[3,],main='1.2.6',
     xlab='x1', ylab='x2',col=ifelse(y1==1, 'blue', 'red'))
xx5<- seq(min(x1), max(x1), length=200)
W5<-down(x1,ww,d,n,y1,10^(-1))
yy5<-(W5[1,]+W5[2,]*xx5)/(-W5[3,])
lines(xx5, yy5, lwd=2, col="black")

d1<-1327
n1<-4
X<-read.table("hy.txt",sep=",")
hy<-X[1:d1,1:n1]
Y<-cbind(X[,n1+1])
XX<-rbind(c(rep(1,d1)),t(hy))
ww<-cbind(rep(0,n1+1))
W5<-Newt(XX,ww,d1,n1,Y,10^(-3))
plot(XX[3,],XX[4,],main='1.2.7',
     xlab='x1', ylab='x2',col=ifelse(Y==1, 'blue', 'red'))
xx5<- seq(min(XX), max(XX), length=200)
yy5<- (W5[1,]+W5[3,]*xx5)/(-W5[4,])
lines(xx5,yy5)