「熊とワルツを」のリスクシミュレーションをRで実装した

最近「熊とワルツを」というプロジェクトにおけるリスク管理についての本を読みました．本の中にでてくるリスク管理シミュレーションをR実装しました，という話です．

熊とワルツを - リスクを愉しむプロジェクト管理

ネタ元のトム・デマルコさんの本は，軽妙な語り口でとても読みやすいので，プロジェクトベースでお仕事する人はみんな一度は目を通してみるといいんじゃないかなと思います*1．偉い人の思いつきで，ケツの決まったプロジェクトを押し付けられそうになったときに，そんなんじゃ崩壊するだろって言い返すためにも*2．そんなわけで，シミュレーションの説明をする前段として，本の中で述べられているリスクについて軽く触れておきます．

5つのコアリスク

この本ではリスク管理について様々な観点から述べられているのですが，その中で主なリスク要素として以下の5つが挙げられています．

1. スケジュールの欠陥
2. 要求の増大
3. 人員の離脱
4. 仕様の崩壊
5. 生産性の低迷

これらのリスク要因の中でも，4番目に挙げられた「仕様の崩壊」については，これが起こるとプロジェクト自体が崩壊して失敗するものとされています．それ以外の4つについては，発生によってスケジュール自体が5%〜20%のような幅のある見積もりで遅延するものと考えられます．

シミュレーションによるスケジュール遅延リスクの見積もり

コアリスクによるスケジュールの遅延について，本の第11-14章あたりでシミュレーションによる見積もり法が紹介されています．このツール，実際に著者のサイトからダウンロードすることができます．しかしそのツールはExcelで作られており*3，正直ちゃんと使うには結構つらみがある仕様になっています．最終更新日も2004年で，今後更新されるようにも思えないので，Rで作り直してみました，というのが今回の趣旨です．

ツールの概要

ツールの中で実際に行われているのは，スケジュール遅延予測分布に基づいた，モンテカルロシミュレーションです．平たくいうと，スケジュールを実施する人が，いくつかのリスク要素についてあらかじめスケジュール遅延度合いを以下の3点で見積もります．

1. スケジュールの最小遅延割合
2. スケジュールの最頻遅延割合
3. スケジュールの最大遅延割合

これを5つのコアリスクについて見積もります．ただし4番目の「仕様の崩壊」については，これが生じるとプロジェクト自体が失敗してしまうクリティカルなリスクということで，別途以下のように見積もります．

• 仕様崩壊が発生する確率
• 仕様崩壊が生じたときに，プロジェクトが失敗するか否か
• プロジェクトが失敗しない際に遅延する月数
• プロジェクトが失敗しない際に遅延する割合

シミュレーションでは，プロジェクトが失敗しないパターンも見積もることが可能で，その場合は遅延月数もしくは遅延割合のどちらかで，スケジュール遅延を表します．この数字は一般にとても大きくなるものと考えられます．

上記の値を設定した上で，乱数を用いて遅延度合いのサンプルを取得し，実際にどのくらいスケジュールが遅延するか（もしくは失敗するか）を計算します．これを繰り返し実施して，スケジュール遅延度合いの分布を作成します．これによってだいたいスケジュールがどのくらい遅延するかの目安を得ることができます．

Rでの実装

ということでコードをみていきます．フルバージョンはgithubリポジトリを参照してください．コードは設定ファイルと実行ファイルの2つに分かれています．設定ファイルは以下のようになります．

# trial iteration count
N_TRIAL = 3000

# project period
START_DATE = as.Date("2015/1/1")
END_DATE   = as.Date("2016/9/30")
SPAN       = END_DATE-START_DATE

# fatal risk: SPECFLAW(failure to reach consensus)
FR = list("prob"=0.15, "fatal"=TRUE, "month"=8, "rate"=0)

# other risk factors
RF = list(1:10, list())
## RF01: SCHEDFLAW(error in original sizing) 
RF[[1]]  = list("min"=0.0, "mode"=0.2, "max"=0.5)
## RF02: TURNOVER(effet of employee turnover) 
RF[[2]]  = list("min"=0.0, "mode"=0.05, "max"=0.1)
## RF03: INFLATION(requirement function growth) 
RF[[3]]  = list("min"=0.0, "mode"=0.1, "max"=0.2)
## RF04: PRODUCTIVITY(effect of productivity variance) 
RF[[4]]  = list("min"=0.0, "mode"=0.05, "max"=0.2)
## RF05: OTHER RISK 
RF[[5]]  = list("min"=0.0, "mode"=0.0, "max"=0.0)
## RF06: OTHER RISK 
RF[[6]]  = list("min"=0.0, "mode"=0.0, "max"=0.0)
## RF07: OTHER RISK 
RF[[7]]  = list("min"=0.0, "mode"=0.0, "max"=0.0)
## RF08: OTHER RISK 
RF[[8]]  = list("min"=0.0, "mode"=0.0, "max"=0.0)
## RF09: OTHER RISK 
RF[[9]]  = list("min"=0.0, "mode"=0.0, "max"=0.0)
## RF10: OTHER RISK 
RF[[10]] = list("min"=0.0, "mode"=0.0, "max"=0.0)

4番目の「仕様の崩壊」については，別計算が必要なので"fatal risk"として別途切り出しています．項目は順にリスクの発生確率，発生時にプロジェクトが崩壊するか否か，崩壊しない場合の遅延月数，遅延割合です．「仕様崩壊」がないとする場合には，prob=0としてください．それ以外の4つ+予備の5つのリスク要素について，最小値，最頻値，最大値を割合で入力します．その要素の影響が特にない場合には，3つとも0に設定してください．

