Rで機械学習（SVM）

人間の直観的な認識力を手軽に模倣するフレームワークである機械学習がもてはやされるようになってもう10年以上経つと思いますが、世間では第3次人工知能ブームでまだまだ注目されているようです。扱いやすいパッケージが拡充されているのはもちろん『情報処理』2015年5月号に載った「機械学習のための数学」のように理屈の方も簡潔で分かりやすい紹介が増えてきました。もう何か付け足すべき事など無さそうですが、Rで説明用のコードを書く必要があったので、ひっそりと紹介したいと思います。
機械学習と言っても手法のバリエーションは豊富なのですが、SVM(Support Vector Machine)が試すには手頃なようです。Karatzoglou, Meyer and Hornik (2005)をはじめとしてパッケージの解説が豊富です。ブログでRでの著名パッケージの使い方も良く紹介されていますね。

1. 学習データと評価データの作成

機械学習は、第一段階で分類情報有りの学習データから分類器を構築し、第二段階で構成した分類器に分類情報無しの他のデータの分類作業をさせます。今回は学習データと、分類器の性能評価のための評価データを用意します。
人間の直観力を示していそうな男女の体格データを用いました。遺伝子をチェックしないと厳密には性別は判別できませんが、ほとんどの場合は外見だけで性別を判別しているものです。残念ながら十分な量の男女の体格データをネット上で発見できなかったので、今回はブートストラップで増やしたデータを用います。

#
# 以下のページから男性と女性の胸囲と臀囲のサイズを拝借
#
# http://homepage3.nifty.com/orangejuice/body4.html
# http://homepage3.nifty.com/orangejuice/body6.html
#
man <- data.frame(
  chest = c(84.8, 87.3, 85.2, 83.4, 83.9, 78.3),
  hip = c(91.3, 92.7, 89.7, 88.1, 88.3, 82.7)
)

woman <- data.frame(
  chest = c(85.0, 89.1, 83.4, 84.1, 81.3, 81.8, 82.1),
  hip = c(95.5, 96.6, 93.3, 92.1, 90.1, 89.6, 89.8)
)

#
# bootstrapでデータ作成
#
library(boot)

create_data <- function(n){
  #
  # データ合成関数
  # df:データフレーム, i:ランダムに選ばれる添字番号
  #
  statistic <- function(df, i){
    # iだけをデータフレームから抜き出す
    s <- df[i, ]
    # 平均をとって戻す
    c(mean(s$chest), mean(s$hip))
  }

  #
  # 男性と女性のそれぞれのデータからサンプル作成
  #
  r_m_boot <- boot(man, statistic, n/2)
  r_w_boot <- boot(woman, statistic, n/2)

  #
  # 男性データと女性データを合成して戻す
  #
  df <- data.frame(
    gender = factor(c(rep("M", n/2), rep("W", n/2))),
    chest = c(r_m_boot$t[,1], r_w_boot$t[,1]),
    hip = c(r_m_boot$t[,2], r_w_boot$t[,2])
  )
}

set.seed(20160614)
s_materials <- create_data(10) # 学習データ・サイズは10
e_materials <- create_data(100) # 評価データ・サイズは100

2. 学習データから分類器を構築

作業を表現すると堅苦しいわけですが、実際はe1071パッケージのおかげで3行で済みます。

library("e1071")
set.seed(722)

#
# 性別（gender）と、胸囲（chest）と臀囲（hip）の関係を学習
#
svp <- svm(gender~chest+hip, data=s_materials)

オプションを工夫すると精度が上がる*1わけですが、今回は作業の流れを追うためにしません。

3. 評価データで構築された分類器の性能評価

データが単純かつSVM向きなせいで、正答率97.0%と高性能になっています。

#
# 性別の推定値をans列、本当の性別との比較結果をchk列に入れる
#
e_materials$ans <- predict(svp, e_materials)
e_materials$chk <- e_materials$ans == e_materials$gender

#
# 正答率を見てみる
#
sprintf("正答率: %.1f%%", 100*sum(e_materials$chk)/length(e_materials$chk))

4. 評価データで構築された分類器の性能評価

以下のようにプロットすると分類器が男性だと判断する水色領域と、女性だと判断するピンク色領域の上に、黒xとoが真の男性データと、赤xとoが女性データを表示してくれます*2。

# 学習データでプロット
plot(svp, s_materials, chest~hip, col=c("lightblue1", "pink1"))

上の学習データの当てはまりが良いのは当然なので、評価データを見て見ましょう。なお、実用では正解が分かりませんが、今回の評価データにはgender列に正解が書いてあるので答え合わせになります。

# 評価データでプロット
plot(svp, e_materials, chest~hip, col=c("lightblue1", "pink1"))

