読者です 読者をやめる 読者になる 読者になる

餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

最適化レベルで結果が変わってしまうC言語のダメなコード

C/C++

前世紀に書いたgccで-O1と-O3で実行結果が異なるコードです。インライン展開が作用しました。今のバージョンのgccでどうなるかは確認していません。なお、警告は出ると思いますし、処理系によってはコンパイルが通らない事もあると思います。

#include<stdio.h>

int *ptr;

void test(void){
  int a;
  a = 10;
  ptr = &a;
}

void test2(void){
  int b;
  b = 23;
}

main(){
  test();
  printf("Result1: %d\n",*ptr);
  test();
  test2();
  printf("Result2: %d\n",*ptr);
}

Rの関数をラッピング

R言語

Rは対話型インターフェイスなので、長々と関数に引数を指定するのは手間隙なので、グローバル・オプションを設定したくなります。しかし、発展途上のライブラリが多いため、グローバル・オプションをつけると上手く動作しなくなったり、そもそもグローバル・オプションを参照しないライブラリがあったりします。こういう都合の悪いときに、やや強引にデフォルト値を設定する細工を考えてみました。

R言語では既に宣言されている関数や変数を再定義する事ができ、さらに再定義した関数内からもともと宣言されている関数を呼び出す事ができます。これを使って、read.tableのfileEncodingにデフォルト値を設定してみましょう。

read.table <- function(...){
  print("ラッパーを通しています(理解したらここは消してください)")
  # 可変長引数をリストに展開
  args <- list(...)
  # fileEncodingが存在しなければ、UTF-8をセット
  if(!("fileEncoding" %in% names(args))){
    args <- append(args, list("fileEncoding"="UTF-8"))
  }
  # ホンモノのread.tableを呼び出す
  # package:utils はホンモノがあるEnvironment
  do.call(get("read.table", env=as.environment("package:utils")), args)
}

これを .Rprofile *1に書いておけば、options(encoding="UTF-8")なしにデフォルトでUTF-8指定ができます。write.tableでも、他の関数でも同様にできます。なお、ラッピングする関数の環境(Environment)が分からない場合は、find("関数名")としてください。
しかし、こういうコーディング、MS-DOS時代の割り込み処理を思い出しますね!

*1: Windowsの場合は、%USERPROFILE%\Documents\.Rprofile になります。

Rの文字コードのデフォルト値の設定

Windows版のRはテキストファイルの入出力にCP932と言うかShiftJISを使うわけですが、LinuxMacintosh OS XWindowsを併用して使っている人は、UTF-8をデフォルトにしたいかも知れません。read.tableのfileEncodingオプションをいちいち指定するのは面倒です。
解決方法は単純で、起動時に読み込む .Rprofile *1に、options(encoding="UTF-8") を書いておけば、UTF-8をデフォルトにします。設定値を確認したい場合は options()$encoding を表示してください。
なお、オブジェクトをそのまま保存するsave()とload()はコード変換が行なわれないので、使わない方が無難です。読み書き速度が問題になる場合は、readr パッケージを使うと改善できます。

*1: Windowsの場合は、%USERPROFILE%\Documents\.Rprofile になります

マクロ経済学の動的計画法説明用コードのMatlabからRへの移植

R言語

f:id:uncorrelated:20161224234007p:plain

一橋大学の阿部教授が院生向けのレクチャー・ノートをアップロードしていたのを見つけ、拝読して勉強させて頂いていたのですが、説明用のコード(P.16--20)がMatlabだったのでRに移植してみました。プロプライエタリMatlab嫌いな人に役立つかも知れないので、公開します。

単に移植するだけではなく、ループ周りを中心に無駄を省いています。理由は、行列にしておいた方が変数の中身を説明しやすいのと、コード全体の見通しが良くなるのと、速度面で有利になるからです。資本量の区切り幅inK=0.4ぐらいにして、行列の中身を表示させると、何を演算しているのか理解しやすいと思います。

alpha <- 0.25 # production parameter
beta <- 0.9 # subjective discount factor
gamma <- -2 # preference parameter
delta <- 1 # 1 -depreciation rate
minK <- 0.2 # minimum value of the capital grid
maxK <- 1.8 # maximum value of the capital grid
inK <- 0.001 # size of capital grid increments
nK <- ceiling((maxK-minK)/inK + 1) # number of grid points

