0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@kozakai-ryouta(亮太 小酒井)

共分散構造分析

0
Last updated at Posted at 2019-11-12

# p.98

library(dplyr)

n=9;

alpha=array(0,dim=c(n,n))

alpha[2,1]=0.1;alpha[3,2]=0.2

alpha[4,1]=0.5;alpha[5,1]=0.5

alpha[6,2]=0.8;alpha[7,2]=0.8

alpha[8,3]=0.4;alpha[9,3]=0.3

vec=rpois(n,sample(10:20,1))

val=alpha%*%vec

u=vec-val

sigma_u=diag(rpois(n,sample(10:30,1)))

t=solve(diag(1,n)-alpha)

vec

t%*%u

sigma_t=t%*%sigma_u%*%t(t)

A_star=solve(sqrt(diag(diag(sigma_t))))%*%alpha%*%sqrt(diag(diag(sigma_t)))

B=solve(sqrt(diag(diag(sigma_t))))%*%sqrt(diag(diag(sigma_t)))

t_star=solve(sqrt(diag(diag(sigma_t))))%*%vec

u_star=solve(sqrt(diag(diag(sigma_u))))%*%u

A_star%*%t_star+B%*%u_star

lt_star=solve(diag(1,n)-A_star)

G=cbind(array(0,dim=c(6,3)),diag(1,6))

sigma_v=G%*%t%*%sigma_u%*%t(t)%*%t(G)

v=G%*%t%*%u

sigma_ru=solve(sqrt(diag(diag(sigma_u))))%*%sigma_u%*%solve(sqrt(diag(diag(sigma_u))))

# sigma_r=G%*%solve(diag(1,n)-A_star)%*%B%*%sigma_ru%*%t(B)%*%t(solve(diag(1,n)-A_star))%*%t(G)



# p.113

n=9

sigma_vec=rpois(n,sample(10:20,1))

sigma=diag(sigma_vec)

alpha=rep(1,n)

x=rpois(n,sample(10:20,1))

S=t(t(x))%*%t(x)


cost=function(a,sig){
  
mat=array(0,dim=c(n,n))

for(i in 1:n){
for(j in 1:n){  

mat[i,j]=(S[i,j]-a[i]*a[j]-sig[i,j])^2    
  
}
}

return(sum(mat))

}


eta=0.0001;h=0.01

ite=30000

for(j in 1:ite){
  
vec=alpha

for(i in 1:length(alpha)){
  
sub=vec;sub[i]=sub[i]+h

alpha[i]=alpha[i]-eta*(cost(sub,diag(sigma_vec))-cost(vec,diag(sigma_vec)))/h
  
}
  
vec=sigma_vec

for(i in 1:length(alpha)){
  
sub=vec;sub[i]=sub[i]+h

sigma_vec[i]=sigma_vec[i]-eta*(cost(alpha,diag(sub))-cost(alpha,diag(vec)))/h
  
}

print(cost(alpha,diag(sigma_vec)))
  
}

predict=t(t(alpha))%*%t(alpha)

diag(predict)=diag(predict)+sigma_vec

GFI=1-sum(diag(solve(predict)%*%(S-predict))^2)/sum(diag(solve(predict)%*%S)^2)






# 潜在クラスモデル　p.220 応用

class=5;n=10

p_mat=array(0.1,dim=c(class,n))

w=rep(0.5,class)

u=array(0,dim=c(class,n))

for(j in 1:class){
  
vec=rpois(2,sample(10:20,1))  
  
u[j,]=sample(0:1,n,replace = T,prob=vec/sum(vec))  
  
}


likelihood=function(u_vec,p_vec,w_vec){
  
return(sum((u_vec*log(p_vec)+(1-u_vec)*log(1-p_vec))*c(w_vec/sum(w_vec))))
  
}


ite=5000

h=0.0001;eta=0.0001

for(j in 1:ite){
  
vec=w

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

sub=w;sub[i]=sub[i]+h    
  
dw=(likelihood(u,p_mat,sub)-likelihood(u,p_mat,vec))/h

if(w[i]+eta*dw>0){

w[i]=w[i]+eta*dw

}

}
 
vec=p_mat
 