それ以外に，プロジェクトの開始日時と終了日時を設定可能です．実際のプロジェクト開始日と終了日を適宜変更して入れてください．またシミュレーションの実行回数は，とりあえず3000回としてありますが，任意に変更可能です*4．本体のriskology.Rを実行すると，結果として以下のように結果が得られます．また "img/histgram_of_delays.png" にヒストグラムが画像で保存されます．ヒストグラムを見るとわかるように，今回の設定だとプロジェクトは最低100日以上，平均して250日程度遅延することがみてとれます．

> sprintf("iteration count: %d", N_TRIAL)
[1] "iteration count: 3000"
> sprintf("project cancel count with fatal risk: %d", cancelled)
[1] "project cancel count with fatal risk: 463"
> sprintf("cancel rate: %02.1f%%", cancelled/N_TRIAL*100)
[1] "cancel rate: 15.4%"
> print("Summary of delay days:")
[1] "Summary of delay days:"
> summary(delays)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  108.0   217.0   253.0   254.2   291.0   426.0 
> 

そんなわけで，ちゃんとスケジュール遅延の可能性を分布で把握して，炎上しないように楽しくプロジェクトをこなしていきましょう*5．

最後にそしてシミュレーション本体はこちらです．

source("conf.R")

# risk distribution sampler
sampler = function(rf_, total=10) {
  range = rf_[["max"]]-rf_[["min"]]
  alpha = total*(rf_[["mode"]]-rf_[["min"]])/range
  beta = total*(rf_[["max"]]-rf_[["mode"]])/range
  sample = rbeta(1, alpha, beta)
  sample*range
}

# calculate fatal risk
fatal_risk = function(fr) {
  if (fr[["prob"]] < runif(1)) {
    return(0)
  } else {
    if (fr[["fatal"]]) {
      return(-1)
    } else {
      month = fr[["month"]]*30
      rate = (END_DATE-START_DATE)*fr[["rate"]]
      return(max(month, rate))
    }
  }
}

# calculate delay factor
delay_factor = function(RF) {
  df = 0
  for (i in 1:length(RF)) {
    df = df+sampler(RF[[i]])
  }
  df
}

# calcurate delay
cancelled = 0
delays = c()
for (i in 1:N_TRIAL) {
  delay = 0
  # check fatal risk
  fatal = fatal_risk(FR)
  if (fatal == -1) {
    cancelled = cancelled+1
    next
  } else {
    delay = delay+fatal
  }
  # add other risk factors
  for (i in 1:length(RF)) {
    delay = delay+(END_DATE-START_DATE)*sampler(RF[[i]])
  }
  delays = c(delays, as.integer(delay))
}

# draw histgram
png("img/histgram_of_delays.png", width=800, height=600)
hist(delays,
     breaks=max(10, as.integer(N_TRIAL/30)),
     main="Histgram of simulated delay days",
     xlab="Delay days",
     ylab="Frequency")
dev.off()

# print result
sprintf("iteration count: %d", N_TRIAL)
sprintf("project cancel count with fatal risk: %d", cancelled)
sprintf("cancel rate: %02.1f%%", cancelled/N_TRIAL*100)
print("Summary of delay days:")
summary(delays)

*1:デマルコさんの本だと，「ピープルウェア」も名著だと思います．

*2:世の中は世知辛いので，わりとこういう炎上プロジェクトの2つや3つ，わりとそこらへんに転がっているものです．

*3:みんなが大好きなExcelです．とはいえコードベースでの管理ができないので，バージョン管理システムとの相性の悪さは天下一品です．

*4:もとのExcelだと500回になっていました．

*5:まぁ偉い人がやれといったらどうにもならないことも，ままありますけどね．そんときはそんときということで．

dplyrを使ったdata.frameの前処理を関数化する

表題の通り，dplyr使って前処理する際に，それを関数化する方法のメモ．ユースケースとしては，ちょっとだけ条件変えてデータフレーム自体を何度か出し直すってとき用．

関数

まるっとフィルタ用のconditionと，select用のvariablesを引数で渡します*1．ポイントは，filterではなくfilter_とアンスコが付いていること．これは，この処理がNon-Standard Evaluation*2で評価してはいけないからです．

preprocess = function(df, condition, variables) {
  df %>%
    filter_(condition) %>%
    select_(variables)
}

呼び出し側

呼び出し側は，条件文をquote()に入れて渡します．selectで使う変数の方は，c()で包んでリストにしておく必要あり．

> res = preprocess(iris,
+                  quote(Species=='virginica'),
+                  quote(c(Sepal.Length, Sepal.Width)))
> nrow(res)
[1] 50
> head(res)
  Sepal.Length Sepal.Width
1          6.3         3.3
2          5.8         2.7
3          7.1         3.0
4          6.3         2.9
5          6.5         3.0
6          7.6         3.0

これで前処理スッキリ，再利用ブラボー．

おまけ

ちょっと違うアプローチでいうと，以下のようなのもあります*3．

直に指定

preprocess = function(x, column, fn) {
  fn(x[,column])
}

> preprocess(iris, 'Sepal.Length', max)
[1] 7.9

subsetを使う

preprocess = function(df, condition) {
  e = substitute(condition)
  r = eval(e, df, parent.frame())
  subset(df, r)
}

> res = preprocess(iris, Sepal.Length > 7.0)
> nrow(res)
[1] 12
> head(res)
    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
103          7.1         3.0          5.9         2.1 virginica
106          7.6         3.0          6.6         2.1 virginica
108          7.3         2.9          6.3         1.8 virginica
110          7.2         3.6          6.1         2.5 virginica
118          7.7         3.8          6.7         2.2 virginica
119          7.7         2.6          6.9         2.3 virginica
> ||<

*1:参考: Passing strings as arguments in dplyr verbs - Stack Overflow

*2:NSEについては R - NSEとは何か - Qiitaとか， dplyrのなんたら_eachを効率的に使う - 東京で尻を洗うとか， NSE（Non-Standard Evaluation）について - 東京で尻を洗うの記述が詳しいです．大元は Non-standard evaluationなので，こちらを読むのが良さげです．

