分類機：SVM（メディア情報処理の基礎数理　尾関和彦先生）

```{r}

#分類機

library(dplyr)

N=100

x=sample(c(1:10),N,replace=TRUE)

y=ifelse(x>5,1,-1)

w=0

eta=0.1

b=0

R=max(abs(x))

for(i in 1:N){

    if((y[i]*(w*x[i]+b)<=0)){
  
      w=w+eta*y[i]*x[i]
  
      b=b+eta*y[i]*R^2

      }

  }

sum(sign(w*x+b)==y)

```

```{r}
#dual perceptron

b=0

m=10000

x=sample(c(1:20),m,replace=TRUE)

y=ifelse(x>5+rnorm(1,mean=0,sd=4),1,-1)

alpha=rep(0,m)

R=max(abs(x))

for(i in 1:m){

    if(y[i]*(sum(alpha*y*x*x[i])+b)<=0){
  
        alpha[i]=alpha[i]+1
    
        b=b+y[i]*R^2
  }
  
}

w=sum(alpha*y*x)

sum(sign(w*x+b)==y)

data_x_y=data.frame(x=x,y=y)

```


```{r}

library(dplyr)

N=1000

x1=sample(c(1:10),N,replace=TRUE)

x2=sample(c(1:10),N,replace = TRUE)

y=ifelse(3*x1-2*x2+5>0,1,-1)

w=rep(0,2)

eta=0.1

b=0

R=max(abs(c(x1,x2)))

for(i in 1:N){

    if((y[i]*(sum(w*c(x1[i],x2[i]))+b)<=0)){
  
      w=w+eta*y[i]*c(x1[i],x2[i])
  
      b=b+eta*y[i]*R^2

      }

  }

sum(sign(cbind(x1,x2)%*%w+b)==y)

data=data.frame(x1=x1,x2=x2,y=y,predict=cbind(x1,x2)%*%w+b)

data=data %>% mutate(sign=ifelse(predict>0,1,-1))



plot(data$x1[data$y==1],data$x2[data$y==1],xlab="x1",ylab="x2",col=2,xlim=c(min(x1),max(x1)),ylim=c(min(x2),max(x2)),main="samples")

par(new=T)

plot(data$x1[data$y==-1],data$x2[data$y==-1],xlab="x1",ylab="x2",col=3,xlim=c(min(x1),max(x1)),ylim=c(min(x2),max(x2)),main="samples")





plot(data$x1[data$sign==1],data$x2[data$sign==1],xlab="x1",ylab="x2",col=2,xlim=c(min(x1),max(x1)),ylim=c(min(x2),max(x2)),main="support vector machine")

par(new=T)

plot(data$x1[data$sign==-1],data$x2[data$sign==-1],xlab="x1",ylab="x2",col=3,xlim=c(min(x1),max(x1)),ylim=c(min(x2),max(x2)),main="support vector machine")


```

```{r}

#perceptron

 data(iris)
 pairs(iris[,1:4], col=iris$Species)

 euclidean.norm = function(x) {sqrt(sum(x * x))}

distance.from.plane = function(z,w,b) {
sum(z*w) + b
}

 classify.linear = function(x,w,b) {
 distances =
 apply(x, 1, distance.from.plane, w, b)
 return(ifelse(distances < 0, -1, +1))
 }

 w = vector(length = ncol(x)) # initialize w
 b = 0 # Initialize b
 k = 0 # count updates
 R = max(apply(x, 1, euclidean.norm))
 ite=100
 
 x=iris[,1:4];y=as.integer(iris[,5])
 
 learning.rate=0.001
 
 for(j in 1:ite){
 
  w_pre=w 
   
 yc <- classify.linear(x,w,b)
 for (i in 1:nrow(x)) {
 if (y[i] != yc[i]) {
 w <- w + learning.rate * y[i]*x[i,]
 b <- b + learning.rate * y[i]*R^2
 k <- k+1
 
  }
 } 
 s = euclidean.norm(w)
 w=w/s;b=b/s;updates=k
 
 print(sum(abs(w-w_pre)))
 
 }
 
 res=data.frame(predict=as.matrix(x)%*%as.numeric(w),y=y)
 
 






```