hipが90から91ぐらいの女性データが男性に分類されてしまっていますが、全体としては悪くない当てはまりのようです。

5. まとめ

画像処理にも使われるのですが、システム構築としては特徴量の取得と整備が重要で、機械学習自体はちょっと前からコモディティ化しつつあるのが、実際に練習してみると良く分かりますね。
なお、機械学習分野でやられている手の込んだものではない、統計学で習うような素朴なロジスティック回帰で同じことが出来ないか試してみたのですが、分類自体は似たような精度で出来たものの、「数値的に 0 か 1 である確率が生じました」と推定量が不安定になったり*3、初期パラメーターによっては収束しなかったりするので、思いつきで実装するのはやめた方が良さそうです。

*1:例えば引数に、method = "C-classification", kernel = "radial", cost = 10, gamma = 0.1を加えると、評価データへの当てはまりが良くなります。なお、e1071パッケージには最適パラメーターを模索するための関数tune.svmが用意されていて、これで試行錯誤ができます。

*2: xはデータポイントかつサポートベクターとなります。

*3:完全分離可能なデータだと推定量が無限大になるので出るエラーだそうです。

結婚市場からの離脱を考慮した夫婦の社会的階層格差

@kabutoyama_taro氏に、「ネット界隈が想起する"上昇婚"のイメージ」は・・・数学的にありえないと言われた*1のですが、ネット界隈の仮定に基づけばあり得るので、シミュレーションをして例示したいと思います。数式を書いて検証すると煩雑になりますし、「ありえない」ことを否定するには具体例を示せば十分なので*2。

1. セットアップ

仮定は以下のようにします。社会的階層は順番で表現していますし、そう無理な仮定は置いていないと思います。ネット界隈の人が考えている一方で、@kabutoyama_taro氏が見落としているのは、ラベルを貼った【結婚市場からの離脱】と言う想定です。所謂、理想が高くて結婚できないケースです。

• それぞれN人の男女がいる
• 一夫一妻制
• それぞれ1からN番目まで社会的階層による序列がついている
• 序列番号の小さい順に男性が結婚していく。つまり、序列番号iの男性が結婚すると、序列番号i+1の男性が求婚をして回る
• 男性は可能な限り社会的序列とは独立に定まる容姿が優れた女性と結婚したがる
• 【結婚市場からの離脱】序列jの女性は、求婚してくる序列iの男性に対して、確率P(i/j)で結婚を承諾する

2. ソースコード

検証用に明示します。Nは300にしました。P(i/j)はプロビット分析などに準じて正規の累積分布関数を用いましたが、どのような分布でも大差は無いです。

set.seed(20160526)
N <- 300
beauty <- runif(N)
status_f <- rep(FALSE, N)
married <- numeric(N)

P <- function(a){
# iが小さいほど結婚を承諾される可能性が高いようにしている
  1 - pnorm(a)
}

for(i in 1:N){
  best_b <- 0
  best_f <- 0
  for(j in 1:N){
    if(!status_f[j] && P(i/j)>runif(1) && best_b < beauty[j]){
      best_f <- j
      best_b <- beauty[j]
    }
  }
  if(0 < best_f){
    status_f[best_f] <- TRUE
    married[i] <- best_f
  }
}

gap <- (1:N)[married>0] - married[married>0]

hist(gap, breaks=11, main=sprintf("それぞれ%d人の男女から%d組が婚姻", N, length(gap)), xlab="男性の序列番号 - 女性の序列番号", ylab="頻度", xlim=c(-N, N))
abline(v = mean(gap), lty=2)

3. 結果

歪みは思ったような感じではありませんでしたが、平均値は-52.37959なので夫婦の階層格差が保存されています。また、男女それぞれ59人が結婚できていません。理由はあれこれ考えるまででも無く、【結婚市場からの離脱】が効いています*3。i/j<1ならば1、i/j>1ならば0を返す指示関数を確率関数にとれば、シミュレーションするまでもありません。
なお、男性が容姿ではなく社会的階層を女性に求めるように設定を変えると、夫婦の階層格差とその分散は一気に縮まる事になります*4。お見合いで家柄をあわせるような世界は、数理的に意味があったわけですね。

*1: Twitterで私が示した図の符号が間違って逆転していたので、修正しています。

*2: 現実にどうなっているかは、本稿の関心対象外なので悪しからず。

*3: 全員が結婚するのであれば、左右対称の分布になります。

*4: ソースコードは省略しますが、プロットすると以下のような結果になります。

Rの並列処理（doParallel/foreach）で繰り返し処理を中断する方法

