潜在構造分析（多変量解析入門Ⅱ 河口至商先生）

とりあえずやってみました。結果はあっていますが、修正が必要だと感じています。

# 潜在構造分析

data=data.frame(pattern=1:8,r1=c(1,1,1,0,1,0,0,0),r2=c(1,1,0,1,0,1,0,0),r3=c(1,0,1,1,0,0,1,0),total=c(61,44,89,20,71,49,39,131))

P=p=array(0,dim=c(2,2))

P[1,1]=sum(data$total[data$r3==1])/sum(data$total)
P[1,2]=sum(data$total[data$r3==1 & data$r2==1])/sum(data$total)
P[2,1]=sum(data$total[data$r3==1 & data$r1==1])/sum(data$total)
P[2,2]=sum(data$total[data$r3==1 & data$r1==1 & data$r2==1])/sum(data$total)

p[1,1]=1
p[1,2]=sum(data$total[data$r2==1])/sum(data$total)
p[2,1]=sum(data$total[data$r1==1])/sum(data$total)
p[2,2]=sum(data$total[data$r1==1 & data$r2==1])/sum(data$total)

func=function(s){
  
  matrix=P-p*s
  
  value=matrix[1,1]*matrix[2,2]-matrix[1,2]*matrix[2,1]
  
  return(value)
  
}

func_data=data.frame(s=seq(-1,1,0.001),values=0)

for(j in 1:nrow(func_data)){
  
func_data$values[j]=func(func_data$s[j])
  
}

# newton-raphson

values=c(-100,100);solution=c()

for(k in 1:2){

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=values[k]
for(j in 1:nrow(newton_raphson)){
f<-func(value)
df<-(func(value+0.0001)-func(value))/(0.0001)

  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value

  newton_raphson$value[j]=value
}

solution=c(solution,newton_raphson$value[nrow(newton_raphson)])

}

solution=sort(solution,decreasing = T)

x=array(0,dim=c(2,2))

for(k in 1:length(solution)){
  
matrix=P-solution[k]*p 

matrix=matrix[2,];

matrix[1]=-matrix[1];

x[k,]=matrix
  
}

x=t(x);

x_data=data.frame(num=1:nrow(x),x)

x_data=x_data[order(x_data$num,decreasing=T),]

x=x_data[,2:ncol(x_data)]

E_x=diag( 1/solve(x)[,1],2)

pai2=E_x%*%solve(x)

y=array(0,dim=c(2,2))

for(k in 1:length(solution)){
  
matrix=t(P)-solution[k]*t(p) 

matrix=matrix[2,];

matrix[1]=-matrix[1];

y[k,]=matrix
  
}

y=t(y);

y_data=data.frame(num=1:nrow(y),y)

y_data=y_data[order(y_data$num,decreasing=T),]

y=y_data[,2:ncol(y_data)]

E_y=diag( 1/solve(y)[,1],2)

pai1=E_y%*%solve(y)

M=solve(t(pai1))%*%p%*%solve(pai2)

result=data.frame(潜在クラス=c("好意的","非好意的"),クラスの確率=diag(M),r1=pai1[,2],r2=pai2[,2],r3=solution)

# 潜在構造分析(改良)

data=data.frame(pattern=1:8,r1=c(1,1,1,0,1,0,0,0),r2=c(1,1,0,1,0,1,0,0),r3=c(1,0,1,1,0,0,1,0),total=c(61,44,89,20,71,49,39,131))

P=p=array(0,dim=c(2,2))

P[1,1]=sum(data$total[data$r3==1])/sum(data$total)
P[1,2]=sum(data$total[data$r3==1 & data$r2==1])/sum(data$total)
P[2,1]=sum(data$total[data$r3==1 & data$r1==1])/sum(data$total)
P[2,2]=sum(data$total[data$r3==1 & data$r1==1 & data$r2==1])/sum(data$total)

p[1,1]=1
p[1,2]=sum(data$total[data$r2==1])/sum(data$total)
p[2,1]=sum(data$total[data$r1==1])/sum(data$total)
p[2,2]=sum(data$total[data$r1==1 & data$r2==1])/sum(data$total)

pisx=eigen(solve(p)%*%P)$values
X=eigen(solve(p)%*%P)$vectors

(P-pisx[1]*p)%*%X[,1]
(P-pisx[2]*p)%*%X[,2]

E_x=diag(c(1/solve(X)[,1]))
pai2=E_x%*%solve(X)


pisy=eigen(solve(t(p))%*%t(P))$values
Y=eigen(solve(t(p))%*%t(P))$vectors

(t(P)-pisy[1]*t(p))%*%Y[,1]
(t(P)-pisy[2]*t(p))%*%Y[,2]

E_y=diag(c(1/solve(Y)[,1]))
pai1=E_y%*%solve(Y)

M=solve(t(pai1))%*%p%*%solve(pai2)

result=data.frame(潜在クラス=c("好意的","非好意的"),クラスの確率=diag(M),r1=pai1[,2],r2=pai2[,2],r3=pisx)