異常検知と変化検知　講談社(機械学習プロフェショナルシリーズ)  　井出剛・杉山将先生

```r

#混合正規分布の忘却アルゴリズム p.63

library(mvtnorm)
library(dplyr)

data(iris)

X=iris[,1:4]
class=iris[,5]
summary=iris%>%group_by(Species)%>%summarise(Sepal.Length=mean(Sepal.Length),Sepal.Width=mean(Sepal.Width),Petal.Length=mean(Petal.Length),Petal.Width=mean(Petal.Width))



r=1　　 #スムージングパラメータ
beta=0.1　#忘却率
K=3　　　#想定されるクラスタ数

pis=rep(1,K)/K

sample_no=sample(c(1:nrow(X)),3)
mu=X[sample_no,]
sig=list()
for(j in 1:K){sig=c(sig,list(diag(rep(1,ncol(X)))))}

ite=nrow(X)
axvec=c()
for(l in 1:ite){
pispre=pis;sigpre=sig;mupre=mu
vec=as.numeric(X[l,])
qk=c()
for(j in 1:K){
qk=c(qk,pis[j]*dmvnorm(vec,as.numeric(mu[j,]),sig[[j]])) 
}

#異常度
ax=-log(sum(qk))
axvec=c(axvec,ax)

qk=qk/sum(qk)
pis=pis*(1-beta)+beta*qk

error=sum(abs(pispre-pis))
for(j in 1:K){
mu[j,]=mu[j,]*(1-beta)+beta*qk[j]*vec  
sig[[j]]=sig[[j]]*(1-beta)+beta*qk[j]*(t(t(c(vec)))%*%t(c(vec)))  
error=error+sum(abs(mu[j,]-mupre[j,]))+sum(abs(sig[[j]]-sigpre[[j]]))
}
print(error)
}

pis_est=(pis+r)/(K*r+sum(pis))
mu_est=mu/pis
sig_est=sig
for(j in 1:K){
sig_est[[j]]=sig_est[[j]]/pis[j]-t(t(as.numeric(mu[j,])))%*%t(as.numeric(mu[j,]))
}

iris_dat=data.frame(iris,singularity=axvec)


```
```r
#p.79

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),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))

X=data[,2:4]
X=matrix(unlist(X),nrow(X))
X=t((t(X)-apply(X,2,mean))/apply(X,2,sd))

m=apply(X,2,mean)
mu=m/sqrt(sum(m^2))
ax=1-X%*%mu

mmo=2*mean(ax)^2/(mean(ax^2)-mean(ax)^2)
smo=(mean(ax^2)-mean(ax)^2)/(2*mean(ax))
qchisq(0.9,mmo,smo)

```

