Windows版のRはテキストファイルの入出力にCP932と言うかShiftJISを使うわけですが、LinuxやMacintosh OS XとWindowsを併用して使っている人は、UTF-8をデフォルトにしたいかも知れません。read.tableのfileEncodingオプションをいちいち指定するのは面倒です。
解決方法は単純で、起動時に読み込む .Rprofile *1に、options(encoding="UTF-8") を書いておけば、UTF-8をデフォルトにします。設定値を確認したい場合は options()$encoding を表示してください。
なお、オブジェクトをそのまま保存するsaveRDSとreadRDSはコード変換が行なわれないので、使わない方が無難です。読み書き速度が問題になる場合は、readr パッケージを使うと改善できます。
マクロ経済学の動的計画法説明用コードのMatlabからRへの移植
一橋大学の阿部教授が院生向けのレクチャー・ノートをアップロードしていたのを見つけ、拝読して勉強させて頂いていたのですが、説明用のコード（pp.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) # Utility 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）
人間の直観的な認識力を手軽に模倣するフレームワークである機械学習がもてはやされるようになってもう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ぐらいの女性データが男性に分類されてしまっていますが、全体としては悪くない当てはまりのようです。
結婚市場からの離脱を考慮した夫婦の社会的階層格差
@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)
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% と書いてしまうポカをしていました。
プレビュー画面で特定情報を出力しないためのBloggerテンプレートの書き方
Bloggerで外部JavaScriptを読み込ませるとプレビュー画面だけ妙に時間がかかる謎現象が発生するので、その回避のために、プレビュー画面だけそこの部分の出力をしないことにしました。Bloggerの記事のプレビューを爆速にするカスタマイズ | クラウド番外地を参照して、以下のようにします。ちょっと変えてありますが、not (A and B) は、(not A) or (not B) なのでほぼ同じです。
<b:if cond='data:blog.pageType != "item" 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を使ってコンパイル。
コンパイルができたら、実行。
ベルマン方程式によって簡単なマクロ経済モデルを解いたところ、倍以上、速くなりました。なお、言語が違うわけで当然ですが、公式と挙動が完全に一致するわけでは無いようです。例えばomp atomicでローカル変数の更新が上手くいかないので、volatileをつけたメンバー変数を更新させるなど工夫が要ると思います。