for(l in 1:nrow(p_mat)){
for(m in 1:ncol(p_mat)){  
  
sub=vec;sub[l,m]=sub[l,m]+h

dp=(likelihood(u,sub,w)-likelihood(u,vec,w))/h
   
if(p_mat[l,m]+eta*dp<(1-h) & p_mat[l,m]+eta*dp>h){

p_mat[l,m]=p_mat[l,m]+eta*dp

}

} 
}  
 

print(likelihood(u,p_mat,w))
  
}



# p.123 トービット変数　応用

n=1000

score=sort(rpois(n,10))

# score=(score-mean(score))/sd(score)

t_b=max(score[1:floor(n*0.05)])

t_a=min(score[(n-floor(n*0.05)):n])

x=score[score>=t_b & score<=t_a]

N_b=score[score<t_b];N_a=score[score>t_a]



log_lik=function(mu,sigma){
  
return((log(pnorm(t_b,mu,sigma))*length(N_b)+sum(log(dnorm(x,mu,sigma)+10^(-10)))+length(N_a)*log(1-pnorm(t_a,mu,sigma))) ) 
  
}



mu=50

sigma=50

h=10^(-3)

ite=50000

log_lik_vec=c()

log_lik_val_pre=-10^(10);log_lik_val=-10^(10)+1

for(j in 1:ite){
  
log_lik_val_pre=log_lik_val  
  
dmu=(log_lik(mu+h,sigma)-log_lik(mu,sigma))/h

dsigma=(log_lik(mu,sigma+h)-log_lik(mu,sigma))/h
  
H=array(0,dim=c(2,2))

# H[2,1]=(log_lik(mu+h,sigma+h)-log_lik(mu+h,sigma)-log_lik(mu,sigma+h)+log_lik(mu,sigma))/(h^2)
  
# H[1,2]=(log_lik(mu+h,sigma+h)-log_lik(mu+h,sigma)-log_lik(mu,sigma+h)+log_lik(mu,sigma))/(h^2)

dh11=(log_lik(mu+h,sigma)-2*log_lik(mu,sigma)+log_lik(mu-h,sigma))/(2*h^2)

dh22=(log_lik(mu,sigma+h)-2*log_lik(mu,sigma)+log_lik(mu,sigma-h))/(2*h^2)

# eigen_values=eigen(H)$values

# eigen_values[eigen_values==0]=10^(-8)

# sol_H=eigen(H)$vectors%*%diag(1/eigen_values)%*%t(eigen(H)$vectors)

vec=c(mu,sigma)+0.0001*c(dmu,dsigma)

# if(vec[2]>0){

mu=vec[1];sigma=vec[2]

log_lik_vec=c(log_lik_vec,log_lik(mu,sigma))

# }

log_lik_val=log_lik(mu,sigma)

# print(log_lik(mu,sigma))

# }else{
  
# print("stop") 
  
# }

}

plot(log_lik_vec)



# p.125 

library(mvtnorm);library(dplyr)

rm(list=ls())

n=50

xy=rmvnorm(n,c(10,30),matrix(c(10,20,20,10),ncol=2))

colnames(xy)=c("x","y")

xy=data.frame(xy)

tbx=xy[order(xy$x),];tbx=tbx[1:floor(n*0.05),];Nbx=nrow(tbx);sam_tbx=tbx;tbx=max(tbx[,1])

tby=xy[order(xy$y),];tby=tby[1:floor(n*0.05),];Nby=nrow(tby);sam_tby=tby;tby=max(tby[,2])

tax=xy[order(xy$x),];tax=tax[(n-floor(n*0.05)):n,];Nax=nrow(tax);sam_tax=tax;tax=min(tax[,1])

tay=xy[order(xy$y),];tay=tay[(n-floor(n*0.05)):n,];Nay=nrow(tay);sam_tay=tay;tay=min(tay[,2])

Nbb=length(xy$x[xy$x<tbx & xy$y<tby])

Nab=length(xy$x[xy$x>tax & xy$y<tby])

Nba=length(xy$x[xy$x<tbx & xy$y>tay])

Naa=length(xy$x[xy$x>tax & xy$y>tay])