*3:参考: r - Pass a data.frame column name to a function - Stack Overflow

BW変換（Burrows Wheeler Transform）を書いた

「高速文字列解析の世界」を読んでて，BW変換のところがよくわからなかったので，実際に書いてみました．

高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学)

概要

BW変換ってなんぞやを一言でいうと，「変換すると，似たような記号がたくさん並ぶようになる，文字列の可逆変換」です．具体例をみせると，以下はWikipediaのBW変換の最初のパラグラフを，実際にBW変換してみた結果です．

元の文章*1

Compression techniques work by finding repeated patterns in the data and encoding the duplications more compactly. The Burrows Wheeler transform (BWT, also called block-sorting compression) rearranges a character string into runs of similar characters. This is useful for compression, since it tends to be easy to compress a string that has runs of repeated characters by techniques such as move-to-front transform and run-length encoding. More importantly, the transformation is reversible, without needing to store any additional data. The BWT is thus a free method of improving the efficiency of text compression algorithms, costing only some extra computation.

変換後の文章*2

mee....ssssyn,amehodksoadreotgra,eleeedheylatefsgnTseesstdssygd) fgsdsof,yorasyntfgngg, ssygrtes,s. nTsnyeonkegaysn ( $ WW BB t r trrrp n c rrrrt lhhhe hehddeemtcp i inu eeionn aaaaneenoen anoeodn hhrrhhbhmrcehlvr ppttltternh se itlrrtltttvguurrrrrmrt ooooef sss -ennnnnnnnnnlnctccc tTTttttWTtcctttslfcms sftrtdrdvdd stssstttnnh rwdrcauibl-eapabanttrrro i oioooooiohoooioouooieeeaaie iiiiiiiiiaehhorouuaaaoaiatsttttlhcc scccCccciiiiiii irftmMgwfffpschmrrm eeummmmmmmmiiaoeetaaartttooo f ppppp ttoooooefpraueeeoo weusrninanniaeidiemrunnnr ssssl -eeeeeo axauinaaru aa ctcc g iei rsaaianc n -sx ssasqqfrrrdBh opoeo oeenbclbsll

ということで，変換すると同じような記号がダーっと並んでるのが見て取れると思います．これの何が嬉しいかというと，同じような文字が連続していると，文字列を圧縮するときに圧縮効率が良くなります*3．

実装

ということで，BWTをするのと，元に戻すのを実際に実装してみました．とはいえ理解のための実装なので，パフォーマンスについては度外視してます．あと，終端文字として$を使っているので，文章中に$があるとうまく変換できません．

これを実際に実行してみると，以下のように可逆に変換されました．

python bwt.py abracadabra
ard$rcaaaabb
abracadabra

本の簡潔データ構造とかちゃんと読んでから，もう一回パフォーマンスチューニングしますかね...

参考

詳しい説明やパフォーマンス改善については，以下の2つの記事を読んでもらえるとよいと思います．特にLFマッピングの部分は，コードだけだと意味がわからない筈...

• Burrows Wheeler Transform と Suffix Array - naoyaのはてなダイアリー
• Burrows Wheeler TransformとLF mapping | Preferred Research

*1:実はハイフンがうまく処理できなかったので，「Burrows-Wheeler transform」を「Burrows Wheeler transform」にこっそりと変えてあります．

*2:blogのフォーマットにキレイに収まるように，半角スペースを3個くらい付け加えてあります．

*3:wikipediaにあるように実際にbzip2の圧縮アルゴリズムの中で使われています．

データビジネス・分析・開発に関して2014年に読んだ本

年末なのでぼちぼち今年の振り返りをします．ちなみに去年のはこちら．

データブジネス，データ分析，ソフトウェア開発の3カテゴリに分けて，それぞれについて上から読んでよかった順に並んでいます．

データビジネス

"超"分析の教科書

“超"分析の教科書 (日経BPムック)

事例集として，非常によくまとまっていました．幅広い業種で典型的に使われるような手法とか問題とかがコンパクトにまとまっていて，実務でデータ分析をしている人ならみておいて損はないのではないでしょうか．内容は割と平易に書かれているので，データ分析専業じゃない人が読んでも割と読みやすく面白いと思います．

アルゴリズムが世界を支配する

アルゴリズムが世界を支配する (角川ＥＰＵＢ選書)

コンピュータ・アルゴリズムがどういう場面で使われて，それが世の中をどう変化させているのかを，具体的な事例とともに語っている本で，単純にとても面白いです．金融業のアルゴリズムトレードの話がメインではありますが，それ以外にも政府系の機関や医療での使われ方が書かれていて，なるほどーという感じです．もっとフラッシュトレード寄りになったらマイケル・ルイスのフラッシュ・ボーイズなんかがありますね*1．

フラッシュ・ボーイズ 10億分の1秒の男たち

夜の経済学

夜の経済学

社会調査手法を使って，風俗・水商売の業界を探っている本です．扱っているネタはキワモノですが，手法や中身は非常にまともだし，考察も示唆に富んでいてとても面白いです．社会の裏側とかが好きなタイプの人は楽しめる本かと思います．

サイレント・ニーズ

サイレント・ニーズ――ありふれた日常に潜む巨大なビジネスチャンスを探る

きちんとした定性的な調査に基づいたサービスデザインを考える際に，非常に参考になる本だと思います．具体的な事例ベースのところが面白いですね．

快感回路

快感回路---なぜ気持ちいいのか なぜやめられないのか (河出文庫)

学術書を一般向けに書き改めたタイプの本で，まだ結構硬い感じが残っています．それさえ苦にしなければ，非常にエキサイティングな中身の本です．「ハマる」とはどういうことか，何をすれば脳が快感を感じるのか，というのを（キワモノではない）脳科学の研究成果と合わせて論じていて，サービス設計にどう活かすかというのを考えるのの一助になるかなーと思ったりします．

振り込め犯罪結社

振り込め犯罪結社 200億円詐欺市場に生きる人々

なんで人は振り込め詐欺にかかるのか，また詐欺を実行している人たちの組織はどうやって回っているのか，という点について非常に具体的にまとめられていて，組織をまとめるという意味で非常に参考になります．普通に考えて，犯罪組織をまとめるのって，普通の組織をまとめるのよりも数段大変だと思うんですよね．その意味でも彼らのやり方というのはとても勉強になるわけです．

インタフェースデザインの心理学

インタフェースデザインの心理学 ―ウェブやアプリに新たな視点をもたらす100の指針

基本的なUI/UXのセオリーを学ぶという意味で，とても良い本だと思います．

ビジネスモデルの教科書

ビジネスモデルの教科書: 経営戦略を見る目と考える力を養う

ビジネスモデルを類型化してまとめた本．カタログ的な意味で理解のフレームワークとして使うのが良いかなと思います．

A/B Testing: The Most Powerful Way to Turn Clicks Into Customers

A/B Testing: The Most Powerful Way to Turn Clicks Into Customers

ABテストについてまとめられた本って，実は少ないんですよね．英語なのがアレですがこの本は割とよくまとまっている印象．

最強のWebコミュニケーション・シナリオ

最強のWEBコミュニケーションシナリオ

ペルソナ周りを調べるときに参考した本です．ペルソナについてちゃんと扱っている本ってあんまりなくて，探すのに苦労しました．

Lean UX

Lean UX ―リーン思考によるユーザエクスペリエンス・デザイン (THE LEAN SERIES)

どうやってUXを改善していくかというお話．

モバイルフロンティア

モバイルフロンティア よりよいモバイルUXを生み出すためのデザインガイド

モバイルUXについて書かれた本ですが，なんかピンボケしててあんまり頭に入らなかった...

「ヒットする」「タッチパネル」「レベルアップ」のゲームデザイン

「ヒットする」のゲームデザイン ―ユーザーモデルによるマーケット主導型デザイン

「タッチパネル」のゲームデザイン ―アプリやゲームをおもしろくするテクニック

「レベルアップ」のゲームデザイン ―実戦で使えるゲーム作りのテクニック

ゲームデザイン本3冊ですが，全体的にピンと来なかった印象...

データ分析

入門リスク分析

入門リスク分析―基礎から実践

統計の一つの役割というのは，信頼区間や予測区間を出すことによって，取ろうとしている施策がどれくらいのリスクを持つか，というのを定量的に明らかにすることにあると思います．実務だと，割と仮説検定して有意水準を下回っていればOKみたいな感じで道具として使われがちですが，単に検定するだけではなく，分布をちゃんと可視化してリスクを定量化する，ということの意味はもっと認識されて良いと思います．

その一方でそうしたリスク分析に関する本は非常に少なくて，本書は非常にわかりやすく，かつ丁寧にまとめられていてとても良いです．本の中頃には，主要な確率分布のまとめもあって，地味にこれがリファレンス的に結構便利だったりします．惜しむらくは，値段が高いのと，すでに廃刊になっていて中古でしか手に入らない点です... と思ってたらいつの間にかKindle版が出てますね．とはいえ¥14,910ととってもお高いですが...

予測にいかす統計モデリングの基本

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)

