#p.75
m=10;n=20
a=(array(rnorm(m*n),dim=c(m,n)))
beta=rep(0.01,m)
cost=function(sit,bs){
V2=t(a)%*%bs
P_k2=exp(sit*V2)/sum(exp(sit*V2))
return(sum(log(exp(sit*V2)))/sit)
}
ite=100000
sita=0.01
h=0.001
for(l in 1:ite){
log_lik_pre=cost(sita,beta)
V=t(a)%*%beta
P_k=exp(sita*V)/sum(exp(sita*V))
db=a%*%P_k;db2=((a^2)%*%P_k)-(a%*%P_k)^2
dsita=(cost(sita+h,beta)-cost(sita,beta))/h
dsita2=(cost(sita+h,beta)-2*cost(sita,beta)+cost(sita-h,beta))/(h^2)
H=diag(c(dsita2,db2));G=diag(1/c(db2));dg=c(db)
params=c(beta)+h*G%*%dg
beta=params
print(log(sum(exp(sita*V)))/sita)
}
#p.75 applied
m=10;n=20
a=(array(rnorm(m*n),dim=c(m,n)))
beta=rep(0.01,m)
cost=function(sit,bs){
V2=t(a)%*%bs
P_k2=exp(sit*V2)/sum(exp(sit*V2))
return(sum(log((P_k2)))/sit)
}
ite=100000
sita=1
h=0.001
log_lik_pre=-10;log_lik=-9
for(l in 1:ite){
if(abs(log_lik_pre-log_lik)>10^(-5)){
log_lik_pre=cost(sita,beta)
V=t(a)%*%beta
P_k=exp(sita*V)/sum(exp(sita*V))
ave=a%*%P_k;var=((a^2)%*%P_k)-(a%*%P_k)^2
dsita=(cost(sita+h,beta)-cost(sita,beta))/h
dsita2=(cost(sita+h,beta)-2*cost(sita,beta)+cost(sita-h,beta))/(h^2)
db=c();db2=c()
vec=beta
for(j in 1:length(beta)){
vec_sub1=vec;vec_sub1[j]=vec_sub1[j]+h
vec_sub2=vec;vec_sub2[j]=vec_sub2[j]-h
db=c(db,(cost(sita,vec_sub1)-cost(sita,vec))/h)
db2=c(db2,(cost(sita,vec_sub1)-2*cost(sita,beta)+cost(sita,vec_sub2))/(h^2))
}
H=diag(c(dsita2,db2));G=diag(1/c(dsita2,db2));dg=c(dsita,db)
params=c(sita,beta)-h*G%*%dg
sita=params[1];beta=params[-1]
log_lik=cost(sita,beta)
}
print(cost(sita,beta))
}
#p.79
library(dplyr)
rs=10;k=5
C_k_rs=array(sample(seq(0,1,0.01),k*rs,replace=T),dim=c(k,rs))
V_k_rs=array(sample(seq(0,1,0.01),k*rs,replace=T),dim=c(k,rs))
sita=10;E_hat=5
L=function(v,zf){
exp_mat=exp(-sita*v)
sum_exp_mat=apply(exp_mat,2,sum)
p_k_rs=t(t(exp_mat)/c(sum_exp_mat))
q_rs=rep(1/rs,rs);f_k_rs=t(t(p_k_rs)*q_rs)
eta_rs=log(sum_exp_mat)-log(q_rs)-1
H_rs=-apply(p_k_rs,1,sum)*log(apply(p_k_rs,1,sum))
if(zf==0){
return(sita*(E_hat-sum(f_k_rs*C_k_rs))^2+sum(eta_rs*(q_rs-apply(f_k_rs,2,sum))^2))
}
if(zf==1){
z_f=(sum(q_rs*H_rs)-sum(q_rs*log(q_rs)))^2
return(z_f+sita*(E_hat-sum(f_k_rs*C_k_rs))^2+sum(eta_rs*(q_rs-apply(f_k_rs,2,sum))^2))
}
if(zf==2){
return(sita*(E_hat-sum(f_k_rs*C_k_rs))^2)
}
if(zf==3){
return(sum(f_k_rs*C_k_rs)-sum(f_k_rs*log(f_k_rs))/sita)
}
}
op=3
ite=100000;h=0.01;eta=0.01
for(l in 1:ite){
exp_mat=exp(-sita*C_k_rs)
sum_exp_mat=apply(exp_mat,2,sum)
p_k_rs=t(t(exp_mat)/c(sum_exp_mat))
q_rs=rep(1/rs,rs);f_k_rs=t(t(p_k_rs)*q_rs)
eta_rs=log(sum_exp_mat)-log(q_rs)-1
H_rs=-apply(p_k_rs,1,sum)*log(apply(p_k_rs,1,sum))
vec=C_k_rs
for(i in 1:k){
for(j in 1:rs){
vec_sub=vec;vec_sub[i,j]=vec_sub[i,j]+h
C_k_rs[i,j]=C_k_rs[i,j]-eta*(L(vec_sub,op)-L(vec,op))/h
}}
print(L(C_k_rs,op))
#print(sum(f_k_rs*C_k_rs))
}
#p.79 apply
library(dplyr)
rs=10;k=5
C_k_rs=array(sample(seq(0,1,0.01),k*rs,replace=T),dim=c(k,rs))
V_k_rs=array(sample(c(10:20)*10,k*rs,replace=T),dim=c(k,rs))
sita=1;E_hat=150
L=function(v,zf){
exp_mat=exp(-sita*v)
sum_exp_mat=apply(exp_mat,2,sum)
p_k_rs=t(t(exp_mat)/c(sum_exp_mat))
q_rs=rep(1/rs,rs);f_k_rs=t(t(p_k_rs)*q_rs)
eta_rs=log(sum_exp_mat)-log(q_rs)-1
H_rs=-apply(p_k_rs,1,sum)*log(apply(p_k_rs,1,sum))
if(zf==0){
return(sita*(E_hat-sum(f_k_rs*V_k_rs))^2+sum(eta_rs*(q_rs-apply(f_k_rs,2,sum))^2))
}
if(zf==1){
z_f=(sum(q_rs*H_rs)-sum(q_rs*log(q_rs)))^2
return(z_f+sita*(E_hat-sum(f_k_rs*V_k_rs))^2+sum(eta_rs*(q_rs-apply(f_k_rs,2,sum))^2))
}
if(zf==2){
return(sita*(E_hat-sum(f_k_rs*V_k_rs))^2)
}
if(zf==3){
return(sum(f_k_rs*V_k_rs)+sum(f_k_rs*log(f_k_rs))/sita)
}
}
op=3
ite=100000;h=0.01;eta=0.01
for(l in 1:ite){
exp_mat=exp(-sita*C_k_rs)
sum_exp_mat=apply(exp_mat,2,sum)
p_k_rs=t(t(exp_mat)/c(sum_exp_mat))
q_rs=rep(1/rs,rs);f_k_rs=t(t(p_k_rs)*q_rs)
eta_rs=log(sum_exp_mat)-log(q_rs)-1
H_rs=-apply(p_k_rs,1,sum)*log(apply(p_k_rs,1,sum))
vec=C_k_rs
for(i in 1:k){
for(j in 1:rs){
vec_sub=vec;vec_sub[i,j]=vec_sub[i,j]+h
C_k_rs[i,j]=C_k_rs[i,j]-eta*(L(vec_sub,op)-L(vec,op))/h
}}
print(L(C_k_rs,op))
#print(sum(f_k_rs*log(f_k_rs)))
}
#p.197
library(dplyr)
s=5;r=10;a_s=5
q_rs=array(0.001,dim=c(r,s));v_rs=matrix(sample(1:10,r*s,replace = T),nrow=r)/1000
x_a=rep(0.1,a_s);y_a=sample(1:10,a_s,replace = T)
beta=0.01
ite=100000;eta=0.001
for(k in 1:ite){
cost=function(alpha){
return(abs(sum(pnorm(x_a+alpha*(y_a-x_a)))-sum(pnorm(q_rs+alpha*(v_rs-q_rs)))))
}
for(l in 1:ite){
bs=beta-eta*((cost(beta+h)-cost(beta))/h)
if(bs>=0 & bs<=1){
beta=beta-eta*((cost(beta+h)-cost(beta))/h)
}
}
x_a=x_a+beta*(y_a-x_a)
q_rs=q_rs+beta*(v_rs-q_rs)
print(cost(beta))
if(beta<10^(-10)){
print("stop")
}
}
#p.254
n=10;a=10
x_a=sample(1:10,a,replace=T)
p_a_rs=array(0,dim=c(1,n+1))
for(j in 1:a){
p_a_rs_sub=cbind(rep(j,n),array(sample(seq(0.01,1,0.01),n^2,replace=T),dim=c(n,n)) )
p_a_rs=rbind(p_a_rs,p_a_rs_sub)
}
p_a_rs=tail(p_a_rs,nrow(p_a_rs)-1)
p_a_rs=data.frame(p_a_rs) %>% rename(a=X1)
cost=function(lambda,Q_rs_mat){
Q_rs_mat=abs(Q_rs_mat);lambda=abs(lambda)
t=sum(Q_rs_mat)
dlam=c()
for(j in 1:a){
p_a_rs_sub=p_a_rs %>% filter(a==j) %>% select(-a)
dlam=c(dlam,(sum(Q_rs_mat*p_a_rs_sub)-x_a[j])^2)
}
L=t*log(t)-t-sum(Q_rs_mat*log(Q_rs_mat))+sum(Q_rs_mat*log(Q_rs_mat/sum(Q_rs_mat)))+sum((lambda^2)*dlam)
return(L)
}
ite=100
h=0.001
lam=rep(1,a);Q_rs=array(1,dim=c(n,n))
for(l in 1:ite){
vec=Q_rs
dQ=array(0,dim=c(n,n));ddQ=array(0,dim=c(n,n))
for(i in 1:n){
for(j in 1:n){
vec_sub1=vec;vec_sub1[i,j]=vec_sub1[i,j]+h
vec_sub2=vec;vec_sub2[i,j]=vec_sub2[i,j]-h
dQ[i,j]=(cost(lam,vec_sub1)-cost(lam,vec))/h
ddQ[i,j]=(cost(lam,vec_sub1)-2*cost(lam,vec)+cost(lam,vec_sub2))/(h^2)
}}
Q_rs=Q_rs-h*dQ/ddQ
vec=lam
for(j in 1:a){
vec_sub1=vec;vec_sub1[j]=vec_sub1[j]+h
vec_sub2=vec;vec_sub2[j]=vec_sub2[j]-h
dlam=(cost(vec_sub1,Q_rs)-cost(vec,Q_rs))/h
ddlam=(cost(vec_sub1,Q_rs)-2*cost(vec,Q_rs)+cost(vec_sub2,Q_rs))/(h^2)
lam[j]=vec[j]-h*dlam/ddlam
}
print(cost(lam,Q_rs))
}
#p.256
m=10;n=20;K=10
U=array(sample(seq(0.01,1,0.01),m*n,replace=T),dim=c(m,n))
p=sample(seq(0.01,1,0.01),n)
x_a=sample(1:10,m,replace = T)
cost=function(o,lambda,vs,t_vec){
o=abs(o);lambda=abs(lambda);
return(K*sum((t_vec*p-o)^2)+(vs^2)*(t_vec-sum(o))^2+sum((lambda^2)*(U%*%o-x_a)^2))
}
O_s=rep(1,n)
lam=rep(1,m)
v=1;t=10
ite=100000
h=0.01
cost_pre_val=1;cost_val=0.5
for(l in 1:ite){
if(abs(cost_pre_val-cost_val)>10^(-8)){
cost_pre_val=cost_val
vec=O_s
do_s=c();ddo_s=c()
for(j in 1:length(O_s)){
vec_sub1=vec;vec_sub1[j]=vec_sub1[j]+h
vec_sub2=vec;vec_sub2[j]=vec_sub2[j]-h
do_s=c(do_s,(cost(vec_sub1,lam,v,t)-cost(vec,lam,v,t))/h)
ddo_s=c(ddo_s,(cost(vec_sub1,lam,v,t)-2*cost(vec,lam,v,t)+cost(vec_sub2,lam,v,t))/(h^2))
}
O_s=O_s-h*do_s/ddo_s;O_s=abs(O_s)
vec=lam
dlam=c();ddlam=c()
for(j in 1:length(lam)){
vec_sub1=vec;vec_sub1[j]=vec_sub1[j]+h
vec_sub2=vec;vec_sub2[j]=vec_sub2[j]-h
dlam=c(dlam,(cost(O_s,vec_sub1,v,t)-cost(O_s,vec,v,t))/h)
ddlam=c(ddlam,(cost(O_s,vec_sub1,v,t)-2*cost(O_s,vec,v,t)+cost(O_s,vec_sub2,v,t))/(h^2))
}
lam=lam-h*dlam/ddlam
dv=(cost(O_s,lam,v+h,t)-cost(O_s,lam,v,t))/h
ddv=(cost(O_s,lam,v+h,t)-2*cost(O_s,lam,v,t)+cost(O_s,lam,v-h,t))/(h^2)
v=v-h*dv/ddv
dt=(cost(O_s,lam,v,t+h)-cost(O_s,lam,v,t))/h
ddt=(cost(O_s,lam,v,t+h)-2*cost(O_s,lam,v,t)+cost(O_s,lam,v,t-h))/(h^2)
t=t-h*dt/ddt
cost_val=cost(O_s,lam,v,t)
}
print(cost(O_s,lam,v,t))
}