#多クラス分類のサポートベクトルマシン p.148
#Broyden Quasi algorithm
data(iris)
XX=as.matrix(iris[,!(colnames(iris) %in% "Species")])
Y=as.integer(iris$Species)
kernel=XX%*%t(XX)
sig=ifelse(Y==1,1,-1)
n=nrow(XX)
classes=sort(unique(Y))
C=100
f=function(p){
aa=p[1:n];hh=p[-c(1:n)]
f1=1+hh*sig-(kernel%*%c(aa*sig))*sig
f2=aa*sig
return(c(f1,f2))
}
X=rep(0.01,2*n)
ite=10^4
H=diag(f(X))
eta=10^(-3)
for(l in 1:ite){
X_pre=X
X=X-eta*H%*%f(X)
Xvec=X[1:n]
Xvec[Xvec<0]=0;Xvec[Xvec>C]=C
X[1:n]=Xvec
a=X[1:n];h=X[-c(1:n)]
m=sum(a>0 & a<C)
w0=sum(sig[a>0 & a<C]-kernel[(a>0 & a<C),]%*%c(a*sig))/m
g=w0+t(kernel[a>0,])%*%c(a[a>0]*sig[a>0])
s=X-X_pre
y=f(X)-f(X_pre)
H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H
print(sum(abs(f(X))))
}
#多クラス分類のサポートベクトルマシン p.148
data(iris)
X=as.matrix(iris[,!(colnames(iris) %in% "Species")])
Y=as.integer(iris$Species)
kernel=X%*%t(X)
classes=sort(unique(Y))
C=100
amat=array(0,dim=c(length(classes),length(a)))
hmat=array(0,dim=c(length(classes),length(h)))
w0vec=c();accuracy=c()
ite=10000
gvec=array(0,dim=c(length(Y),length(classes)))
ep=0.01
ep=0.0001
for(i in classes){
a=rep(1,nrow(X))
h=rep(1,nrow(X))
sig=Y;sig[sig!=i]=-1;sig[sig==i]=1
for(l in 1:ite){
apre=a;hpre=h
a=a+ep*(1+h*sig-(kernel%*%c(a*sig))*sig)
s=sum(a>0);m=sum(a>0 & a<C)
w0=sum(sig[a>0 & a<C]-kernel[(a>0 & a<C),]%*%c(a*sig))/m
g=w0+t(kernel[a>0,])%*%c(a[a>0]*sig[a>0])
a[a<0]=0;a[a>C]=C
h=h-ep*a*sig
Ldual=sum(a)-sum(c(kernel%*%c(a*sig))*c(a*sig))
print(sum(abs(apre-a))+sum(abs(hpre-h)))
#print(Ldual)
}
amat[i,]=c(a)
hmat[i,]=c(h)
gvec[,i]=c(g)
w0vec=c(w0vec,w0)
accuracy=c(accuracy,sum(ifelse(g>0,1,-1)*sig)/length(sig))
}
res=ifelse(gvec>=0,1,0)
library(dummies)
sum(ifelse(gvec>=0,1,0)*dummy(Y))/length(Y)
#関連ベクトルマシン p.155
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")])
kernel=array(0,dim=c(nrow(X),nrow(X)))
for(j in 1:nrow(X)){
kernel[j,]=apply((t(X)-c(X[j,]))^2,2,sum)/2
}
G=cbind(rep(1,nrow(X)),log(1+kernel))
beta=2
I=diag(1,nrow(X))
sigk=I/beta
q1=t(G[,1])%*%solve(sigk)%*%y
s1=t(G[,1])%*%solve(sigk)%*%G[,1]
alpha1=ifelse(q1^2-s1>=0,s1^2/(q1^2-s1),Inf)
alpha=c(alpha1,rep(Inf,nrow(X)))
n=nrow(X)
ite=1000
for(l in 1:ite){
no=c(1:length(alpha))[alpha!=Inf]
if(length(no)==1){
siga=1/(diag(c(alpha))[no,no]+(beta*t(G)%*%G)[no,no])
alpha_vec=c()
for(j in 1:(nrow(X)+1)){
Qk=beta*t(G[,j])%*%y-(beta^2)*t(G[,j])%*%G%*%(siga*t(G)%*%y)
Sk=beta*t(G[,j])%*%G[,j]-(beta^2)*t(G[,j])%*%G%*%(siga*t(G)%*%G[,j])
qk=ifelse(alpha[j]==Inf,Qk,alpha[j]*Qk/(alpha[j]-Sk))
sk=ifelse(alpha[j]==Inf,Sk,alpha[j]*Sk/(alpha[j]-Sk))
alphak=ifelse(qk^2-sk>=0,sk^2/(qk^2-sk),Inf)
alpha_vec=c(alpha_vec,alphak)
}
alpha=alpha_vec
mua=beta*siga*t(G)%*%y
beta=(n-sum(1-alpha[no]*diag(siga)))/sum((y-G%*%mua)^2)
}else{
siga=solve(diag(c(alpha))[no,no]+(beta*t(G)%*%G)[no,no])
alpha_vec=c()
for(j in 1:(nrow(X)+1)){
Qk=beta*t(G[,j])%*%y-(beta^2)*((t(G[,j])%*%G)[no])%*%siga%*%t(G[,no])%*%y
Sk=beta*t(G[,j])%*%G[,j]-(beta^2)*((t(G[,j])%*%G)[no])%*%siga%*%t(G[,no])%*%G[,j]
qk=ifelse(alpha[j]==Inf,Qk,alpha[j]*Qk/(alpha[j]-Sk))
sk=ifelse(alpha[j]==Inf,Sk,alpha[j]*Sk/(alpha[j]-Sk))
alphak=ifelse(qk^2-sk>=0,sk^2/(qk^2-sk),Inf)
alpha_vec=c(alpha_vec,alphak)
}
mua=beta*siga%*%t(G[,no])%*%y
beta=(n-sum(1-alpha[no]*diag(siga)))/sum((y-G[,no]%*%mua)^2)
print(sum((y-G[,no]%*%mua)^2))
}
}
#混合正規分布モデル p.196
library(dummies)
library(mvtnorm)
data(iris)
X=as.matrix(iris[,!(colnames(iris) %in% "Species")])
n=apply(dummy(iris$Species),2,sum)
w=n/sum(n)
ite=100
r=matrix(sample(1:9,nrow(X)*length(n),replace=T),nrow=nrow(X))
r=r/apply(r,1,sum)
for(l in 1:ite){
mu=t(t(X)%*%r)/n
sig=array(0,dim=c(length(n),ncol(X)))
for(j in 1:length(n)){
sig[j,]=((t(X)-mu[j,])^2)%*%c(r[,j])/n[j]
}
for(j in 1:length(n)){
p=dmvnorm(X,mu[j,],diag(sig[j,]))
vec=p*w[j]
r[,j]=vec
}
pvec=r
r=r/apply(r,1,sum)
}
#潜在クラスモデル p.201
library(dummies)
library(mvtnorm)
data(iris)
X=as.matrix(iris[,!(colnames(iris) %in% "Species")])
class=3
ite=100
r=matrix(sample(1:9,nrow(X)*class,replace=T),nrow=nrow(X))
r=r/apply(r,1,sum)
for(l in 1:ite){
n=apply(r,2,sum)
w=n/sum(n)
mu=t(t(X)%*%r)/n
sig=array(0,dim=c(class,ncol(X)))
for(j in 1:length(n)){
sig[j,]=((t(X)-mu[j,])^2)%*%c(r[,j])/n[j]
}
for(j in 1:class){
p=dmvnorm(X,mu[j,],diag(sig[j,]))
vec=p*w[j]
r[,j]=vec
}
pvec=r
r=r/apply(r,1,sum)
}