入門パターン認識と機械学習　コロナ社　後藤正幸・小林学先生

```{r}



#多クラス分類のサポートベクトルマシン 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))))

}

```


```{r}

#多クラス分類のサポートベクトルマシン 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)


```




```{r}

#関連ベクトルマシン　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))

}  

}



```


```{r}



#混合正規分布モデル 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)

}



```

```{r}



#潜在クラスモデル 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 
matrix_list <- list()
for(j in 1:length(n)){
matrix_list[[j]]=(t(X)-mu[j,])%*%diag(r[,j])%*%t((t(X)-mu[j,]))/n[j]}
sig=matrix_list

for(j in 1:class){
p=dmvnorm(X,mu[j,],sig[[j]])
vec=p*w[j]
r[,j]=vec
}

pvec=r
r=r/apply(r,1,sum)
if(sum(is.nan(r)>0)>0){
num=max(c(1:length(w))[w==max(w)])  
nan_num=c(1:nrow(r))[is.nan(apply(r,1,sum))>0]
rsub=rep(0,length(w));rsub[num]=1
for(j in nan_num){
r[j,]=rsub  
}}
loglik=sum(log(apply(pvec,1,sum)))
print(loglik)
}

km=prod(dim(matrix_list[[1]]))*class+prod(dim(mu))+length(n)+length(w)
AIC=loglik-km


```