交互平均法のコード

西里先生のこの本を見ていると，自分でコードを書いてちゃんと理解したくなった。

[amazonjs asin=”4563052183″ locale=”JP” title=”行動科学のためのデータ解析―情報把握に適した方法の利用”]

ということで，2章にある双対尺度法の基本アルゴリズム，交互平均法についてのコードがこちら。

recipro_ave <- function(dat,eps=1e-10,iter.max=1000,start.vec=NULL,proc=F){
  nc <- ncol(dat)
  nr <- nrow(dat)
  tmp <- 0
  FLG <- FALSE
  iter <- 0
  if(is.null(start.vec)){
    u <- rnorm(nc)
  }else{
    if(length(start.vec)==nc){
      u <- start.vec      
    }else{
      warning("vec size is not correct. Starting val. is changed.")
      u <- rnorm(nc)      
    }
  }
  while(FLG==FALSE){
    #algorithm
    v <- u %*% t(dat)
    v <- v / rowSums(dat)
    av <- (rowSums(dat) %*% t(v))/sum(rowSums(dat))
    v <- v - (av * rep(1,nr))
    gy <- max(abs(v))
    v <- v/gy
    
    u <- v %*% dat
    u <- u / colSums(dat)
    av <- (colSums(dat) %*% t(u))/sum(colSums(dat))
    u <- u - (av *rep(1,nc))
    gx <- max(abs(u))
    u <- u/gx
    
    #process print
    if(proc){
      cat("iter  ",iter,"\n")
      cat("eigen col:",gy," col weight:",v,"\n")
      cat("eigen row:",gx," row weight:",u,"\n")
    }
    #conv. check
    if(abs(tmp-gx)<eps){FLG <- TRUE}else{tmp <- gx}
    if(iter>iter.max){FLG <- TRUE}else{iter <- iter + 1}

  }

  # return
  
  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(colSums(dat))/colSums(dat)%*%t(u^2)) #multiplied constant
  nv.c <- sqrt(sum(rowSums(dat))/rowSums(dat)%*%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){
    warning("Warning::did not converge.")
  }else{
    return(list(eigen.col=gy,eigen.row=gx,cor.ratio=eta2,singular.val=sqrt(eta2),projected.col=pv,projected.row=pu,normed.col=nv,normed.row=nu))
  }

}

関数名はrecipro_aveで，クロス集計表のデータを渡してやれば結果を返すようになっている。
オプションとして，収束基準eps(デフォルトは1e-10),回転上限iter.max(デフォルトは1000),最初につける列重みベクトルstart.vec(デフォルトはNULLで，正規乱数を入れるようになっている),各ステップを出力するかどうかproc(デフォルトはFALSE)が決められるようになっている。
授業などでデモンストレーションするときは，start.vecにc(1,2,3)などを渡してproc=Tとすると，一手ずつ計算結果が出力されるのでよろしいかと。