```{r}

#perceptron2

 data(iris)
 pairs(iris[,1:4], col=iris$Species)

 euclidean.norm = function(x) {sqrt(sum(x * x))}

distance.from.plane = function(z,w,b) {
sum(z*w) + b
}

 classify.linear = function(x,w,b) {
 distances =
 apply(x, 1, distance.from.plane, w, b)
 return(ifelse(distances < 0, -1, +1))
 }

 w = vector(length = ncol(x)) # initialize w
 b = 0 # Initialize b
 k = 0 # count updates
 R = max(apply(x, 1, euclidean.norm))
 
 ite=10000
 
 x=iris[,1:4];y=as.integer(iris[,5])
 
 learning.rate=1
 
 w_box=array(0.1,dim=c(length(unique(y)),ncol(x)))
 
 b_box=rep(0.1,length(unique(y)))
 
for(j in 1:ite){
   
for(l in 1:length(unique(y))){      
   
w=w_box[l,];b=b_box[l] 

yc <- classify.linear(x,w,b)

y_sub=y;y_sub[y_sub!=l]=-1;y_sub[y_sub!=-1]=1

 for (i in 1:nrow(x)) {
   
 if (y_sub[i] != yc[i]) {
   
 w <- w + learning.rate * y_sub[i]*x[i,]
 b <- b + learning.rate * y_sub[i]*R^2
 k <- k+1
 
 }
   
 } 
 s = euclidean.norm(w)
 
 w=w/s;b=b/s;updates=k

 w_box[l,]=as.numeric(w);b_box[l]=as.numeric(b)
 
}

}   
 
 
 


```


```{r}

#gradient algorithm

n=100;m=10

x=array(0,dim=c(n,m));y=c()

for(j in 1:n){
  
x[j,]=rpois(m,sample(1:100,1))

y=c(y,rpois(1,sample(1:100,1)))
  
}




eta=0.0001

K=(as.matrix(x)%*%t(as.matrix(x)))

kernel=function(w){

return(abs(w)/max(abs(K)))

}

alpha=rep(0.01,nrow(x))

ite=100000

C=100

for(j in 1:ite){
  
alpha_pre=alpha  
  
for(l in 1:nrow(x)){
  
value=sum(alpha*y*kernel(as.matrix(x)%*%as.numeric(x[l,])))

alpha[l]=alpha[l]+eta*(1-y[l]*value)  
  
}  
  
W_a=sum(alpha)-sum(kernel(as.matrix(x)%*%t(as.matrix(x)))*(t(t(c(alpha)))%*%t(c(alpha)))*(t(t(c(y)))%*%t(c(y))))/2

#alpha[alpha>C]=C;alpha[alpha<0]=0

print(W_a)  
  
}

 
```

