共分散構造分析

```{r}

#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)

```

```{r}

#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)



```

```{r}


#潜在クラスモデル　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))
  
}

```

```{r}

#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)

```

```{r}

# 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))
  
}



```