CH1R<-function(y=ch1)
{
y=as.matrix(y)
Ny <- ncol(y)
NLV <- 6
N <- nrow(y)
#for(i in N:1){if(max(y[i,  ]) > 10) y <- y[ - i,  ]}
N <- nrow(y)
Y <- sweep(sweep(y, 2, apply(y, 2, mean)), 2, sqrt(apply(y, 2, var)), "/")
INR <- c(0,0, 0, 0, 0, 1,
1,0, 0, 0, 0, 1,
1,1, 0, 0, 0, 1,
1, 1, 1, 0, 0,1,
0, 0, 0, 1, 0,0,
0, 0, 0, 0, 0, 1)
INR <- matrix(INR, NLV, NLV, byrow = T)
OUTR <- c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,
  0, 0, 0, 0, 1, 1, 1, 1, 0, 0,0,0,0,0,0,0,0,0,
  0, 0, 0, 0, 0, 0, 0, 0, 1, 1,0,0,0,0,0,0,0,0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0,1,-1,-1,-1,0,0,0,0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,1,1,0,0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,1,1)
OUTR <- matrix(OUTR, NLV, Ny, byrow = T)
fs <- matrix(1, NLV, 1)

V <- OUTR
W2 <- V
W1 <- V
f <- sqrt(N)/sqrt(apply((V %*% t(Y))^2, 1, sum))
f <- f * fs
W2 <- sweep(V, 1, f, "*")
CRITERIA <- 1e-005
MAXM <- 10
A <- 0
while(MAXM > CRITERIA) {
A <- A + 1
W1 <- W2
X <- Y %*% t(W1)
#print(X[1:3,]);print(dim(X));print(dim(Y));print(Y[1:3,]);print(dim(t(W1)));print(t(W1)[1:3,])
CORX <- cor(X)
TRASH <- (INR == 1) + (t(INR) == 1)
TEM <- matrix(0, NLV, NLV)
diag(TEM) <- diag(CORX)
TRASH <- (TRASH >= 1) * (CORX - TEM)
SIGNCOR <- (TRASH > 0) - (TRASH < 0)
SWS <- X %*% SIGNCOR
for(I in 1:NLV)
    for(J in 1:Ny)
if(OUTR[I, J] != 0)
    V[I, J] <- sum(SWS[, I] * Y[, J])/sum(SWS[, I] * SWS[, I])
f <- sqrt(N)/sqrt(apply((V %*% t(Y))^2, 1, sum))
f <- f * fs
W2 <- sweep(V, 1, f, "*")
MAXM <- max((W2 != 0) * abs((W2 - W1)/(W2 + (W2 == 0))))
}

X <- Y %*% t(W2)
OUTCOE <- matrix(0, Ny, 1)
for(I in 1:NLV)
for(J in 1:Ny)
if(OUTR[I, J] != 0) OUTCOE[J] <- sum(X[, I] * Y[, J])/sum(X[, I] * X[, I])
INCOE <- matrix(0, NLV, NLV)
for(I in 1:NLV) {
IDV <- matrix(0, N, 1)
if(sum(INR[I,  ]) != 0) {
IDVNO <- 0
for(J in 1:NLV) {
if(INR[I, J] != 0) {
  IDV <- cbind(IDV, X[, J])
  IDVNO <- cbind(IDVNO, J)
}
}
IDV2 <- IDV[, 2:ncol(IDV)]
IDVNO <- IDVNO[, 2:ncol(IDV)]
B <- solve(t(IDV2) %*% IDV2) %*% t(IDV2) %*% as.vector(X[, I])
}
K <- 1
for(J in 1:NLV) {
if(INR[I, J] != 0) {
INCOE[I, J] <- B[K]
K <- K + 1
}
}
}

HH<-function(W,j) {k<-NULL;for (i in 1:dim(W)[1])  if (W[i,j]==1) k<-c(k,i);k}
ID <- function(W2, y, i, jk)
{
#i is #col of W2', jk are #rows;W2:10x5
W <- t(W2[i, jk])
if(length(jk) > 1)
XBA <- apply(y[, jk], 2, mean)
else XBA <- mean(y[, jk])
if(length(jk) > 1)
MINX <- apply(y[, jk], 2, min)
else MINX <- min(y[, jk])
if(length(jk) > 1)
MAXX <- apply(y[, jk], 2, max)
else MAXX <- max(y[, jk])
MMI <- sum(W * MINX)
MMA <- sum(W * MAXX)
LOYI <- ((sum(W * XBA) - MMI) * 100)/(MMA - MMI)
LOYI
}
IDEX <- rep(0, NLV)
for(i in 1:NLV) {
jk <- HH(t(OUTR), i)
IDEX[i] <- ID(W2, y, i, jk)# It was IDEX[i] <- ID(W2, Y, i, jk) in 2001
}
#list(W2 = W2, N = N, INCOE = INCOE, OUTCOE = OUTCOE, IDEX = IDEX)
list(N = N, INCOE = INCOE, OUTCOE =t(OUTCOE), IDEX = IDEX)
}
大家讨论一下，对于数据y，有什么要求？
程序的一般性？
对于程序中，不同的问题，需要修改的命令是？
欢迎举例说明用法，和结果的解释。