RでタグチメソッドのMTA法を組んでirisで試してみる

水着よりもいまは浴衣姿を拝みたい.


あっ, どーも僕です.

みなさんこちらの記事見ましたか?
Webデータ分析&データサイエンスで役立つ統計学・機械学習系の分析手法10選 - 道玄坂で働くデータサイエンティストのブログ

Rをはじめて一年あまりだそうです. すごいですね.

わたしなんてもうまる二年Rを使っているのに半分くらいしか数学的にもR的にもわかりませんでした. お恥ずかしい.

....

まっ, それはおいといて.
記事を読んでいて思ったのですが, Webの世界ではあまりタグチメソッドは使われていないようですね. タグチメソッドも結構使えそうな気もするだけどなあ.

ってことで, 今回はタグチメソッドのMTA法を組んでみてirisで試してみました.
MTA法はパターン認識の一種で, 何がいいかはわかっているけど, 何が悪いかわかっていないときに有効です.
なにでパターン認識するかというと正しいと思うものの空間からの距離(違うけど)です. 距離が小さいほどいいものってことですね. まぁ, イメージとしてはマハラノビス距離の拡張です.
数学的に説明するのは大変なのですのでわかりやすい参考書を上げておきますね.

よくわかるMTシステム―品質工学によるパターン認識の新技術

よくわかるMTシステム―品質工学によるパターン認識の新技術


組んだといっても, 0から組んだのでわなくて, 卒論のときMTA法を援用したので, そのとき組んだのを見やすいように書き直しました.
ぱぱぱっとできるかなと思ったのですが, 当時使えた記法が使えなかったり, 余因子行列が半正定値にならなかったりとちょっと大変だったかな.



本題に入ると, 今回はirisの"setosa" で作った単位空間からの距離でシミュレーションしてみました. 結果がこちらです.
縦軸に距離で, 横軸がデータ番号です. いろは正しいものをつけています. きっちり距離がわかれていますね. *1

f:id:aaaazzzz036:20130622103707j:plain


# **************************************************
# MTA法による総合推定値の算出
# **************************************************
# MTA法による総合推定値 
computeMTAEstimateValue <- function (UnitSpaceData, 
                                     signalData, 
                                     Data) 
{
  # UnitSapceData : 単位空間サンプル
  # signalData    : 信号サンプル(最後列が信号値)
  # Data          : 評価サンプル
  
  # 追加パッケージ
  require (LoopAnalyst)
  
  
  # サンプルの正規化
  signal <- signalData [ , ncol (signalData), drop = TRUE]
  l      <- length (signal)
  mu     <- colMeans (UnitSpaceData)
  sigma  <- apply (UnitSpaceData, 2, function (x) sqrt (var (x) * (length (x) - 1) / length (x)))
  modified.sigma <- sigma; modified.sigma [sigma < 1e-9] <- 1
  normalized.UnitSpaceData <- scale (UnitSpaceData, center = mu, scale = modified.sigma)
  normalized.signalData    <- scale (signalData,    center = c(mu, 0), scale = c(modified.sigma, 1))
  normalized.Data          <- scale (Data,          center = mu, scale = modified.sigma)

  # sigma = 0の項目の処理
  if (sum (sigma < 1e-9) != 0) {
    index1   <- which (sigma < 1e-9)
    r       <- crossprod (signal)
    L       <- crossprod (signal, normalized.signalData [, index1])
    ST      <- diag (crossprod (normalized.signalData [ , index1]))
    SB      <- L ** 2 / r
    Ve      <- (ST - SB) / (l - 1)
    beta.zeroVar <- L / r
    ita.zeroVar <- (SB - Ve) / (r * Ve)
    ita.zeroVar [ita.zeroVar < 0] <- 0
    
    # 評価サンプルの計算
    X1 <- sapply (seq_along (index1), function (i) normalized.Data[,index1[i]] / beta.zeroVar[i]) ** 2
    X1 <- X1 %*% ita.zeroVar
  } else {
    X1 <- NULL
    ita.zeroVar  <- NULL
  }
  
  # sigma != 0の項目の処理
  if (sum (sigma >= 1e-9) != 0) {
    index2 <- which (sigma >= 1e-9)
    k      <- length (index2)
    adjointMat.normalized.UnitSpaceData <- round (make.adjoint (cor (normalized.UnitSpaceData [ , index2])), 5)
    # 計算された余因子行列だと距離が負になるので絶対値をつけておく
    # 余因子行列が半正定値行列じゃないのはなぜだああああああああ
    mad   <- sqrt (abs (diag (crossprod (t (normalized.signalData [ , index2]), 
                                     adjointMat.normalized.UnitSpaceData %*% t (normalized.signalData [ , index2])))/ k))
    r     <- crossprod (signal)
    L     <- c (crossprod (signal, mad))
    ST    <- crossprod (mad)
    SB    <- L ** 2 / r
    Ve    <- (ST - SB) / (l - 1)
    beta.notzeroVar <- L / r
    ita.notzeroVar <- (SB - Ve) / (r * Ve)
    ita.notzeroVar [ita.notzeroVar < 0] <- 0
    
    # 評価サンプルの計算
    # 計算された余因子行列だと距離が負になるので絶対値をつけておく
    # 余因子行列が半正定値行列じゃないのはなぜだああああああああ
    D <- sqrt (abs (diag (crossprod (t (normalized.Data [, index2]), 
                                adjointMat.normalized.UnitSpaceData %*% t (normalized.Data [, index2]))) / k))
    X2 <- ita.notzeroVar * ((D / beta.notzeroVar) ** 2)
  } else {
    X2 <- NULL
    ita.notzeroVar  <- NULL
  }
  # 推定値
  X       <- rowSums (matrix (c (X1, X2), nrow = nrow (Data), ncol = 2))
  sum.Ita <- sum (c (ita.notzeroVar, ita.notzeroVar))
  result  <- sqrt (X / sum.Ita)

  return (result)
  
}

# 動作確認 with 適当に変換したiris
myiris <- iris
myiris <- cbind (myiris[, -5], V2 = rep (8:10, each = 50))
myiris <- cbind (myiris, V3 = c(iris[,5]))
USP <- myiris [ 1:25, -6]                        # 単位サンプル
SD  <- myiris [ c(26:30, 76:80, 126:130),   ]    # 信号サンプル
Dat <- myiris [ -c(1:30, 76:80, 126:130), -6]    # 評価サンプル
rtn <- computeMTAEstimateValue (USP, SD, Dat)

# グラフ化
COL <- rep (brewer.pal(3, "Set2"), each = 50)[-c(1:30, 76:80, 126:130)]
def.par <- par(no.readonly = TRUE) # デフォルトの環境を保存
par (mai = c(1.5, 1.5, .5, .5))
plot (rtn , type = "h", col = COL, cex = 3, pch = 19, lwd = 3, 
      cex.axis = 1.8, cex.lab = 1.8, ylab = "MTAEstimateValue", xlab = "Index", 
      tcl = .5)
legend ("bottomright", legend = c(levels (iris[,5])), 
        cex = 2.5, pch = 19, col = unique(COL), 
        bty = "n", title = "Species")
par(def.par)            # デフォルトの環境を復帰



*1:距離をプロットするのはおかしな話でしたので, 棒グラフにしました