サポートベクターマシン オーム社 小野田崇先生

#回帰問題のためのサポートベクトルマシン p.70


#サンプルデータ
data("airquality")
data=na.omit(airquality)
y=data$Ozone
x=as.matrix(data[,colnames(data)!="Ozone"])
n=nrow(x)

#サンプルデータ
data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)
y=data$y*100
x=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])
n=nrow(x)

#サンプルデータ
data("AirPassengers")
Y=as.numeric(t(AirPassengers))
p=6
y=c()
x=array(0,dim=c(length(Y)-p,p))
n=nrow(x)
for(j in 1:nrow(x)){
val=Y[length(Y)-j+1]  
y=c(y,val)
vec=Y[(length(Y)-j):(length(Y)-j-p+1)]
x[j,]=vec  
}

#Newton Algorithm 

ep=0.001
C=10^(-3)
X=sample(1:(2*nrow(x)+1),(2*nrow(x)+1))
eta=10^(-2)
ite=2*10^3
lambda=1

for(l in 1:ite){
X_pre=X  
lam=X[1]
a=X[-c(1)]
a1=a[c(1:nrow(x))]
a2=a[-c(1:nrow(x))]
f1=y-ep
f2=-y-ep
for(i in 1:n){
value=sum(x%*%c(x[i,])*(a1-a2))+(a1-a2)[i]*n/C
f1[i]=f1[i]+value
f2[i]=f2[i]-value
}
df=c(sum(a1-a2),f1,f2)
mat=x%*%t(x);diag(mat)=diag(mat)+n/C
ddf=rbind(cbind(-mat,mat),cbind(mat,-mat))
ddf=cbind(c(rep(1,length(a1)),rep(-1,length(a2))),ddf)
ddf=rbind(c(0,rep(1,length(a1)),rep(-1,length(a2))),ddf)
X=X+eta*solve(ddf+diag(lambda,nrow(ddf)))%*%df
X=abs(X)
L=sum((a1-a2)*y)-ep*sum(a1+a2)-sum(t(mat*c(a1-a2))*c(a1-a2))/2
f=apply(x%*%t(x)*c(a1-a2),2,sum)
f=f[length(f):1]
b=mean(y-f)+ep
f=f+b
diff=sum((y-f)^2)
print(L)
}




plot(c(1:length(y)),y,xlim=c(1,length(y)),ylim=c(min(c(y,f)),max(c(y,f))),xlab="index",ylab="predict",col=1)

par(new=T)

plot(c(1:length(y)),f,xlim=c(1,length(y)),ylim=c(min(c(y,f)),max(c(y,f))),xlab="index",ylab="predict",col=2)

#教師なし学習へのサポートベクターアプローチ p.125
#外れ値を除外した学習データの生成法


#サンプルデータ
data("airquality")
data=na.omit(airquality)
y=data$Ozone
x=as.matrix(data[,colnames(data)!="Ozone"])
n=nrow(x)

#サンプルデータ
data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)
y=data$y*100
x=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])
n=nrow(x)

#サンプルデータ
data("AirPassengers")
Y=as.numeric(t(AirPassengers))[length(Y):1]
p=6
y=c()
x=array(0,dim=c(length(Y)-p,p))
n=nrow(x)
for(j in 1:nrow(x)){
val=Y[length(Y)-j+1]  
y=c(y,val)
vec=Y[(length(Y)-j):(length(Y)-j-p+1)]
x[j,]=vec  
}

#Newton Algorithm 

ep=0.001
C=10^(-3)
X=rep(1/nrow(x),nrow(x))
eta=10^(-7)
ite=2*10^2
lambda=1

for(l in 1:ite){
X_pre=X  
dX=(diag(x%*%t(x)))-apply((x*c(X))%*%t(x),2,sum)
X=X+eta*dX
X=abs(X)/sum(abs(X))
L=sum(diag(x%*%t(x))*X)-sum((x*c(X))%*%t(x*c(X)))
print(L)
r=diag(x%*%t(x))-2*apply((x*c(X))%*%t(x),2,sum)+sum((x*c(X))%*%t(x*c(X)))

}

#外れ値を除いたデータ
z=x[r<10^5,]

#RBFネットワーク　p.151
#最急降下法、モーメント法
#計算時間がかかるため改良したい

data(iris)
X=iris[,1:4]
species=unique(iris$Species)
species_vector=iris$Species

Y=c()
averages=sample(1:100,length(species))
sds=sample(1:100,length(species))
for(j in 1:length(unique(species))){
Y=c(Y,abs(rnorm(sum(species_vector==species[j]),averages[j],sds[j])))  
}
p=length(Y)

#パラメータ
n=3
clusters=kmeans(X,n)$cluster
Cj=array(0,dim=c(n,ncol(X)))
rj=c()
w=rep(1,n)
lam=10^(0)
for(j in 1:n){
Cj[j,]=apply(X[clusters==j,],2,mean)
Xsub=X[clusters==j,]
rj=c(rj,sqrt(max(apply((t(Xsub)-apply(X[clusters==j,],2,mean))^2,2,sum))))
}

b=0.95
m1=0
m2=0

ite=10^4
eta1=10^(-5)
eta2=10^(-4)
eta3=10^(-3)

for(l in 1:ite){
H=c()
for(j in 1:n){
H=c(H,exp(-apply((t(X)-Cj[j,])^2,2,sum)/(rj[j]^2)))  
}
H=matrix(H,nrow=nrow(X))  
f=H%*%w

drj=c()
dCj=array(0,dim=dim(Cj))
for(j in 1:n){
drj=c(drj,-2*(w[j]/(rj[j]^3))*sum(apply((t(X)-Cj[j,])^2,2,sum)*exp(-apply((t(X)-Cj[j,])^2,2,sum)/(rj[j]^2))*(Y-f)))
dCj[j,]=-2*(w[j]/(rj[j]^2))*apply(t(t(X)-Cj[j,])*c(exp(-apply((t(X)-Cj[j,])^2,2,sum)/(rj[j]^2))*(Y-f)),2,sum)
}
m1=b*m1-eta1*dCj
m2=b*m2-eta2*drj
Cj=Cj+m1
drj=drj+m2


w=w*(1-eta3)+eta3*solve(t(H)%*%H+diag(lam,n))%*%t(H)%*%Y
cost=sum((Y-f)^2)
print(cost)
}