状態空間モデル〜パーティクルフィルタあたりの入門書として，非常にわかりやすく良い本だと思います．ただ例のごとく，高校卒業程度の数学力があればいいとか書いておきながら，普通に数式展開ガチガチなので，そういうのが嫌いな人には辛いです．逆にそのあたりが追える人であれば，本当に理解しやすい良書です．

応用のための確率論入門

応用のための 確率論入門

確率論，測度論，ルベーグ積分の周りについてとても丁寧にまとめられていて，わかりやすいです．このあたりをわかりやすく説明している本はとても少ないので，貴重です．測度論周りは，大概ガチの専門書になっていて，実務サイドで統計を扱っている人が理解できる範疇のものがほとんどないのですね...

数学は言葉

数学は言葉―math stories

数式や数学記号も，広い意味では人に意味を伝えるための言語と捉えることができます．ただ扱う範囲が限定されていて，かつ厳密性が求められるというだけの話です．それを中学生でも十分に理解できる形で書いていて，本当にわかりやすいです．この本を読めばε-δ論法もバッチリ理解できます．

これなら分かる最適化数学

これなら分かる最適化数学―基礎原理から計算手法まで

細かい言葉の使い方や数式展開をすぐ忘れちゃうんで，この本を定期的に読み直すんですが，本当にわかりやすいです．

プログラマのための線形代数

プログラミングのための線形代数

線形代数のアルゴリズムを，実際にプログラムに落とし込んでいるので理解しやすくなっています．文章自体もわかりやすく加工という意図が見て取れてよいです．ただレベル感がアンバランスで，全くの初学者だと厳しいかも．先に，意味がわかる線形代数あたりを読んでからのほうが良い気がします．

まずはこの一冊から 意味がわかる線形代数 (BERET SCIENCE)

科学と証拠

科学と証拠-統計の哲学 入門-

統計哲学の入門書，といったところでしょうか．統計ってそもそもが，いろいろとややこしい考え方を前提としている部分があって，そのあたりの問題点や考え方の違いなんかを比較してまとめている本です．まぁ好きな人には良いかなぁという感じの本だと思います．

手を動かしながら学ぶ ビジネスに活かすデータマイニング

手を動かしながら学ぶ ビジネスに活かすデータマイニング

そもそも私自身は対象読者ではないと思うんですが，初級者を抜けた人が読むのに良い本だなぁという感想です．そういう意味でターゲットがニッチな感じはします．

ビジネス活用で学ぶデータサイエンス入門

ビジネス活用事例で学ぶ データサイエンス入門

ストーリーとしては良いと思うんですが，Rコードについての説明が非常に少なく，本自体が分析の初心者向けなのに，R本としては全く期待できないのが辛い本だと思います．

