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

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

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

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件法でデータ取ってる研究全部やり直したいぐらいだよ,ほんとに。

R
ベイズ塾合宿2019参戦

今年も塾合宿の季節がやってまいりました。会場は浜松!   史上最大規模の,33人が参加。二 …

分類できない何か
日本心理学会の引用スタイルjecon_jpa.styを作った

PDFが静的でリッチではないフォーマットだ,という意見はよくわかるけど,記録として残るものは静的であ …

日記
文字数

小学生の頃の作文の授業。「たくさん書けば書くほど良い」というルールだった(と勘違いした)ので、遠足の …