数量化法の基礎 岩坪秀一先生 オーム社

#数量化三類からマルコフ連鎖へ

#要確認

mat=t(matrix(c(1,1,1,1,0,1,1,1,0,1,1,1),ncol=4))

PXY=mat/sum(mat)

PX=diag(apply(mat,1,sum),length(apply(mat,1,sum)))/sum(mat)

PY=diag(apply(mat,2,sum),length(apply(mat,2,sum)))/sum(mat)

PYX=t(PXY)

sqrt_PY=((PY)^(-1/2));sqrt_PY[sqrt_PY==Inf]=0

A_star=sqrt_PY%*%PYX%*%solve(PX)%*%PXY%*%sqrt_PY

lambda=eigen(A_star)$values

eigen_vectors=eigen(A_star)$vectors

z2=eigen_vectors[,2];z3=eigen_vectors[,3]


y2=sqrt_PY%*%z2;y3=sqrt_PY%*%z3

x2=solve(PX)%*%PXY%*%sqrt_PY%*%z2/sqrt(lambda[2])

x3=solve(PX)%*%PXY%*%sqrt_PY%*%z3/sqrt(lambda[3])

A=solve(PY)%*%PYX%*%solve(PX)%*%PXY

#定常確率の確認

times=100

trans=A

for(j in 1:times){

trans=trans%*%A

}

stationary_prob=PY%*%rep(1,ncol(PY))

#カイ二乗検定(2変数の独立性)

floor((sum(diag(solve(PY)%*%PYX%*%solve(PX)%*%PXY))-1)*10^4)/(10^4)==floor((lambda[2]+lambda[3])*(10^4))/10^4

N2=length(lambda)

kai2=length(lambda)*(sum(diag(solve(PY)%*%PYX%*%solve(PX)%*%PXY))-1)

qchisq(1-0.05/2,length(lambda))

#エントロピーを最大にする遷移確率を無効グラフから生成する

A=matrix(sample(c(0,1),6*6,replace = T),ncol=6)

lambda0=max(Re(eigen(A)$values))

eigen_vectors=Re(eigen(A)$vectors[,1])

P=array(0,dim=c(nrow(A),ncol(A)))

for(i in 1:nrow(P)){
for(j in 1:ncol(P)){  

P[i,j]=(1/lambda0)*A[i,j]*eigen_vectors[j]/eigen_vectors[i] 

}
}

entropy=(-P*log(P))

entropy=sum(entropy[is.nan(entropy)==F])


times=1000

stationary=P

for(j in 1:times){

stationary=stationary%*%P  

}