X=xy %>% filter(x>tbx & y>tby & x<tax & y<tay)

X=as.matrix(X)



log_lik=function(mu,sigma){
  
del=0.0001  
  
cost=c(log(pmvnorm(c(tbx,tby),mean=mu,sigma=sigma)[1])*Nbb)

vecs=c()

for(j in 1:nrow(sam_tby)){

vecs=c(vecs,(pmvnorm(lower=c(-Inf,sam_tby[j,2]-del),upper=c(tbx,sam_tby[j,2]),mean=mu,sigma=sigma)[1]/del))  
  
}

vecs[vecs==0]=(10^(-10));vecs=log(vecs);cost=c(cost,vecs)

val=pmvnorm(lower=c(-Inf,tay),upper=c(tbx,Inf),mean=mu,sigma=sigma)[1]

val[val==0]=10^(-10);val=Nba*log(val)

cost=c(cost,val)

vecs=c()

for(j in 1:nrow(sam_tby)){
  
vecs=c(vecs,(pmvnorm(lower=c(sam_tby[j,1]-del,-Inf),upper=c(sam_tby[j,1],tby),mean=mu,sigma=sigma)[1]/del))  
  
}

vecs[vecs==0]=(10^(-10));vecs=log(vecs);cost=c(cost,vecs)

vecs=c()

for(j in 1:nrow(X)){
  
vecs=c(vecs,(dmvnorm(c(X[j,]),mean=mu,sigma=sigma)[1]))  
  
}

vecs[vecs==0]=(10^(-10));vecs=log(vecs);cost=c(cost,vecs)

vecs=c()

for(j in 1:nrow(sam_tay)){
  
vecs=c(vecs,(pmvnorm(lower=c(sam_tay[j,1]-del,tax),upper=c(sam_tay[j,1],Inf),mean=mu,sigma=sigma)[1]/del))  
  
}

vecs[vecs==0]=(10^(-10));vecs=log(vecs);cost=c(cost,vecs)

val=pmvnorm(lower=c(tax,-Inf),upper=c(Inf,tby),mean=mu,sigma=sigma)[1]

val[val==0]=10^(-10);val=Nab*log(val)

cost=c(cost,val)

vecs=c()

for(j in 1:nrow(sam_tax)){
  
vecs=c(vecs,(pmvnorm(lower=c(tax,sam_tax[j,2]-del),upper=c(Inf,sam_tax[j,2]),mean=mu,sigma=sigma)[1]/del))  
  
}

vecs[vecs==0]=(10^(-10));vecs=log(vecs);cost=c(cost,vecs)

val=(pmvnorm(lower=c(tax,tay),upper=c(Inf,Inf),mean=mu,sigma=sigma))

val[val==0]=(10^-10)

cost=c(cost,Naa*log(val))

return(as.numeric(sum(cost)))

}



mu_vec=as.numeric(rep(10,2));sigma=array(10,dim=c(2,2))

eta=0.0001;h=0.001

ite=100

for(j in 1:ite){
  
vec=mu_vec

for(i in 1:length(mu_vec)){
  
mu_sub=vec;mu_sub[i]=mu_sub[i]+h  
  
vec[i]=(log_lik(mu_sub,sigma)-log_lik(vec,sigma))/h  
  
}

dmu=vec
  
mat=sigma;dsigma=mat

mat_sub=mat;mat_sub[1,1]=mat_sub[1,1]+h

dsigma[1,1]=(log_lik(mu_vec,mat_sub)-log_lik(mu_vec,mat))/h

mat_sub=mat;mat_sub[2,2]=mat_sub[2,2]+h

dsigma[2,2]=(log_lik(mu_vec,mat_sub)-log_lik(mu_vec,mat))/h

mat_sub=mat;mat_sub[1,2]=mat_sub[1,2]+h;mat_sub[2,1]=mat_sub[2,1]+h

dsigma[1,2]=(log_lik(mu_vec,mat_sub)-log_lik(mu_vec,mat))/h
 
dsigma[2,1]=(log_lik(mu_vec,mat_sub)-log_lik(mu_vec,mat))/h
 

print(log_lik(mu_vec,sigma))
  
}
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?