```r

#p.91

#ガウス過程回帰による異常検知　

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),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))

X=data[,3:4]
X=matrix(unlist(X),nrow(X))
Y=data$y
N=nrow(data)

dist=dist(X,method="euclidean")
metric=array(0,dim=c(N,N))
metric[lower.tri(metric)]=c(dist)
metric=metric+t(metric)

RBF=function(r,l){
 return(r*exp(-metric/(l^2)) )
}

RBF_part=function(r,l,x,y){
return(r*exp(-sum((x-y)^2)/(l^2)))
}


Mlik=function(r,l){
K_childer=RBF(r,l)
sigma=as.numeric(c(Y)%*%solve(diag(1,N)+RBF(r,l))%*%c(Y))/N
#-N*(log(2*pi/N)+1)-log(det(diag(1,N)+K_childer))/2-N*log(as.numeric(c(Y)%*%solve(diag(1,N)+RBF(R,L))%*%c(Y)))/2
return(-N*log(2*pi*sigma)/2-log(det(diag(1,N)+K_childer))/2-N/2) 
}

#初期パラメータ
R=1
L=1

eta=0.01
ite=3*10^4
h=0.001
for(l in 1:ite){
R=R+eta*(Mlik(R+h,L)-Mlik(R,L))/h
L=L+eta*(Mlik(R,L+h)-Mlik(R,L))/h
#print(Mlik(R,L))
}

sig_hat=as.numeric(c(Y)%*%solve(diag(1,N)+RBF(R,L))%*%c(Y))/N


mu_y=function(x){
K_childer=RBF(R,L)
kvec=c()
for(j in 1:nrow(X)){
kvec=c(kvec,RBF_part(R,L,X[j,],x))
}
K0=RBF_part(R,L,x,x)
return(as.numeric(c(kvec)%*%solve(K_childer+sig_hat*diag(1,N))%*%c(Y)))
}
  
sig_y=function(x){
K_childer=RBF(R,L)
kvec=c()
for(j in 1:nrow(X)){
kvec=c(kvec,RBF_part(R,L,X[j,],x))
}
K0=RBF_part(R,L,x,x)
return(as.numeric(sig_hat+K0-c(kvec)%*%solve(K_childer+sig_hat*diag(1,N))%*%c(kvec)))
}

mu_y_vec=c()
sig_y_vec=c()

for(j in 1:N){
mu_y_vec=c(mu_y_vec,mu_y(X[j,]))
sig_y_vec=c(sig_y_vec,sig_y(X[j,]))
}

#プロット
plot(c(1:N),mu_y_vec,xlab="index",ylab="values",xlim=c(1,N),ylim=c(min(c(Y,mu_y_vec)),max(c(Y,mu_y_vec))),col=2)

par(new=T)

plot(c(1:N),Y,xlab="index",ylab="values",xlim=c(1,N),ylim=c(min(c(Y,mu_y_vec)),max(c(Y,mu_y_vec))),col=3,main="mu_y_vec")

#異常度
axy=-log(2*pi*sig_y_vec)/2+((Y-mu_y_vec)^2)/(2*sig_y_vec)

#リッジ回帰のパラメータの場合
alpha=solve(t(X)%*%X+diag(sig_hat,ncol(X)))%*%t(X)%*%Y

#プロット
plot(c(1:N),X%*%alpha,xlab="index",ylab="values",xlim=c(1,N),ylim=c(min(c(Y,X%*%alpha)),max(c(Y,X%*%alpha))),col=2)

par(new=T)

plot(c(1:N),Y,xlab="index",ylab="values",xlim=c(1,N),ylim=c(min(c(Y,X%*%alpha)),max(c(Y,X%*%alpha))),col=3,main="Ridge")


```

```r

#確認中

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),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))

X=data[,2:4]
X=matrix(unlist(X),nrow(X))
C=X%*%t(X)

#ランチョス法

s=6
M=ncol(C)
a=rep(0,M)
beta=1
r=rep(1,M)/sqrt(M)

alphavec=c()
betavec=c()
avec=array(0,dim=c(s,M))
for(j in 1:s){
a_pre=a
a=r/beta
avec[j,]=a
alpha=c(a)%*%C%*%c(a)
alphavec=c(alphavec,alpha)
r=C%*%c(a)-as.numeric(alpha)*c(a)-beta*c(a_pre)
beta=sqrt(sum(r^2))
betavec=c(betavec,beta)
}
A=t(avec)

Ts=array(0,dim=c(s,s))
diag(Ts)=alphavec
for(j in 1:(s-1)){
Ts[j,(j+1)]=betavec[j] 
Ts[(j+1),j]=betavec[j]
}
C%*%A
A%*%Ts



```


```r

#p.134 疎なガウス型グラフィカルモデルの学習

library(mvtnorm)

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),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))

X=iris[,1:4]
X=matrix(unlist(X),nrow(X))
S=cor(X)
M=ncol(S)

#正則化パラメータ
p=0.1

lam=S
inv_lam=S
eta=0.01

ite=100
for(l in 1:ite){
for(i in 1:M){
W=inv_lam[-i,-i]
s=inv_lam[i,-i]
sii=inv_lam[i,i]
W0=W;diag(W0)=0

beta=rep(1,M-1)
beta_pre=beta+1
while(sum(abs(beta-beta_pre))>10^(-8))({
beta_pre=beta  
Aj=s-W%*%beta  
beta=beta*(1-eta)+eta*ifelse(Aj>(p),(Aj-p)/diag(W),ifelse(Aj<(-p),(Aj+p)/diag(W),0)) 
#print(sum(abs(beta-beta_pre)))
})

w=W%*%beta
sig=sii+p  
lambda=as.numeric(1/(sig-c(w)%*%W%*%c(w)))
l=-solve(W)%*%w/as.numeric(sig-c(w)%*%W%*%c(w))
vec=rep(1,M)
vec[-i]=l;vec[i]=lambda
lam[i,]=vec;lam[,i]=vec
vec=rep(1,M)
vec[-i]=w;vec[i]=sig
inv_lam[i,]=vec;inv_lam[,i]=vec
loglik=log(p*exp(-p*sum(abs(lam)))/2)+sum(log(1+dmvnorm(X,mean=rep(0,M),lam)))
print(loglik)
}
}  

```

