Help us understand the problem. What is going on with this article?

# KL数量化

More than 1 year has passed since last update.

n=20;m=10

X=array(0,dim=c(n,m))

for(j in 1:ncol(X)){

X[,j]=rpois(n,lambda=sample(10,1))

}

ave_x=apply(X,2,mean)

for(j in 1:m){

X[,j]=(X[,j]-ave_x[j])

}

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

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

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

d[i,j]=sum(abs(X[i,]-X[j,])^2)

s[i,j]=sum(abs(X[i,]-X[j,])^2)

}
}

e=-s

H=array(0,dim=c(nrow(e),ncol(e)))

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

if(j!=i){

H[j,i]=e[i,j]+e[j,i]

}
}
}

alpha=-max(H)

H=H+alpha;diag(H)=0

for(j in 1:ncol(H)){

diag(H)[j]=-sum(H[j,])

}

x_vec=eigen(H)\$vectors[,1];y_vec=eigen(H)\$vectors[,2]

lambda=eigen(H)\$values[1];mu=eigen(H)\$values[2]

ite=100

for(t in 1:ite){

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

vec=array(0,dim=c(2,1))

for(j in 1:nrow(X)){
for(k in 1:nrow(X)){

matrix[1,1]=matrix[1,1]+sum((x_vec[j]-x_vec[k])^4)

matrix[1,2]=matrix[1,2]+sum((x_vec[j]-x_vec[k])^2)

matrix[2,1]=matrix[2,1]+sum((x_vec[j]-x_vec[k])^2)

matrix[2,2]=nrow(X)*(nrow(X)-1)

vec[1,1]=vec[1,1]+s[j,k]*sum((x_vec[j]-x_vec[k])^2)

vec[2,1]=vec[2,1]+s[j,k]

}
}

KL=c(solve(matrix,vec))

x1=c(sqrt(KL[1]))*c(x_vec)

L1=KL[2]

r=s-L1;

guzai=array(0,dim=c(length(x1),length(x1)))

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

guzai[i,j]=x1[i]-x1[j]

}
}

A=array(0,dim=c(length(x1),length(x1)))

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

A[i,j]=r[i,j]-3*guzai[i,j]^2

}
}

A_mat=A

for(l in 1:length(x1)){

A_mat[l,]=A_mat[l,]*x1

}

diag(A_mat)=0

A_mat_sub=A_mat

A_sub=A;diag(A_sub)=0

diag(A_mat)=-apply(A_sub,1,sum)*x1

term=apply(r*guzai,1,sum)-apply(guzai^3,1,sum)

A_mat[1,]=x1;term[1]=0

dx=solve(A_mat,term)

x2=x1*(1+dx)

print(max(abs(dx)))

x_vec=x2

}

n=20;m=10

X=array(0,dim=c(n,m))

for(j in 1:ncol(X)){

X[,j]=rpois(n,lambda=sample(10,1))

}

ave_x=apply(X,2,mean)

for(j in 1:m){

X[,j]=(X[,j]-ave_x[j])

}

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

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

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

d[i,j]=sum(abs(X[i,]-X[j,])^2)

s[i,j]=sum(abs(X[i,]-X[j,])^2)

}
}

e=-s

H=array(0,dim=c(nrow(e),ncol(e)))

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

if(j!=i){

H[j,i]=e[i,j]+e[j,i]

}
}
}

alpha=-max(H)

H=H+alpha;diag(H)=0

for(j in 1:ncol(H)){

diag(H)[j]=-sum(H[j,])

}

x_vec=eigen(H)\$vectors[,1];y_vec=eigen(H)\$vectors[,2]

lambda=eigen(H)\$values[1];mu=eigen(H)\$values[2]

ite=100