エンジニアのための データ可視化[実践]入門

エンジニアのための データ可視化[実践]入門 ~D3.jsによるWebの可視化 (Software Design plus)

D3.jsの本ではなく，あんちべさんの可視化お作法講義がまとまっている本です．企画さんに読んでもらいたい一冊．

データ分析支えるスマホゲーム開発

データ分析が支えるスマホゲーム開発 ~ユーザー動向から見えてくるアプリケーションの姿~

具体的なシステム構築の事例が割と参考になります．こういうのを本にして外に出していくスタンスは非常に良いなーと思います．

実践機械学習システム

実践 機械学習システム

とりあえず手を動かしましょう，というスタンスの本で，最近やったマンションポエムの自然言語処理あたりでも結構参考にしました．ただ中身が結構散漫で，興味のあること，やってみたいことがある場合に，この本の内容が参考になるのであればその章だけ読む，というスタンスが良い本かなーと思います．

データサイエンス講義

データサイエンス講義

サーっと目を通しただけですが，結構ずっしり書いてある本な気がします．下の本と割と被っている印象です．

戦略的データサイエンス入門

戦略的データサイエンス入門 ―ビジネスに活かすコンセプトとテクニック

上の本と結構被っている印象です．と思ったら訳者が結構被ってますね...

日本語入力を支える技術

日本語入力を支える技術　?変わり続けるコンピュータと言葉の世界 (WEB+DB PRESS plus)

安定の技術評論社のシリーズですが，この本は結構レベル感がアンバランスで，割と読者に求める要件がシビアというか，自然言語処理も機械学習もプログラミングもそこそこの事前知識がないと，読むのが結構辛いんじゃないかなぁと思います．上記の基準をクリアできていたら，それでもわかりやすく書かれた良書だと思うんですが...

現代思想2014年6月号

現代思想 2014年6月号 特集=ポスト・ビッグデータと統計学の時代

著者による当たりと外れの差が大きい特集号でした．ザーッと一読してみると，いろいろなスタンスや考え方があって面白いなぁとは思います．

ベイズ統計の理論と方法

ベイズ統計の理論と方法

買ったはいいものの，難しくて積ん読化してます．

ソフトウェア開発

ZooKeeperによる分散システム管理

ZooKeeperによる分散システム管理

お仕事で分散システムをよく使うわけですが，だいたい裏でZooKeeperが動いていて，実は中身を知らなかったわけです．分散処理制御の基本概念がとてもコンパクトにまとまっていて良いです．

Hadoopオペレーション

Hadoopオペレーション ―システム運用管理ガイド

読み直してみて，非常に細かくいろんなエッセンスが詰まってるなぁと，改めて思いました．

サーバ／インフラエンジニア養成読本 ログ収集〜可視化編

サーバ／インフラエンジニア養成読本 ログ収集〜可視化編 [現場主導のデータ分析環境を構築！] Software Design plus

fluentd/Kibana/ElasticSearchなど最近はやりのログ・可視化系の技術についてまとめられていてお得だと思います．

体系的に学ぶ 安全なWebアプリケーションの作り方

体系的に学ぶ 安全なWebアプリケーションの作り方　脆弱性が生まれる原理と対策の実践

定番の徳丸本です．今更ですが読みました．

教養としてのコンピュータ・サイエンス

教養としてのコンピュータ・サイエンス (Information Science&Engineering)

私自身，計算機科学を大学で専攻していたわけではないので，結構基礎的な事柄に抜けが多いのですね．この本は入門書なので当たり前かもしれませんが，かなり平易に書かれていて高校生でも十分に理解できる内容になっています．それでいて基本的なエッセンスは書かれているので良書だと思います．

ふつうのLinuxプログラミング

ふつうのLinuxプログラミング Linuxの仕組みから学べるgccプログラミングの王道

上記と同じで，Linuxの基本的な概念をちゃんと理解しなおそうと思って読みました．丁寧にわかりやすく書かれていて，ほどほどのプログラミング経験がある人なら難なく読みこなせ，かつLinuxの思想や仕組みを把握しやすくなると思います．

アンダースタンディング・コンピュテーション

アンダースタンディング コンピュテーション―単純な機械から不可能なプログラムまで

コンパイラや構文解析，有限オトーマトン，チューリングマシンなどコンピュータサイエンス的な概念を，一からRubyベースで作って理解するという，基礎から勉強し直すときに役立つ本です．ただ結構筆致が回りくどく，分かりにくい部分も結構あるなぁという印象です．

Web+DB PRESS Vol.82

WEB+DB PRESS Vol.82

APIデザインの特集を読みたかったので．コンパクトによくまとまっていて良いです．

Web+DB PRESS Vol.77

WEB+DB PRESS Vol.77

AWSの特集が読みたかったので．一から丁寧に解説されていて良いです．

JavaScriptエンジニア養成読本

JavaScriptエンジニア養成読本［Webアプリ開発の定番構成Backbone.js＋CoffeeScript＋Gruntを1冊で習得！］

下の本を読んでから，こっちを読んで楽で綺麗でよいなぁと思いました．素のJSとか見たくもないです．

パーフェクトJavaScript

パーフェクトJavaScript (PERFECT SERIES 4)

一時期気が触れてJSの勉強をしようと思って読んでちょっと書いて嫌になりました．

インフラエンジニアの教科書

インフラエンジニアの教科書

わかりやすいんですが，中身が薄くて，この一冊だけだとかなり足りない印象．

エキスパートPythonプログラミング

エキスパートPythonプログラミング

この一年で，だいぶPythonを使うようになったので，文法的なところとか知らないところを抑えるために読みました．

その他

駆け出しマネージャーの成長論

駆け出しマネジャーの成長論 - 7つの挑戦課題を「科学」する (中公新書ラクレ)