```r

#p.152

data(iris)
iris=matrix(unlist(iris),nrow=nrow(iris))
iris=t((t(iris)-apply(iris,2,mean))/apply(iris,2,sd))

X_train=iris[1:25,-5]
X_test=iris[26:59,-5]
N=nrow(X_train)

t=5
ht=sample(seq(0,1,0.01),t)

K=4
X_train_no=sample(1:K,nrow(X_train),replace=T)
X_test_no=sample(1:K,nrow(X_test),replace=T)

eta=10^(-1)

pthi=function(x,Xtrain,h){
return(exp(-apply((t(Xtrain)-c(x))^2,2,sum)/(2*h^2)))
}

J=function(sit,Xtrain,Xtest,h){
cost=c() 
for(i in 1:nrow(Xtest)){cost=c(cost,sum(sit*pthi(Xtest[i,],Xtrain,h)))} 
cost=mean(cost)

cost2=c()
for(i in 1:nrow(Xtrain)){cost2=c(cost2,sum(sit*pthi(Xtrain[i,],Xtrain,h)))}
cost2=mean(log(cost2))

return(cost-cost2)
}

Jvec=c()

for(j in 1:t){
Jkvec=c()
for(k in 1:K){
  
X_train_dash=X_train[!(X_train_no==k),]
X_test_dash=X_test[!(X_test_no==k),]
sita=rep(5,nrow(X_train_dash))  
  
gradient=rep(0,nrow(X_train_dash))
for(i in 1:nrow(X_test_dash)){gradient=gradient+pthi(X_test_dash[i,],X_train_dash,ht[j])}
gradient=gradient/nrow(X_test_dash)

sita_pre=sita+1
while(abs(J(sita,X_train_dash,X_test_dash,ht[j])-J(sita_pre,X_train_dash,X_test_dash,ht[j]))>10^(-4))({
sita_pre=sita
gradient2=rep(0,nrow(X_train_dash))
for(i in 1:nrow(X_train_dash)){
gradient2=gradient2+pthi(X_train_dash[i,],X_train_dash,ht[j])/sum(pthi(X_train_dash[i,],X_train_dash,ht[j])*sita)  
}
dJ=gradient-gradient2/(nrow(X_train_dash))
sita=sita-eta*dJ
#print(abs(J(sita,X_train_dash,X_test_dash,ht[j])-J(sita_pre,X_train_dash,X_test_dash,ht[j])))
})

Jkvec=c(Jkvec,J(sita,X_train_dash,X_test_dash,ht[j])) 
}  
Jvec=c(Jvec,mean(Jkvec))  
}



h_hat=ht[c(1:t)[Jvec==min(Jvec)][1]]
sita=rep(5,N)

gradient=rep(0,nrow(X_train))
for(i in 1:nrow(X_test)){gradient=gradient+pthi(X_test[i,],X_train,ht[j])}
gradient=gradient/nrow(X_test)

sita_pre=sita+1
while(abs(J(sita,X_train,X_test,h_hat)-J(sita_pre,X_train,X_test,h_hat))>10^(-4))({
sita_pre=sita
gradient2=rep(0,nrow(X_train))
for(i in 1:nrow(X_train)){
gradient2=gradient2+pthi(X_train[i,],X_train,h_hat)/sum(pthi(X_train[i,],X_train,h_hat)*sita)  
}
dJ=gradient-gradient2/(nrow(X_train))
sita=sita-eta*dJ
print(abs(J(sita,X_train,X_test,h_hat)-J(sita_pre,X_train,X_test,h_hat)))
})

#異常度
axy=c()
for(i in 1:nrow(X_test)){
axy=c(axy,-log(sum(pthi(X_test[i,],X_train,h_hat)*sita)))  
}




```