```{r}

#kernel ridge regression

n=50;m=10

x=array(0,dim=c(n,m));y=c()

for(j in 1:n){

x[j,]=rpois(m,sample(1:100,1))

y=c(y,rpois(1,sample(1:100,1)))

}

y=y/100;x=x/100

train_no=c(1:floor(n*0.7))

x.train=x[train_no,];x.test=x[-train_no,]

y.train=y[train_no];y.test=y[-train_no]


K=(as.matrix(x.train)%*%t(as.matrix(x.train)))

kernel=function(w){

return(abs(w)/max(abs(K)))

}





mesh=0.1

sigma_data=data.frame(sigma=seq(mesh,10^4,mesh),GCV=0)

for(j in 1:nrow(sigma_data)){
  
sigma=sigma_data$sigma[j]  
  
alpha=(2*sigma)*solve(kernel(K)+diag(sigma,nrow(x.train)))%*%y.train

w=apply(x.train*c(alpha),2,sum)/(2*sigma)  
  
sigma_data$GCV[j]=mean((x.train%*%w-y.train)^2)  
  
}

plot(sigma_data$sigma,sigma_data$GCV)

sigma=min(sigma_data$sigma[sigma_data$GCV==min(sigma_data$GCV)])


alpha=(2*sigma)*solve(kernel(K)+diag(sigma,nrow(x.train)))%*%y.train

w=apply(x.train*c(alpha),2,sum)/(2*sigma)

K=(as.matrix(x.test)%*%t(as.matrix(x.test)))

kernel=function(w){

return(abs(w)/max(abs(K)))

}


f=c();V=c();

for(j in 1:nrow(x.test)){

#f=c(f,t(y)%*%solve(kernel(K)+diag(sigma,nrow(x)))%*%kernel(as.matrix(x)%*%c(x[j,])))
  
f=c(f,sum(w*x.test[j,]))  

V=c(V,kernel(sum(x.test[j,]^2))-t(kernel(x.test%*%c(x.test[j,])))%*%solve(kernel(K)+diag(sigma,nrow(x.test)))%*%kernel(x.test%*%c(x.test[j,])))


}



plot(c(1:nrow(x.test)),y.test,xlab="index",ylab="values",xlim=c(1,nrow(x.test)),ylim=c(min(y.test,f),max(y.test,f)),col=2,type="o",main="gauss process reg")

par(new=T)

plot(c(1:nrow(x.test)),f,xlab="index",ylab="values",xlim=c(1,nrow(x.test)),ylim=c(min(y.test,f),max(y.test,f)),col=3,type="o",main="gauss process reg")

plot(c(1:nrow(x.test)),f+qnorm(1-0.05/2)*V,xlab="index",ylab="values",xlim=c(1,nrow(x.test)),ylim=c(min(y,f+qnorm(0.05/2)*V),max(y,f+qnorm(1-0.05/2)*V)),col=2,type="o")

par(new=T)

plot(c(1:nrow(x.test)),f+qnorm(0.05/2)*V,xlab="index",ylab="values",xlim=c(1,nrow(x.test)),ylim=c(min(y.test,f+qnorm(0.05/2)*V),max(y.test,f+qnorm(1-0.05/2)*V)),col=2,type="o")

par(new=T)

plot(c(1:nrow(x.test)),f,xlab="index",ylab="values",xlim=c(1,nrow(x.test)),ylim=c(min(y.test,f+qnorm(0.05/2)*V),max(y.test,f+qnorm(1-0.05/2)*V)),col=3,type="o")


EV=sum((y.test-f)^2);RV=sum((f-mean(f))^2)

R2=RV/(EV+RV);R=sqrt(R2)



```

```{r}

#gradient algorithm p.143



data=data.frame(iris)

#x=array(0,dim=c(n,m));y=c()

#for(j in 1:n){
  
#x[j,]=rpois(m,sample(1:100,1))

#y=c(y,rpois(1,sample(1:100,1)))
  
#}

x=as.matrix(data[,!(colnames(data) %in% "Species")])

y=as.integer(data$Species)

n=dim(x)[1];m=dim(x)[2]


eta=0.0001

K=(as.matrix(x)%*%t(as.matrix(x)))

kernel=function(w){

return(abs(w)/max(abs(K)))

}

alpha=rep(0.01,nrow(x))

ite=100000

lagrangean=1000;lagrangean_pre=1001

C=100;C_coef=max(diag(kernel(K)))^2

for(j in 1:ite){
  
if((lagrangean_pre-lagrangean)>0){  
  
alpha_pre=alpha;f=c() 

lagrangean_pre=lagrangean
  
for(l in 1:nrow(x)){
  
value=sum(alpha*y*kernel(as.matrix(x)%*%as.numeric(x[l,])))

f=c(f,value)

alpha[l]=alpha[l]+eta*(1-y[l]*value-(y[l]^2)*alpha[l]/C_coef)  
  
}  

#alpha[alpha<0]=0
  
W_a=sum(alpha)-(sum(kernel(as.matrix(x)%*%t(as.matrix(x)))*(t(t(c(alpha)))%*%t(c(alpha)))*(t(t(c(y)))%*%t(c(y))))+n/C_coef)/2

#alpha[alpha>C]=C;alpha[alpha<0]=0

w=apply(x*c(y*alpha),2,sum);guzai=alpha/C_coef

#print(W_a) 

#print(sum(w^2)+C_coef*sum(guzai^2))

b=(sum(alpha)-sum(alpha*guzai)-sum(alpha*y*(t(w)%*%t(x))))/sum(alpha*y)
 
lagrangean=sum(w^2)/2+sum(guzai^2)*C_coef/2-sum(alpha*(y*t(w)%*%t(x)+y*b-1+guzai))

print(lagrangean)

}else{
  
print("stop")  
  
}
 
}

res=data.frame(guzai=guzai,predict=c(t(w)%*%t(x)+b),y=c(y)) 


```
```{r}

#support vector reg p.155



data=data.frame(iris)

#x=array(0,dim=c(n,m));y=c()

#for(j in 1:n){

#x[j,]=rpois(m,sample(1:100,1))

#y=c(y,rpois(1,sample(1:100,1)))

#}

x=as.matrix(data[,!(colnames(data) %in% "Species")])

y=as.integer(data$Species)

n=dim(x)[1];m=dim(x)[2]


eta=0.0001

K=(as.matrix(x)%*%t(as.matrix(x)))

kernel=function(w){

return(abs(w)/max(abs(K)))

}

alpha=rep(0.01,nrow(x))

ite=100000;epsilon=0.1



C=100;C_coef=max(diag(kernel(K)))^2

for(j in 1:ite){

f=c() 

for(l in 1:nrow(x)){

value=sum(alpha*kernel(as.matrix(x)%*%as.numeric(x[l,])))

f=c(f,value);

dw=(y[l]-epsilon-sum(alpha*kernel(as.matrix(x)%*%as.numeric(x[l,])))-alpha[l]/C_coef)  

alpha[l]=alpha[l]+eta*dw

}  

#alpha[alpha<0]=0

W_a=sum(y*alpha)-epsilon*sum(abs(alpha))-sum(t(t(alpha))%*%t(alpha)*kernel(K))/2-sum(diag(t(t(alpha))%*%t(alpha)))/(2*C_coef)

#print(W_a)

#b=mean(y-f)

b=mean(y-f)-epsilon-mean(alpha)/C_coef

w=apply(x*c(y*alpha),2,sum)

loss=sum((y-f-b)^2)

print(loss)

}

#C_coef*sum(w^2)+


plot(c(1:n),c(y),xlab="index",ylab="values",type="o",col=2,xlim=c(1,n),ylim=c(min(c(y,f)),max(c(y,f))))

par(new=T)

plot(c(1:n),c(f),xlab="index",ylab="values",type="o",col=3,xlim=c(1,n),ylim=c(min(c(y,f)),max(c(y,f))))


```