いちおうお仕事だと駆け出しマネージャーなので，こういう本も読みます．中原先生は，ブログを読んでいても語り口が誠実でとても好感を持っています．この本も愛に溢れていて良いです．

現役・三井不動産グループ社員が書いた! やっぱり「ダメマンション」を買ってはいけない

現役・三井不動産グループ社員が書いた! やっぱり「ダメマンション」を買ってはいけない

今年は趣味で不動産データの分析をずーっとやってましたので，こういう本も読んだりします．裏側見せます系の本大好きなんですよねー．

おわりに

今年は分析に関する本とか読んでない気がします．来年はちゃんと論文読むか，もう少し違う分野のデータ関連本*2を読むか，とかをしようかなと思っていたりします．あと，そもそも本よりもWeb上の文献を参考にしている量のほうが多いと思うし，いい加減本のレビューをまとめをするのも結構意味が薄いなぁと思い始めているたりもします．

*1:私自身はこの本は未読ですが，マイケル・ルイスだと世紀の空売りとかライアーズ・ポーカーあたりが滅茶苦茶スリリングで面白かったです．

*2:医療統計とか，経済系とか，異分野の手法を知っておいたほうがいい気がしている今日この頃です．

PythonのnimfaでNMFを試す

PythonでNMFやるには，nimfaというパッケージを使えばよいらしいです．とりあえず使うだけなら，適当なnumpy行列vecを用意して，以下のように関数に投げてあげます．

factor = nimfa.mf(vec, seed='random_vcol', method='nmf', rank='5', max_iter=10)
res = nimfa.mf_run(factor).basis()

とりあえずシードはランダムで，手法はベーシックなnmf．何次元に削減するかをrankで指定して，イテレーション回数を決めればOKです．

nmfは関連手法が山ほどあって，ざっと以下のようになります．説明文は基本的に意訳です．正直意訳があってるかも自信はないので，こちらから元論文を読みましょう*1．

手法	概要
BD	ギブスサンプラーを使ったベイジアンNMF
BM	バイナリのMF
ICM	Iterated conditional modesを用いたNMF*2
LFNMF	局所特徴量を用いたフィッシャーのNMF
LSNMF	最小二乗法を用いたNMF
NMF	通常のNMF（更新式としてユークリッド/KL情報量，損失関数としてFrobenius/divergence/connectivityが指定可能）
NSNMF	non-smoothなNMF
PMF	確率的NMF
PSMF	確率的スパースMF
SNMF	最小二乗制約に基づく非負性を用いたスパースNMF
SNMNMF	スパース正則化ネットワークNMF
PMFCC	制約クラスタリングによる罰則MF

ということで，実際に回してみました．前回のマンションポエムデータ，218サンプル*1797次元のデータを10次元に圧縮します．これを10回繰り返して得られた10の圧縮行列に対して，適当なidを選んで類似度上位10件のデータを抜き出し，その一致度を集計しました．

結果は以下の通りで，時間がかかるものほど一致率も高く安定的な結果ということのようです．とはいえ安定しているから結果が良いかというと，パッとみた感じそんなに手法ごとに精度の差が歴然としているかといわれると，若干首をひねらざるを得ない感はあります*3．にしても，全体的に安定性はイマイチです... そしてBDとかPMFとかは，とにかく重すぎてパパッと結果も帰ってこず．まぁギブスサンプラーとか使ってたら思いに決まってるわけですが...

手法	一致度	NMFの算出にかかった時間 [秒]
NMF	16.7%	0.45秒
LSNMF	30.4%	0.94秒
BMF	17.8%	0.33秒
SNMF	40.7%	47.95秒

もう少しちゃんとパラメタとかチューニングしないといけないなぁという思いしかない．

*1:自分は読んでないです．

*2:この資料をみる限りICMは画像の復元手法として用いられているものみたいですね．

*3:単に類似性を目でみて判断しただけなので，完全にただの印象ではあります．

マンションポエムで新築マンションをクラスタリング

今回は自然言語処理の話です．それも若干不自然な言語のマンションポエムが対象になります．マンションポエムというのは，工場萌え*1の著者大山さんが提唱している，マンション広告に入っている詩的なコピーのことです*2．具体的にはこんな感じのやつです*3．

PLATINUM SHIP ここは、東京の暮らしの新しき起点。

そこは、時空をかける東京。 TOKYO NON DISTANCE

データ

さて，そんなマンションポエムですが，実はデータが公開されておりまして，先述の大山さんがGoogle Mapにまとめていらっしゃいます*4．このデータ，KMLという三次元地理情報を扱うためのXMLベースのマークアップ言語で，Pythonで適当にパースしてあげれば扱いやすいデータに落とし込むことができます．

ということで早速加工してみたんですが，結構データの抜けや欠けがあって，実は割りと地道に手を動かさないといけないことが判明しました．緯度経度とマンションポエムについては，全てのマンションについてあるのですが，1/3くらいは分譲時期や物件価格が入っていません．特に物件価格は注視したいところなので，いちいちググって手で入力するという作業を繰り返しました．マンションの公式サイトは，大半のものが全物件成約した段階で消されてしまうため，過去に分譲された物件は定番のマンションコミュニティのほか，こちらやこちらといった，個人で実際の広告や価格表をまとめている方のものを参考に大まかな価格を入れていきました*5*6．

あと全然関係ないですが，マンションコミュニティといえばこの覆面座談会，めっちゃ面白いです．

そんなこんなで，おおむね2005年以降に建設されたマンション218件を対象として，マンションポエムを分析しました．

いろいろと前処理

やったこと自体は非常にベーシックで，手作業で作成した辞書を使ってポエムを形態素解析して，名詞/動詞/形容詞だけを取り出して*7，ストップワード*8を削除した上で，LDAにかけてトピックを抜き出しました．

形態素解析

定番の形態素解析エンジンMeCabを使って，ざっくりと形態素に分けます．Pythonバインディングを使って，以下のようにざーっと処理を行っていきます．

