工学系のための最適設計法ー機械学習を活用した理論と実践ー 共立出版 北山哲士先生

#p.52 逐次二次計画法
#過去のレートのデータからM=15000をC=10000からできるだけ近づけるようにする
#パラメータマイナス(=sell)でパラメータプラス(=buy)

data("EuStockMarkets")
mat=as.matrix(EuStockMarkets)/1000


initial_rate=EuStockMarkets[1,]/1000
C=10000
M=15000

g=function(x){
return((sum(x*initial_rate)-C)^2)  
}

f=function(x){
return(mean((M-mat%*%x)^2))  
}

L=function(x,lam){
return(f(x)+lam*g(x))  
}

h=0.001
Mc=10^2
tau=1-10^(-5)
alpha=0.001

param=rep(0.1,4)
z=rep(0.1,4)
LAM=10
ite=3*10^4

B=diag(1,length(param))

for(l in 1:ite){
  
if(l>1){
dL_pre=dL  
}

dL=c()
dg=c()
df=c()
for(j in 1:length(param)){
param_sub=param
param_sub[j]=param_sub[j]+h
dL=c(dL,(L(param_sub,LAM)-L(param,LAM))/h)
dg=c(dg,(g(param_sub)-g(param))/h)
df=c(df,(f(param_sub)-f(param))/h)
}

if(l==1){
q=dL  
s=param
}else{
q=dL-dL_pre
s=param-param_pre
}
q=t(t(c(q)))
s=c(s)

pthi=ifelse(as.numeric(s%*%q)>=0.2*as.numeric(s%*%B%*%s),1,0.8*as.numeric(s%*%B%*%s)/as.numeric(s%*%(B%*%s-q)))

q=q*pthi+(1-pthi)*B%*%s

B=B-(B%*%s)%*%t(B%*%s)/as.numeric(s%*%B%*%s)+q%*%t(q)/as.numeric(s%*%q)

J1=cbind(B,dg)
J2=c(c(dg),0)
J=rbind(J1,J2)
r=c(df,g(param))

U=svd(J)$u
sig=svd(J)$d
V=t(svd(J)$v)
inv_J=t(V)%*%diag(c(1/sig))%*%t(U)

w=inv_J%*%(-r)
param_pre=param
parameters=c(param,LAM)+alpha*w
param=parameters[1:4];parameters=parameters[-c(1:4)]
LAM=parameters[1]
print(mean(mat%*%param))

}

#RBF p.128
data("airquality")

dat=airquality
dat=na.omit(dat)

Y=dat$Ozone
X=cbind(dat[,colnames(dat) %in% c("Temp","Solar.R","Wind")])
X=matrix(as.numeric(unlist(X)),nrow=nrow(X))
Kernel=X%*%t(X) #線形カーネル

lam=0.5
param=solve(t(Kernel)%*%Kernel+diag(rep(lam,ncol(Kernel))))%*%t(Kernel)%*%Y
pred=Kernel%*%param


plot(c(1:length(Y)),pred,xlim=c(1,length(Y)),ylim=c(min(c(Y,pred)),max(c(Y,pred))),col=2,type="p",ylab="",xlab="")

par(new=T)

plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(c(Y,pred)),max(c(Y,pred))),col=3,type="p",ylab="",xlab="")

#カーネルサポートベクトル回帰 p.134
data("airquality")

dat=airquality
dat=na.omit(dat)

Y=dat$Ozone
X=cbind(dat[,colnames(dat) %in% c("Temp","Solar.R","Wind")])
X=matrix(as.numeric(unlist(X)),nrow=nrow(X))
Kernel=X%*%t(X) #線形カーネル
omega=(t(t(c(Y)))%*%t(c(Y)))*Kernel

C=10
mat=omega+diag(1/C,ncol(omega))
mat=rbind(c(0,rep(1,ncol(omega))),cbind(1,mat))

svd_mat=svd(mat)
svd_u=svd_mat$u
svd_d=svd_mat$d
svd_v=svd_mat$v

#もとの行列と同じ
svd_u%*%diag(c(svd_d))%*%t(svd_v)

#一般化逆行列の計算
B=(svd_v)%*%diag(c(1/svd_d))%*%(t(svd_u))
param=B%*%c(0,Y)

b=param[1]
beta=param[-1]
mean(beta)

pred=b+(omega+diag(1/C,ncol(omega)))%*%beta

plot(c(1:length(Y)),pred,xlim=c(1,length(Y)),ylim=c(min(c(Y,pred)),max(c(Y,pred))),col=2,type="p",ylab="",xlab="")

par(new=T)

plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(c(Y,pred)),max(c(Y,pred))),col=3,type="p",ylab="",xlab="")

#カーネルサポートベクトル判別 p.134
data("airquality")

dat=airquality
dat=na.omit(dat)

Y=dat$Ozone
Y=ifelse(Y-mean(Y)>0,1,-1)
X=cbind(dat[,colnames(dat) %in% c("Temp","Solar.R","Wind")])
X=matrix(as.numeric(unlist(X)),nrow=nrow(X))
Kernel=X%*%t(X)　#線形カーネル
omega=(t(t(c(Y)))%*%t(c(Y)))*Kernel

C=10
mat=omega+diag(1/C,ncol(omega))
mat=rbind(c(0,Y),cbind(Y,mat))

svd_mat=svd(mat)
svd_u=svd_mat$u
svd_d=svd_mat$d
svd_v=svd_mat$v

#もとの行列と同じ
svd_u%*%diag(c(svd_d))%*%t(svd_v)

#一般化逆行列の計算
B=(svd_v)%*%diag(c(1/svd_d))%*%(t(svd_u))
param=B%*%c(0,rep(1,length(Y)))

b=param[1]
beta=param[-1]


pred=b+(omega+diag(1/C,ncol(omega)))%*%beta