# 横(列)の当期資本と、縦(行)の来期資本で定まる効用を表す行列
U <- matrix(rep(-Inf, nK^2), nK, nK)

K <- t(replicate(nK, ((1:nK)-1)*inK + minK)) # 当期資本。nK×nKに拡張。当期資本は列ごとに変えているのに注意。
K_prime <- replicate(nK, ((1:nK)-1)*inK + minK) # 来期資本k_{t+1}。nK×nKに拡張。来期資本は行ごとに変えているのに注意。
I <- K_prime - delta*K # 投資量。横(列)の当期資本と、縦(行)の来期資本で定まる。
C <- ((1-beta)/(alpha*beta))*(K^(alpha)) - I # (79)式を(80)式に代入した式から消費量を計算
U[C>0] <- ((C^(gamma+1))/(gamma+1))[C>0] # 効用を計算して行列に格納。消費がゼロ以下の場合は、初期値-Infを保存している

#
# whileループのための変数を初期化
#
V <- rep(0, nK) # value functionの代わりになるベクトル
Decis <- rep(0, nK) # Vに対応する来期資本を表すグリッド位置
iter <- 0 # ループ回数
metric <- 10 # 更新前後のvalue functionの差分の最大要素の値が入る。初期値は終了条件より大きければ、何でも良い。

while(metric > 0.001){
  # (82)式の最大化問題に対応する来期資本を探す
  # beta*Vは1列分のベクターだが、自動的にnK列分に複製拡大して足されるのに注意
  r <- U + beta*V

  # それぞれの列の最大値を選んで、value functionとする。
  # 当期資本に対応する列の中の最大値は、最大化問題の解となる来期資本
  tmpV <- apply(r, 2, max) # value functionの値
  tmpDecis <- apply(r, 2, which.max) # value functionに対応する来期資本を表すグリッド位置

  metric <- max(abs(V - tmpV)) # 終了条件と比較する更新によって生じた差分を計算

  V <- tmpV;
  Decis <- tmpDecis;
  iter <- iter + 1 # ループ回数をカウント

# 計算状況を表示
  vfor8 <- V[1] # value of V for k=kmin=0.2
  ufor8 <- (Decis[1] - 1)*inK + minK # value of control for k=0.2
  print(sprintf("ループ回数:%d, VF変化量:%.4f, vfor8:%.4f, ufor8:%.4f", iter, metric, vfor8, ufor8))
}

# policy functionと対応する各期の効用を計算
policy <- (Decis - 1)*inK + minK # policy function
U_op <- rep(-Inf, nK) # Uity under the optimal policy
K <- ((1:nK)-1)*inK + minK
I <- policy - delta*K;
C <- (1-beta)/(alpha*beta)*(K^alpha) - I
U_op[C>0] <- ((C^(gamma+1))/(gamma+1))[C>0]

# 各期の効用の割引価値を計算
betam <- beta*rep(1, nK)
value <- U_op/(rep(1, nK)-betam)

# 結果をプロット
x_at <- seq(1, length(value), length.out=17)
x_labels <- sprintf("%.2f", seq(minK, maxK, length.out=length(x_at)))
y_unit <- 5
ylim <- c(floor(min(value)/y_unit)*y_unit, ceiling(max(value)/y_unit)*y_unit)
y_at <- sprintf("%.2f", seq(ylim[1], ylim[2], length.out=5))
plot(value, main="DETERMINISTIC GROWTH MODEL: VALUE FUNCTION", xlab="Capital", ylab="Present Value", axes=FALSE, ylim=ylim)
axis(1, at=x_at, labels=x_labels, las=2)
axis(2, at=y_at)

線形近似の方は省略しました。

Rで機械学習(SVM)

R言語

人間の直観的な認識力を手軽に模倣するフレームワークである機械学習がもてはやされるようになってもう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"))

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

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

f:id:uncorrelated:20160614202854p:plain

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

5. まとめ

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

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

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

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

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

R言語

f:id:uncorrelated:20160526172931p:plain

@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. 結果

f:id:uncorrelated:20160526170750p:plain

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

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

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

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

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

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

R言語

確認してみたのですが、試行錯誤してみた結果、結論は「ありません」でした*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の日記