tagger = MeCab.Tagger()
node = tagger.parseToNode(sentence) # 対象の文字列を読み込んで解析する
words = []
while node:
    features = node.feature.split(",") # 解析結果をカンマで区切って取り出す
    if features[0] in ['名詞', '動詞', '形容詞']: # 名詞と動詞と形容詞だけが対象
        word = node.surface.decode('utf-8') if features[6] == '*' else features[6] # 活用語の場合は基本形を取り出す
        words.append(word.lower())
    node = node.next

辞書

ここが地味に一番だるかったところで，上記の形態素解析をうまく動かすためには，ただしく形態素に分けられないといけないわけです．しかしマンションポエムはかなり妙な日本語を使っており，一般的な日本語辞書だけではうまく分割できない例が多々あります．例えば「駅徒歩圏」ですが，普通に形態素解析してしまうと「駅」「徒歩」「圏」に分かれてしまいます．これは「駅徒歩圏」ではじめて意味をなす単語なので，分割されてしまっては困ります．他にも「庭苑」「緑景」「美邸」「逸邸」といった造語に近いような単語が頻発するので，これを逐一辞書に登録していく必要があります．今回の単語をまとめた辞書は，githubに公開しているので，もし使いたい方がいらしゃればご自由にどうぞ．

同様に「神保町」や「幕張ベイタウン」みたいな地名も，ちゃんと辞書登録してあげないと無残に分割されてしまいますので，これもセコセコと手作業で登録します．地味ですがこれが一番精度に効いてくる気がしています．今回は面倒なので，動詞に関しては辞書登録するのをサボってしまいましたが*9，例えばマンションポエム頻出語として「手に入れる」があって，これが「手」「に」「入れる」に分割されてしまうと，全然ダメだったりします...

ストップワード

上記処理を経ても，意味のない言葉や頻出語は残ってしまうので，これをストップワードとして処理します．辞書と同様にこちらにありますので，興味のある方はご覧ください．

モデル作成とクラスタリング

Latent Dirichlet Allocation (LDA) によるトピックモデルの作成

上記前処理を経て形態素に分けられたマンションポエムを，bug-of-words*10のベクトルに直して，LDAにかけてトピックを抜き出します．このあたりはgensimパッケージを使ってやってしまいました．パッケージ自体の使い方は，ゆうくさんの記事やsucroseさんの記事をご覧ください．またLDAの理論的な解説については，例によってあらびきさんの解説を読んでいただくのが良いと思います．

この処理の意味するところをざっくりとまとめると，たくさんの文書から，文書間に共通するいくつかの潜在的なトピックを抜き出す，というようなことをしています*11．このトピックを使ってクラスタリングを行います．コード的には以下のような感じになります．

# 辞書の作成
dictionary = corpora.Dictionary(words)
dictionary.filter_extremes(no_below=10, no_above=0.4)
# コーパスの作成
corpus = [self._dictionary.doc2bow(w) for w in words] # bag of words
corpus = models.TfidfModel(corpus)[corpus] # apply tf-idf
corpora.MmCorpus.serialize(corpus_path, corpus) # save serialize data
# LDAモデルの作成
model = models.LdaModel(corpus=corpus,
                        num_topics=20,
                        id2word=dictionary)

Gaussian Mixture Model (GMM) によるクラスタリング

上記のLDAによるトピックは，実質的には次元削減と似た効果を持ちます．各ポエムにおける文章の出現数をまとめたベクトル（今回の例だと約1700個の単語があったため，1700次元になります）から20個程度のトピックを抜き出すことは，1700次元を20次元に圧縮することとほぼ同義です．そこで，LDAにより得られた20次元のベクトルを用いて，マンションポエムのクラスタリングを行います．別にk-meansでもいいんですが*12，今回はGMMによるクラスタリングを行いました．GMMの理屈についてはsleipnir002さんの解説をご覧ください．いくつかクラスタ数を試した結果，4つが良さげであることがわかりました*13．

可視化

坪単価 × 都心からの距離

ということで，作成した4つのクラスタについて，あとはRにかませてざっと集計&可視化してみました．dplyr使ってパイプっぽく処理して，クラスタごとのサンプル数，平均坪単価，都心からの距離平均（km）*14をまとめたのが以下になります．クラスタ1が最大派閥で90サンプル，価格はクラスタ0が平均坪単価277万でトップ*15，距離もクラスタ0と3が10km程度でした．

> # summarize
> data %>% 
+   group_by(cluster_id) %>% 
+   summarize(n=n(),
+             price.mean=mean(price),
+             distance.mean=mean(distance))
Source: local data frame [4 x 4]

  cluster_id  n price.mean distance.mean
1          0 35   277.0857      10.62453
2          1 90   257.3667      13.93625
3          2 47   228.8511      13.91467
4          3 46   255.6522      10.02251

続いて散布図に全サンプルをプロットしてみました．比較的綺麗に，都心からの距離と坪単価が逆相関しているのがみて取れます．一部に坪単価600万とか800万とかあるのは，いわゆる億ションなので，本当は外れ値にしちゃった方が良いかもですけどね... クラスタごとにみても，そこまではっきりとではないですが，まぁそこそこは綺麗に区分けされているかなぁと思います．

クラスタごとのワードクラウド

続いて，今回の素性ベクトルを作るもととなったマンションポエム形態素を可視化するために，ワードクラウドを作成します．これはTagexedoというサービスを使いました*16．現状日本語で自由な文字列からワードクラウドを作れるのは，このサービスしかないっぽいです．

ということで，各クラスタの形態素をワードクラウド化してみた結果が以下の通りです．

クラスタ0

もっとも平均坪単価が高く，都心からの距離が近いだけあって，「都心」「高台」「空」「空間」などの，都心の閑静な高級住宅地や超高層のタワーマンションなんかを表していそうなワードが並んでいます．他には「世界」「美」「上質」といったいかにもなワードも散りばめられています．