```{r}


#support vector reg 2 p.155



data=data.frame(iris)

#x=array(0,dim=c(n,m));y=c()

#for(j in 1:n){

#x[j,]=rpois(m,sample(1:100,1))

#y=c(y,rpois(1,sample(1:100,1)))

#}

x=as.matrix(data[,!(colnames(data) %in% "Species")])

y=as.integer(data$Species)

n=dim(x)[1];m=dim(x)[2]


eta=0.001

K=(as.matrix(x)%*%t(as.matrix(x)))

kernel=function(w){

return(abs(w)/max(abs(K)))

}


alpha=rep(0.01,nrow(x))

ite=1000000;epsilon=0.01

h=0.01

C=100;C_coef=max(diag(kernel(K)))^2

for(j in 1:ite){

f=c() 

cost=function(a){

a=a-mean(a)  
  
W_a=sum(y*a)-epsilon*sum(abs(a))-sum(t(t(a))%*%t(a)*kernel(K))/2-sum(diag(t(t(a))%*%t(a)))/(2*C_coef)

return(W_a)

}

vec=alpha

for(l in 1:nrow(x)){

value=sum(alpha*kernel(as.matrix(x)%*%as.numeric(x[l,])))

f=c(f,value);

vec_sub=vec;vec_sub[l]=vec_sub[l]+h

dw=(cost(vec_sub)-cost(vec))/h

alpha[l]=alpha[l]+eta*dw

}  

#alpha[alpha<0]=0

alpha=alpha-mean(alpha)

W_a=cost(alpha)

#print(W_a)

#b=mean(y-f)

b=mean(y-f)-mean(alpha)/C_coef

w=apply(x*c(y*alpha),2,sum)

loss=sum((y-f-b)^2)

print(loss)

}




plot(c(1:n),c(y),xlab="index",ylab="values",type="o",col=2,xlim=c(1,n),ylim=c(min(c(y,f)),max(c(y,f))))

par(new=T)

plot(c(1:n),c(f),xlab="index",ylab="values",type="o",col=3,xlim=c(1,n),ylim=c(min(c(y,f)),max(c(y,f))))




```