双対尺度法を自作してみた

前回の続き。
交互平均法で出てくるのは第一成分だけなので，複数の次元を出すために残差行列を計算してさらに交互平均法をする，という作業を繰り返す必要がある。

実装してみたコードがこちら。このコードの中に交互平均法のアルゴリズムも含まれています。

DualScaling <- function(dat,eps=1e-10,iter.max=1000){

  nc <- ncol(dat)
  nr <- nrow(dat)
  cname <- colnames(dat)
  rname <- rownames(dat)
  mg.col <- colSums(dat)
  mg.row <- rowSums(dat)
  mg.t <- sum(mg.row)
  dm <- min(nc,nr)-1
  normed.col <- matrix(ncol=nc,nrow=dm)
  projected.col <- matrix(ncol=nc,nrow=dm)
  normed.row <- matrix(ncol=nr,nrow=dm)
  projected.row <- matrix(ncol=nr,nrow=dm)
  singular <- c()
  
  dat.tmp <- dat
  
  # 0-exp
  for( i in 1:nrow(dat)){
    for( j in 1:ncol(dat)){
      dat.tmp[i,j] <- dat.tmp[i,j]- (mg.row[i]*mg.col[j]/mg.t)
    }  
  }
  
  # 1-n
  for(i in 1:dm){
    FLG <- FALSE
    iter <- 0
    u <- rnorm(nc)
    tmp <- 0
    while(FLG==FALSE){
      #algorithm
      v <- u %*% t(dat.tmp)
      v <- v / mg.row
      av <- (mg.row %*% t(v))/mg.t
      v <- v - (av * rep(1,nr))
      gy <- max(abs(v))
      v <- v/gy
      
      u <- v %*% dat.tmp
      u <- u / mg.col
      av <- (mg.col %*% t(u))/mg.t
      u <- u - (av *rep(1,nc))
      gx <- max(abs(u))
      u <- u/gx
      
      #conv. check
      if(abs(tmp-gx)<eps){FLG <- TRUE}else{tmp <- gx}
      if(iter>iter.max){FLG <- TRUE}else{iter <- iter + 1}
      
    }   
  
  
    if(u[which.max(abs(u))]<0) u *-1
    if(v[which.max(abs(v))]<0) v *-1
    
    eta2 <- gx * gy  # Correlation ratio
    nu.c <- sqrt(sum(mg.col)/mg.col%*%t(u^2)) #multiplied constant
    nv.c <- sqrt(sum(mg.row)/mg.row%*%t(v^2))
    nu <- nu.c %*% u #normed u
    nv <- nv.c %*% v #normed v
    pu <- nu * sqrt(eta2)
    pv <- nv * sqrt(eta2)
    
    if(iter>iter.max){
      cat("Warning::did not converge.")
    }else{
      singular[i] <- sqrt(eta2)
      normed.col[i,] <- nu
      normed.row[i,] <- nv
      projected.col[i,] <- pu
      projected.row[i,] <- pv        
    }
    
    # residual Matrix
    dat.tmp <- dat.tmp-(mg.row %*% t(mg.col)/mg.t)*((sqrt(eta2)*t(nv)%*%nu))
  }
  
  # Correct dims
  delta_k <- singular^2 / sum(singular^2) * 100
  cum_delta_k <- c()
  cum_delta_k[1] <- delta_k[1]
  for(i in 2:dm)
    cum_delta_k[i] <- cum_delta_k[i-1]+delta_k[i]
  
  dm.cr <- which.max(cum_delta_k)
  if(dm.cr<dm){
    singular <- singular[-(dm.cr+1:dm)]
    delta_k <- delta_k[-(dm.cr+1:dm)]
    cum_delta_k <- cum_delta_k[-(dm.cr+1:dm)]
    normed.col <- normed.col[-(dm.cr+1:dm),] 
    normed.row <- normed.row[-(dm.cr+1:dm),]
    projected.row <- projected.row[-(dm.cr+1:dm),]
    projected.col <- projected.col[-(dm.cr+1:dm),]
  }
  
  # var names
  colnames(normed.col) <- cname
  colnames(projected.col) <- cname
  colnames(normed.row) <- rname
  colnames(projected.row) <- rname
  
  return(list(ndims=dm.cr,sing=singular,eig=singular^2,deltaK=delta_k,cum_dk=cum_delta_k,nc=normed.col,nr=normed.row,pc=projected.col,pr=projected.row))
}

ちなみに，青木先生のサイトにも双対尺度法のコードがあって，こちらはRの固有値分解関数などをつかっているので，より安心安全にお使いいただけます。こっちのコードは交互平均法による習作ですから。

あと，西里先生の研究チームによる双対尺度法のRパッケージ,dualScaleというのもある。

dualScale: Dual Scaling Analysis of Multiple Choice Data

これは多肢選択法にdual scalingをするdsMC,強制分類法をするdsFCという二つの関数しかないけれども，どちらもとても魅力的な手法なので，パッケージで供給されているのはとてもありがたい。今後このパッケージがどんどん発展してくれるといいですね。

dsMCをつかって5件法でデータ取ってる研究全部やり直したいぐらいだよ，ほんとに。