クラスタ1

こちらはもっとも都心から遠いクラスタだけあって，「住」「街」「暮」「生活」といった，そこでの具体的な地に足のついた生活をイメージさせるような単語が並んでいます．他にも「家族」「快適」「自然」「伝統」といったワードが並びます．「伝統」「継承」といった言葉は，おおむね古いランドマークがあるような場所なので（例えば浅草や鎌倉のような場所ですね），都心から離れた高級住宅地も含まれている感じがします．

クラスタ2

最も坪単価が安いだけあって，物件ではなく「家族」が一番大きな言葉となっています．家族生活に適した郊外の住宅地，というイメージでしょうか．トップ3が「家族」「暮らし」「住まい」となっているのが，いかにもという感じです．「レジデンス」については，以下のような文言が典型的な使われ方です．

風と緑のレジデンス「×××」*17駅徒歩8分、賑わいの先に広がる住宅地。駅近ながら住宅地の空気に包まれた、家族の暮らしがここから始まる。都市の利便性を手中にしながら、開放的で落ち着きに満たされたこの地。家族を優しく包み込む、ゆとりを育む、この地ならではのレジデンスを創り上げます。住宅地の最前席に住まう駅徒歩圏では限られた、第一種住居地域。

クラスタ3

最後は「tokyo」「東京」「都心」とひたすら首都をアピールしていますね．「金町」「豊洲」といった地名が並んでいるところからも，都心からの距離が近く，それでいて極端な高値ではないゾーンの物件が中心になっていることが伺えます．臨海副都心や，都心近くで再開発した東京のイメージの薄い場所などが，敢えて東京と強調する，ということなのもしれません．

終わりに

まずは，改めて大山さんにお礼を申し上げます．元となるデータがあることの意義は本当に大きいです*18．

そして今回のネタに関連して，いくつか検討したところとか，まだできそうなこととかがあるので，そのうち続きを書こうかなと考えています．このネタを書くために19日から冬休みに入ってて，ようやく今日これが仕上がりました．とはいえ，まだ一般的な冬休みの初めにすらさしかかっていないので，冬休み使ってもう1つ2つできないかなーとか... 本当はこのデータを加工して簡単なWebサービスにしようかと思ってたので，気が向いたらそれやります．あとデベロッパーとか物件の戸数とか階数とかのデータをまとめたら，もう少し面白いことできないかなと思ったり...

*1:私自身，この本を読んで工場夜景クルーズに行ったりもしました．夜の工場いいですよね．FFでいうと7のミッドガルみたいな感じで素敵です．

*2:日経にも取り上げられたりしています．

*3:個人的にマンション広告が好きで，いつの日からかマンション建設予定地に光の柱が立つようになったのを微笑ましくみていたりします．

*4:このデータ作成にかかった多大な労力に，感謝の意を表しつつ利用させていただきます．

*5:物件の広さに影響されるとアレなので，実際にはマンションの価格ではなく，坪単価の形でまとめています．

*6:もちろん同じ物件でも，高層階の東南角地と低層階の形が悪い部屋とでは価格が異なるわけで，複数の坪単価がある場合には，概ね真ん中くらいの数字を使うようにしています．

*7:今回のマンションポエムだと，助詞や副詞がほとんど意味をなさなそうなので無視した，というだけの話で，全ての場合で削除するべきだ，ということではありません．

*8:あまりに一般的すぎるので，実際に自然言語処理を行うときの妨げとなるワードのことです．googleでも，英単語の"a"や"the", "for"なんかは検索の際にストップワードとして無視されています．

*9:動詞や形容詞だと，すべての活用形をちゃんと登録してあげないといけないのですが，これが1語につき10くらいの活用形を登録する必要があって大層面倒なのです．

*10:単語の並び方は考慮せず，単に文章に単語が含まれているかどうかだけをみて，その回数を数字とした形です．

*11:因子分析で，各パラメタの背後にある共通因子を抜き出す，というのと同じイメージかなと思います．まぁモデルの数理的な形は全く違いますが...

*12:とはいえ実際のところはk-meansは精度悪いんで，あんま使いたくはないってことなんですけどね...

*13:5個以上だと似たようなクラスタがいくつもできたり，サンプル数が極端に少ないものができたり，でイマイチだったので...

*14:ここでは都心として皇居の緯度経度（35.6825° N, 139.7521° E）を用いています．ある程度の目安にしかなりませんが，まぁそこは目を瞑ってください．

*15:坪単価が想像できない人のために補足すると，坪単価277万のマンションというのは，親子4人で住むような70平米のマンションで5876万円!! もする価格です．新築マンションってタカイデスネー

*16:利用するためにはMS Silverlightが必要で，しかもChromeだと動きません...

*17:具体的な物件名は特定できない形にしておきます

*18:とはいえ草の根でデータまとめるの本当に大変なので不動産屋さんとかデータ提供してくれたら嬉しいのになぁとか思っています．まぁ企業側に何のメリットもないから無理でしょうけど...

Yosemiteでのmecabのパス

mecabの辞書作成をしようとしたら，なぜかgooglecodeにある説明とパスが違ってて，探すのに多少苦労したのでメモ．Yosemitenのせいなのか，何か自分がインストールのやり方を間違えたのかはわからないのだけど．

# mecabのディレクトリ
/usr/local/Cellar/mecab

# ipadic辞書のディレクトリ
/usr/local/Cellar/mecab/0.996/lib/mecab/dic/ipadic

そんなわけで，辞書作成コマンドはこんな感じ．権限も管理者権限で入ってたので，sudoつけないとダメだった．

sudo /usr/local/Cellar/mecab/0.996/libexec/mecab -d /usr/local/Cellar/mecab/0.996/lib/mecab/dic/ipadic -u my_dictionary.dic -f utf-8 -t utf-8 my_dictionary.csv