for(t in 1:ite){

matrix=array(0,dim=c(3,3))

vec=array(0,dim=c(3,1))

for(j in 1:nrow(X)){
for(k in 1:nrow(X)){

matrix[1,1]=matrix[1,1]+sum((x_vec[j]-x_vec[k])^4)

matrix[1,2]=matrix[1,2]+((y_vec[j]-y_vec[k])^2)*(x_vec[j]-x_vec[k])^2

matrix[1,3]=matrix[1,3]+sum((x_vec[j]-x_vec[k])^2)

matrix[2,1]=matrix[2,1]+((y_vec[j]-y_vec[k])^2)*(x_vec[j]-x_vec[k])^2

matrix[2,2]=matrix[2,2]+sum((y_vec[j]-y_vec[k])^4)

matrix[2,3]=matrix[2,3]+sum((y_vec[j]-y_vec[k])^2)

matrix[3,1]=matrix[3,1]+sum((x_vec[j]-x_vec[k])^2)

matrix[3,2]=matrix[3,2]+sum((y_vec[j]-y_vec[k])^2)

matrix[3,3]=n*(n-1)

vec[1,1]=vec[1,1]+s[j,k]*sum((x_vec[j]-x_vec[k])^2)

vec[2,1]=vec[2,1]+s[j,k]*sum((y_vec[j]-y_vec[k])^2)

vec[3,1]=vec[2,1]+s[j,k]

}
}

KL=c(solve(matrix,vec))

x1=c(sqrt(KL[1]))*c(x_vec)

y1=c(sqrt(KL[2]))*c(y_vec)

L1=KL[3]

r=s-L1;

guzai=array(0,dim=c(length(x1),length(x1)))

eta=array(0,dim=c(length(y1),length(y1)))

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

guzai[i,j]=x1[i]-x1[j]

}
}

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

eta[i,j]=y1[i]-y1[j]

}
}

A=array(0,dim=c(length(x1),length(x1)))

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

A[i,j]=r[i,j]-3*guzai[i,j]^2-eta[i,j]^2

}
}

A_mat=-A

for(l in 1:length(x1)){

A_mat[l,]=A_mat[l,]*x1

}

diag(A_mat)=0

A_sub=A;diag(A_sub)=0

diag(A_mat)=apply(A_sub,1,sum)*x1

A_sub_x=A_mat

A=array(0,dim=c(length(y1),length(y1)))

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

A[i,j]=2*eta[i,j]*guzai[i,j]

}
}

A_mat=A

for(l in 1:length(x1)){

A_mat[l,]=A_mat[l,]*y1

}

diag(A_mat)=0

A_sub=-A;diag(A_sub)=0

diag(A_mat)=apply(A_sub,1,sum)*y1

A_sub_y=A_mat

A_sub_xy1=cbind(A_sub_x,A_sub_y)

A_sub_xy1[1,]=c(rep(1,length(x1)),rep(0,length(y1)))

A=array(0,dim=c(length(x1),length(x1)))

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

A[i,j]=r[i,j]-3*eta[i,j]^2-guzai[i,j]^2

}
}

A_mat=-A

for(l in 1:length(x1)){

A_mat[l,]=A_mat[l,]*y1

}

diag(A_mat)=0

A_sub=A;diag(A_sub)=0

diag(A_mat)=apply(A_sub,1,sum)*y1

A_sub_x=A_mat

A=array(0,dim=c(length(y1),length(y1)))

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

A[i,j]=2*eta[i,j]*guzai[i,j]

}
}

A_mat=A

for(l in 1:length(x1)){

A_mat[l,]=A_mat[l,]*x1

}

diag(A_mat)=0

A_sub=-A;diag(A_sub)=0

diag(A_mat)=apply(A_sub,1,sum)*y1

A_sub_y=A_mat

A_sub_xy2=cbind(A_sub_y,A_sub_x)

A_sub_xy2[1,]=c(rep(0,length(x1)),rep(1,length(y1)))

A_sub=rbind(A_sub_xy1,A_sub_xy2)

term1=apply(guzai^3,1,sum)+apply(guzai*eta^2,1,sum)-apply(r*guzai,1,sum)

term2=apply(eta^3,1,sum)+apply(eta*guzai^2,1,sum)-apply(r*eta,1,sum)

term1[1]=0;term2[1]=0

term=c(term1,term2)

dx=solve(A_sub,term)

x2=x1*(1+dx[1:length(x1)])

y2=y1*(1+dx[(length(x1)+1):(length(x1)+length(y1))])

print(max(abs(dx)))

x_vec=x2;y_vec=y2

}

Why not register and get more from Qiita?
1. We will deliver articles that match you
By following users and tags, you can catch up information on technical fields that you are interested in as a whole
2. you can read useful information later efficiently
By "stocking" the articles you like, you can search right away