確認してみたのですが、試行錯誤してみた結果、結論は「ありません」でした*1。仕方が無いので、OpenMP 4.0以前の代替措置*2と同じ事をしてみたいと思います。速度的には改善になるはずです。

library("doParallel")

cl <- makeCluster(detectCores())
registerDoParallel(cl)

f <- TRUE
# 処理を if(f){ ... } で囲み、f==TRUEの間だけ計算するようにする
r <- foreach(i=1:100, .combine='c') %dopar% if(f){
# √計算するだけの処理
  v <- sqrt(i)
  if(v > 6){
# 目標達成をしたら f にFALSEを入れる
    f <- FALSE
  }
  v
}

# 平行処理終了
stopCluster(cl)

# 結果を確認する。f <- FALSE をコメントアウトして比較すると違いが分かるはず。
r

1から100までのルート計算をしていくループですが、ループは止められないけれども、6より大きい結果が得られたら計算自体はスキップすると言う代替です。ループ内の処理が重い場合は、高速化になります。

%dopar% if(f) を %:% when(f) %dopar% に書き換えたくなるのですが、平行処理スコープから親環境の変数は書き換え不可なので、狙った通りには動きません。f <- FALSE を assign("f", FALSE, envir=.GlobalEnv) に変えてみたりしたのですが、狙ったとおりには動きません。どうも平行処理スレッドごとに .GlobalEnv もコピーして保持しているようで、シェアードしないみたいですね。

*1: 初出と記事内容が大きく変わっていましたが、大きく勘違いしていたので訂正しています。動作検証時に %dopar% を %do% と書いてしまうポカをしていました。

*2: OpenMP並列forループの中断処理 - yohhoyの日記

プレビュー画面で特定情報を出力しないためのBloggerテンプレートの書き方

Bloggerで外部JavaScriptを読み込ませるとプレビュー画面だけ妙に時間がかかる謎現象が発生するので、その回避のために、プレビュー画面だけそこの部分の出力をしないことにしました。Bloggerの記事のプレビューを爆速にするカスタマイズ | クラウド番外地を参照して、以下のようにします。ちょっと変えてありますが、not (A and B) は、(not A) or (not B) なのでほぼ同じです。

<b:if cond='data:blog.pageType != &quot;item&quot; or data:blog.url != data:blog.homepageUrl'>

<script src='...' type='text/javascript'/>

</b:if>

何分もかかっていたのが、あっさり表示されるようになりました。とても謎です。

JavaでOpenMPもどき

気づくとマルチコア化が一般化したため、数値演算をするときにカジュアルにマルチスレッド対応にしたい状況も随分と増えてきました。計算時間が数分の一になる程度ですが、何時間もかかる計算もあります。GoやFortran 2008のように並列化容易なような言語もありますが、C言語など枯れた言語を使い続けたいと言うニーズもあります。この用途で人気なのはOpenMPです。

OpenMPはコメントアウトを利用したプリプロセッサとして実現しています。元ソースコード中のコメントにOpenMPへの命令を書いておき、OpenMPにそれに従って元ソースコードをSMP対応ソースコードに修正してから、コンパイルして実行ファイルを得るわけです。コンパイラが対応していると、コンパイルオプションを指定するだけで利用できます。

浮動小数点演算では最速のケースも多々あるJavaは数値演算でも使われるケースが多いでしょうが、OpenMPは公式にはサポートしていません。しかし、OpenMP互換プリプロセッサomp4jと言うのがあるので試しに使ってみました。チュートリアルの通り手軽に使える反面、エラーメッセージから問題究明に慣れが必要そうでした。

使うには、まずソースコードのforループなどの前に、コメントアウトした omp で始まる命令を書きます。単に命令を追加すれば良いだけではなく、平行処理するループ中ではbreakできないなど縛りがあるので、注意してください。

