library(dplyr)
m=1000
X=sample(20:400,m,replace=T)
C=t(t(X))%*%t(X)+diag(1,m)
y=sample(1:1000,m,replace=T)
V=sample(y,30);E=c(1:1000)
dat_y=data.frame(num=1:length(y),y=y,sign=ifelse(y %in% V,1,0))
b_var=sum(dat_y$y[dat_y$sign==1])/length(V)
y=y-b_var
a=1;b=1;r=1
ite=30000;ite2=300
a_vec=c();b_vec=c();r_vec=c();cost_vec=c()
for(j in 1:ite){
u=as.numeric(sum(diag(C%*%solve((b+r)*diag(1,m)+a*C)))/length(E)+((b^2)*t(y)%*%C)%*%((solve((b+r)*diag(1,m)+a*C))^2)%*%t(t(y))/length(E))
v=as.numeric(sum(diag(solve((b+r)*diag(1,m)+a*C)))/length(E)+((b^2)*t(y))%*%((solve((b+r)*diag(1,m)+a*C))^2)%*%t(t(y)))
b=as.numeric(sum(diag(solve((b+r)*diag(1,m)+a*C)))/length(V)+(b*t(y))%*%((r*diag(1,m)+a*C)^2)%*%((solve((b+r)*diag(1,m)+a*C))^2)%*%t(t(y))/length(V))
a_pre=a;r_pre=r
for(i in 1:ite2){
a_pre=a;r_pre=r
a=a*sqrt(sum(diag(C%*%solve(r*diag(1,m)+a*C)))/(u*length(E)))
r=r*sqrt(sum(diag(solve(r*diag(1,m)+a*C)))/(v*length(E)))
#print(paste0("a:",a,"; r:",r))
}
a_vec=c(a_vec,a);b_vec=c(b_vec,b);r_vec=c(r_vec,r)
#print(paste0("a:",a,"; b:",b,"; r:",r))
cost=sqrt(det(b*(r*diag(1,m)+a*C))/(((2*pi)^length(V))*det((b+r)*diag(1,m)+a*C)))*exp(-b*t(y)%*%(r*diag(1,m)+a*C)%*%solve((b+r)*diag(1,m)+a*C)%*%t(t(y)))
#print(cost)
cost_vec=c(cost_vec,cost)
}
print(paste0("a:",a,"; b:",b,"; r:",r))
x=b*solve((b+r)*diag(1,m)+a*C)%*%t(t(y))
More than 3 years have passed since last update.
Register as a new user and use Qiita more conveniently
- You get articles that match your needs
- You can efficiently read back useful information
- You can use dark theme
List of users who liked
00