ガウス型基底関数（正則化項付きGIC)

library(dplyr)

data=data.frame(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))

data[,colnames(data) %in% c("x1","x2")]=data[,colnames(data) %in% c("x1","x2")]/10

clusters=3;X=as.matrix(data[,colnames(data) %in% c("x1","x2")]);Y=data$y

data=data%>%mutate(clusters=kmeans(data,clusters)$cluster)

params=data%>%group_by(clusters)%>%summarise(x1=mean(x1),x2=mean(x2))

hs=data%>%group_by(clusters)%>%summarise(x1=sum((x1-mean(x1))^2),x2=sum((x2-mean(x2))^2),n=n())

hs=hs%>%mutate(h=(x1+x2)/n)

n=nrow(data)

B=array(0,dim=c(nrow(X),clusters))

for(j in 1:clusters){

vec=params%>%filter(clusters==j)%>%select(-clusters)

vec=as.numeric(vec)

pthi_vec=exp(-apply((t(X)-vec)^2,2,sum)/(2*hs$h[hs$clusters==j]))

B[,j]=pthi_vec

}

B=cbind(rep(1,nrow(B)),B)



GIC=function(ws,sig,lam){

K=diag(1,ncol(B));lam_mat=diag(c(Y-B%*%ws));sig=sig^2

R=cbind(t(B)%*%B+n*lam*sig*K,t(B)%*%lam_mat%*%rep(1,nrow(B)))

R=rbind(R,cbind((t(rep(1,nrow(B)))%*%lam_mat%*%B/sig),n/(2*sig)))

Q11=t(B)%*%(lam_mat^2)%*%B/sig-lam*K%*%ws%*%t(rep(1,nrow(B)))%*%lam_mat%*%B

Q12=t(B)%*%(lam_mat^3)%*%rep(1,nrow(B))/(2*sig^2)-t(B)%*%lam_mat%*%rep(1,nrow(B))/(2*sig)

Q21=t(rep(1,nrow(B)))%*%(lam_mat^3)%*%B/(2*sig^2)-t(rep(1,nrow(B)))%*%lam_mat%*%B/(2*sig^2)

Q22=t(rep(1,nrow(B)))%*%(lam_mat^4)%*%rep(1,nrow(B))/(4*sig^3)-n/(4*sig)

Q=rbind(cbind(Q11,Q12),cbind(Q21,Q22))

value=n*(log(2*pi)+1)+n*log(sig)+2*sum(diag(solve(R)%*%Q))

return(value)

}



w=rep(1,ncol(B))

lambda=1;sigma=1

ite=10^9

h=0.01;eta=10^(-4)

for(l in 1:ite){

vec=w

for(j in 1:length(w)){

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

w[j]=w[j]-eta*(GIC(vec_sub,sigma,lambda)-GIC(vec,sigma,lambda))/h

}

lambda=lambda-eta*(GIC(w,sigma,lambda+h)-GIC(w,sigma,lambda))/h

sigma=sigma-eta*(GIC(w,sigma+h,lambda)-GIC(w,sigma,lambda))/h

#print(GIC(w,sigma,lambda))

}


plot(c(1:length(Y)),Y,xlab="index",ylab="values",xlim=c(1,length(Y)),ylim=c(min(c(Y,B%*%w)),max(c(Y,B%*%w))),col=2,type="o")

par(new=T)

plot(c(1:length(Y)),B%*%w,xlab="index",ylab="values",xlim=c(1,length(Y)),ylim=c(min(c(Y,B%*%w)),max(c(Y,B%*%w))),col=3,type="o")