// omp parallel for
for (int nProductivity = 0; nProductivity < nGridProductivity; nProductivity++) {

次に、omp4j-1.2.jarを使ってコンパイル。

java -jar omp4j-1.2.jar RBC_Java.java

コンパイルができたら、実行。

java RBC_Java

ベルマン方程式によって簡単なマクロ経済モデルを解いたところ、倍以上、速くなりました。なお、言語が違うわけで当然ですが、公式と挙動が完全に一致するわけでは無いようです。例えばomp atomicでローカル変数の更新が上手くいかないので、volatileをつけたメンバー変数を更新させるなど工夫が要ると思います。

Rでテキストファイルを読み込んで分解

RでCSVファイルを上手く読み込めないときは、read.tableにあわせてテキスト・ファイルを書き換えた方が手軽な気がしますが、それが都合で出来ないときは、テキスト・ファイルとしてCSVファイルを読み込んでから、カンマで分離することもできます。

con <- file("example.csv", "r", blocking = FALSE)
lines <- readLines(con)
close(con)
lst <- strsplit(lines, ",")
row <- 1
col <- 1
sprintf("%d行%d列 %s", row, col, lst[[row]][col])
# as.numeric()で数値型に直すなど、適時型変換が必要

なお、strsplitは正規表現が使えます。
特定の列のデータをベクトルにしたい場合は以下のようにすると手軽です。

sapply(lst, "[[", 1) # 1列目を取り出す

お気軽にAR(1)の構造転換点を区間推定してみる

学術的な意義は無いのですが、時系列データの構造転換点を求めたいときがあります。記述的に転換点を仮定して推定することが多いわけですが、主観的になりやすく言い合いになることがあるからです。統計学者に殺されそうな荒業ですが、尤度を使って構造転換点を探し、さらに区間推定してみましょう。
時系列データであれば何でも良いのですが、AR(1)を考えます。

1. モデル

は期の観測値、は切片項、は推定する構造転換点、はtがTより大のときに1、他は0をとる構造転換ダミー、とは推定する係数、は誤差項です。
なお、AR(1)にダミーを入れていいのかは謎です。

2. データセット

モデルにそって、データセットを作成します。

set.seed(20150214)
n <- 100 # 観測数
y <- numeric(n) # 観測値
C <- 3 # 切片項
gamma <- 4 # 構造転換ダミー
T <- round(n/2) # 構造転換点
beta <- 0.8 # β
epsilon <- rnorm(n) # ε
y[1] <- 15 + epsilon[1]
for(i in 2:n){
  y[i] <- C + beta * y[i-1] + epsilon[i] + (i>T)*gamma
}

# 生成データを確認する
plot(y, type="l", main="データセット", xlab="t", ylab="y")

3. 構造転換点が分かっている場合

構造転換点Tが固定で分かっているのであれば、簡単に推定できます。

d <- (2:n>50)*1 # tが50以下は0、51以上は1になる
r <- lm(y[-1] ~ y[-n] + d)
summary(r)

しかし、実際にはT=50は分かりません。40〜60のどれを仮定しても、それらしい結果になります。

4. 対数尤度を最大にする点を探す

Tを変えて説明変数の数は増えたりしないので、対数尤度を最大化する点を探しましょう。対数尤度関数を書いてニュートン・ラフソン法などで探したくなりますが、構造転換点は不連続なため上手く推定できません。ループして力技で最大値を求める方が確実です。

logLik <- numeric(n) # 対数尤度の保存用ベクトル
logLik[1] <- logLik[n] <- NA
ml_scp <- 0 # 尤度が最大の点
r <- lm(y[-1] ~ y[-n]) # 尤度が最大の推定結果
for(scp in 2:(n-1)){
  d <- (2:n>scp)*1
  tmp <- lm(y[-1] ~ y[-n] + d)
  logLik[scp] <- logLik(tmp) # 対数尤度は蓄えておく
  if(logLik(r) < logLik(tmp)){
    r <- tmp
    ml_scp <- scp
  }
}
summary(r)
sprintf("最尤推定された構造転換点: %d", ml_scp)

5. 構造転換点の区間推定を行なう

構造転換点が分かれば十分なときも多いわけですが、信頼性を疑われるので区間推定が行なえるように標準誤差を求めたいです。しかし、対数尤度の集合は離散データなので、このままでは最尤推定値の標準誤差は求まりません。

そこでスプラインで補間して、連続な対数尤度関数モドキを作ってしまいましょう。

t <- 2:(n-1)
sp <- smooth.spline(t, logLik[t])
f <- function(p){
  predict(sp, p)$y
}

プロットするとそこそこの精度で近似できていることが分かります。全体としては丸くなるので、区間推定の範囲は広くなりそうですが、狭くなるよりは良いでしょう。

plot(logLik, type="p", main="対数尤度", xlab="t", ylab="y")
lines(f(1:n))

作った対数尤度関数モドキから、構造転換点を再推定します。

r_ml <- nlm(function(p){
  -f(p)
}, ml_scp, hessian=TRUE)
SEs <- sqrt(diag(solve(r_ml$hessian))) # 標準誤差

これで構造転換点の標準誤差が出ました。区間推定もしてみましょう。

interval <- function(beta,se,range,nof){
  a <- 1 - range/2
  sprintf("%.3f(%d%%信頼区間%.3f〜%.3f)",beta,as.integer((1-range)*100),beta-se*qt(a, nof),beta+se*qt(a, nof))
}
interval(r_ml$estimate[1], SEs[1], 0.05, summary(r)$df[2])

厳密な方法ではありませんが、tは49から51ぐらいを見とけば良いとなります。