横須賀の某Prisonで働く研究者?のブログ

NAISTを出て横須賀の某社の研究所で研究者をしています.本ブログの内容は業務内容や所属企業とは一切関係ありません.

【論文紹介】Black Box FDR (ICML 2018)

こんにちは,こんばんは,接点QBです. 今回はICML 2018 に投稿された Black Box FDR (ansey et al., 2018) という論文の解説をしたいと思います.

概要

データを集めるために色々な実験をしますが,その際に実験ごとにサンプルの条件が異なっていることがあります. たとえば,癌細胞に作用する薬の効果を測定するための実験を考えます(図1). 30個の実験結果があるとします. 各実験においては,同一の遺伝子を持つ癌細胞に対して処置を行う(薬を投与)という試行と処置を行わないという試行を複数回行って統計量を算出します. しかし,実験が30個の実験全てで同一の遺伝子を持つ癌細胞を用いているわけではありません. つまり,実験ごとに用いている癌細胞の遺伝子は異なっていて,かつ遺伝子の差異が実験結果(薬が効くかどうか)に影響を及ぼす可能性があります.

f:id:setten-QB:20180901123034p:plain
左側が30個の実験,左側が各実験に用いた細胞の遺伝子に関する情報

このように複数の実験を行った場合で,実際に処置に対して効果があったのかを判定し,かつその実験結果に影響を与える変数を決定するための方法が Black Box FDR (BB-FDR) です.

問題設定

$n$ 個の実験結果が与えられていて,各実験から検定統計量 $z_1, \cdots, z_n$ が計算されているとします. さらに,各実験ごとに$z_i,\, i=1,\cdots, n$ に影響を与えると考えられる変数 $X_{i1},\cdots, X_{im}$ が存在します.

BB-FDR の目的は,(1) 処置に対して統計的に有意な反応を示した実験(つまり,棄却域に落ちる$z_i$)を同定し,(2) 実験結果に影響を与えている$X_{im}$を求めることです. Tansey et al. (2018) では (1) と (2) の2段階に分けて手法の設計が行われています.

各実験に対して以下の仮説検定を考えます:

$H_0$(帰無仮説):処置は効果がない      $H_1$(対立仮説):処置は効果がある

実験 $i$ で$H_0$が棄却される場合を$h_i=1$,$H_0$が棄却されない場合を$h_i=0$で表し,帰無分布を$f_0$,対立分布を$f_1$で表します.

Stage 1: 統計的に有意に処置効果が認められる実験結果の同定

検定統計量$z_i$に対して,以下のようなモデルを入れます:

\begin{align} & z_i \sim h_i f_1(z) + (1-h_i) f_0(z_i) \\ & h_i \sim {\rm Bernoulli}(c_i) \\ & c_i \sim {\rm Beta}(a_i, b_i) \tag{1} \\ & a_i, b_i \sim G_{\theta, i}(X). \end{align}

$G_{\theta, i}$ は$\theta$をパラメータに持つブラックボックス関数で,ここではDNNを用います.

$\theta$を最適化するために以下の最適化問題SGDで解きます:

\begin{equation} \min_{\theta} -\sum_i \log p_{\theta} (z_i) + \lambda \left \| G_{\theta}(X) \right \|^2_F \end{equation}

ただし, \begin{align} p_{\theta}(z_i) =& \int_{0}^{1} \sum_{h_i \in {0, 1}} p\left( z_i, h_i, c_i\right)\, dc_i \\ = & \int_{0}^{1} \sum_{h_i \in {0, 1}} p(z_i | h_i, c_i) p(h_i) p(c_i)\, dc_i \\ =& \int_{0}^{1} \left \{ c_i f_1(z_i) + (1-c_i) f_0(z_i)\right \} {\rm Beta}(a_i, b_i) \, dc_i \\ \end{align}

であり,$\left \| \cdot \right\|$ はフロベニウスノルムを表します. この最適化問題を解くことで$\hat{\theta}$が得られると,以下のように対立仮説の下での各検定統計量の事後確率を計算することが出来ます:

\begin{align} \hat{w}_i :=& p_{\hat{\theta}}(h_i=1 | z_i \sim f_1)\\ =& \int_{0}^{1} \frac{c_i f_1(z_i) {\rm Beta}(c_i | G_{\hat{\theta}, i}(X))}{c_i f_1 (z_i) + (1-c_i) f_0(z_i)}\, dc_i . \end{align}

この$\hat{w}$を降順にソートして,最初の$q$個の仮説検定において帰無仮説を棄却します. つまり,以下の最適化問題を解きます: \begin{equation} \max_{q} q \quad \textrm{subject to} \quad \frac{\sum_{i=1}^{q} \left( 1-\hat{w}_i\right)}{q} \leq \alpha . \end{equation} 要するに,$q$ 個の実験における type I erro の確率の平均を $\alpha$ 以下に抑えて,その下で検出力を最大化するというネイマン・ピアソン流の検定の構成方法をとっています.

Stage 2: important variables の同定

モデルに $G_{\theta, i}(X)$ を用いることで,以下のトレードオフが発生します:

  • メリット:変数$X$と検定統計量の関係を表すための関数のクラスが大きくなる(仮説集合のサイズが大きくなる)
  • デメリット:stage 1 の検定の検出力が(linear modelを使った場合に比べて)小さくなる

ただ,クラスが大きくはなりますが black box 的なDNNから重要な特徴量を同定するということは簡単ではありません. 特に,今回のように false discovery rate をコントロールしながらの変数選択は著者らが探した範囲では手法が提案されていないそうです.

このような問題にチャレンジするために,BB-FDR では conditional randomization tests (CRTs) (Candes et al., 2018) を用います. CRTでは,特徴行列 $X$ の第 $j$ 列に対応する変数(特徴量) $X_{\cdot j}$ を他の変数 $X_{\cdot -j}$ のみを用いてモデリングするという方法です. CRTによって,条件付き分布 $\mathbb{P}(X_{\cdot j} | X_{\cdot -j})$ が仮説 $X_{\cdot j} \mathop{\perp\!\!\!\!\perp} Z|X_{\cdot -j}$ を検定する際の妥当な帰無分布を表します. ただし,$Z$は検定統計量です. 対応する $p$ 値は条件付き分布からのサンプリングによって計算されます: \begin{align} p_j =& \mathbb{E}_{\tilde{X}_{\cdot j} \sim \mathbb{P}(X_{\cdot j} | X_{\cdot -j}) } \left[ \mathbb{I} \left[ t(\mathbf{z}, X) \leq t (\mathbf{z}, (\tilde{X}_{\cdot j}, X_{\cdot -j})) \right] \right] \\ =& {\rm Pr} \left( \left. t(\mathbf{z}, X) \leq t (\mathbf{z}, (\tilde{X}_{\cdot j}, X_{\cdot -j})) \right| \tilde{X}_{\cdot j} \sim \mathbb{P}(X_{\cdot j} | X_{\cdot -j} ) \right) \end{align} ここで,$t$ は検定統計量です. 一旦全ての変数に対して $p$ 値を求めれば,あとは標準的な Benjamini-Hochberg 法を適用して重要な変数を同定すれば良いことになります.

BB-FDR は,どの変数が帰無分布から生成された $z_i$ の事後確率の変化に関連するのかを検定します. これは,事後確率の負のエントロピーを検定統計量として用います: \begin{align} t(\mathbf{z}, X) = \sum_{i} \hat{w}_i \log \hat{w}_i + \sum_{i} \left( 1 - \hat{w}_i \right) \log \left( 1 - \hat{w}_i \right) . \end{align} 処置効果を予測するために有用な変数は,事後確率のエントロピーを減少させます.

評価実験

$P(X), P(h=1 | X)$ の組み合わせ3パターンと,検定統計量$z$ の確率分布3パターンで計6パターンでシミュレーションを行って評価しています. $X=(X_1, \cdots, X_m),\, m=50$のうち実際に目的変数に影響しているのは25個の変数だけで,false discovery rate の閾値は$10\%$ に設定されています. また,$n$ は $[100, 5000]$ で動かします.

比較手法は

  1. Benjamini-Hochberg method (Benjamini & Hochberg, 1995)
  2. NeuralFDR (Xia et al, 2017)
  3. Eq. (1) で $c_i$ の分布に fully-Bayesian logistic regression を使用したモデル

の3つです.

Stage 1 の評価

全体的に 3. fully-Bayesian logistic regression を使った手法が高い性能を指名しています. しかし,$n$ が大きくて,$z$ の帰無分布と対立分布が近いときは BB-FDR が最も高い性能を示しています. また,手法3 は学習に数時間必要で,BB-FDR は数分で学習が完了したそうです.

Stage 2 の評価

Benjamini-Hochberg method と NeuralFDR は重要な変数の同定に関しては扱っていないので, ここでは 手法3 との比較だけになります. $z$ の帰無分布と対立分布が十分に遠い場合には両手法の間に大きな差はありません. $n$ が大きくなると,false discovery rate・検出力ともにBB-FDR が手法3 を上回っています.

まとめと所感

今回の記事では,

  • 複数の実験結果が得られているとき各実験の結果を検定し
  • 実験ごとに複数の(結果に影響しているか怪しい)変数のうち,実際に結果に影響を及ぼしている変数を同定する

手法である Black Box FDR を紹介しました. 一番最初の例で出されているように,ニーズとしては生物学や医学系が多いのかなと思います. あとは社会科学系でも使えそうかな?

モデルを記述する部分にDNNを使っていますが,話のジャンルとしてはどちらかといえば伝統的な数理統計学やメタアナリシスに近いと感じました.

参考文献

  • Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the royal statistical society. Series B (Methodological), 289-300.
  • Candes, E., Fan, Y., Janson, L., & Lv, J. (2018). Panning for gold:‘model‐X’knockoffs for high dimensional controlled variable selection. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(3), 551-577.
  • Tansey, W., Wang, Y., Blei, D. M., & Rabadan, R. (2018). Black Box FDR. International Conference of Machine Learning 2018 (ICML 2018).
  • Xia, F., Zhang, M. J., Zou, J. Y., & Tse, D. (2017). Neuralfdr: Learning discovery thresholds from hypothesis features. In Advances in Neural Information Processing Systems (pp. 1541-1550).

高次元小標本データ (HDLSS data) に関するあれこれ

いつもお世話になっております.接点QBです.←社会人になった雰囲気出てます?www

ここ数年続いているDNNブームですが,その傍らでこういう思いをした人もたくさんいらっしゃるのではないでしょうか?

「そんなにデータねぇよ!!!」

そうです.基本的にDNNを使うには大量の学習データが必要になりますが,十分な量のデータがあることはそれほど多くないのではないでしょうか? DNNが覇権を取っている分野と言えば,画像・音声・自然言語あたりがメジャーかと思いますが,この辺の分野は先人たちが大きなデータセットを公開してくれていたり, 割と簡単に(専門の人,気を悪くさせてしまったらゴメンナサイ)収集できるのではないかと思います.

しかし,例えば医療分野などはそうもいきませんし,他にも大量のデータを集めることが困難な分野は多々あります. そして更に困るのが,データのサンプル数 $n$ よりもデータの次元 $d$ が圧倒的に大きい — $d>>n$ — というシチュエーションにも遭遇します. 今回は,このように $d>>n$ というシチュエーション (High Dimension, Low Sample Size; HDLSS) でデータを扱う手法をいくつかザックリと紹介したいと思います.

学部生の頃の私の専門がこの分野だったので,少し理論屋っぽい話が多くなるかもしれませんが, あまり数理的な詳細には立ち入らずに,ザックリと紹介していけたらと思います. また,DNN界隈の人の感覚で言うと「古い」話が多いので,「最新手法こそ至高」というタイプの方には微妙な記事かもしれません.

そもそも何でHDLSS だとマズイの?

そもそも,HDLSS だと何がマズイのかという話ですが,以下のような問題がよく言われます:

  1. 次元の呪い
  2. オーバーフィッティング
  3. 共分散行列周りの問題

次元の呪いについてはカステラ本を始めとする機械学習に関する本なら必ず載っていると思いますので, ここでは詳細は割愛します. まあ,「次元が大きくなるほど空間的な問題で分類や回帰が難しくなる問題」程度に考えておいてもらえればよいかと思います.

オーバーフィッティングについてもググればすぐに出てくるので割愛します.

共分散行列周りの問題についてだけ少し解説しておきます. 古典的な機械学習手法では,サンプルから標本共分散行列 — 共分散行列の推定値 — やその固有値固有ベクトル,さらにその逆行列を計算することがあります. このようなとき,$d>>n$ という状況下では

  • 標本共分散行列の逆行列 $S^{-1}$ が存在しない
  • 共分散行列の固有値推定が不安定

といった問題がよく取り上げられます(印象ですが…). とりあえず,上記のような色々な問題がHDLSS data での機械学習には存在するということを認識してもらえればOKです.

ちなみに,統計的推測のお話だと,古典的な統計的推測理論は $n\to \infty$ かつ $d$: fix という枠組みで作られていたため, HDLSS でよく設定される $d\to \infty$, $n$: fix という漸近理論の設定では原則として使えないという問題もあります. もう少し具体的に言うと,正規性を仮定できないことが普通といった状況になるということです.

本論

では,HDLSS データに対してどのような手法を使えばよいのかと言うと,基本的なアプローチは以下のようなものがあります:

  • 特徴選択による特徴ベクトルの次元削減
  • PCA等で新しい特徴量を作成
  • 特徴ベクトルはそのまま使って,標本共分散行列の計算や固有値問題の解き方を工夫する

2つ目と3つ目は背反ではないことも多いですが,大雑把には上記の3つがメジャーです.

特徴選択に関する手法あれこれ

これに関しては2つの手法を紹介します.

Gradient Boosted Feature Selection

Xu et al. (2014) によって提案された手法で,中身はタイトルの通り. これについてはTJO氏@TJO_datasci詳細な解説 があるので,そちらを参照するのが良いかと思います.

HSIC Lasso

HSIC Lasso (Yamada et al., 2012) は,大雑把に言うと Hilbert-Schmidt information criterion を基準として用いて特徴選択を行う手法です. 再生核ヒルベルト空間上でLassoに対応する特徴選択をすることで,特徴量と目的変数の間や特徴量同士の間に非線形な関係がある場合でも,良い特徴選択が可能になるようです. 現在はこれの改良版 も出ているようです (Yamada et al., 2018). また,実装は著者の山田先生 (@myamada0)のページ に公開されています.

DNP

DNP (Liu et al., 2017) では DNN の学習中に,勾配情報を用いて特徴選択とDNNの学習を同時に行う方法が提案されています. 「gradient の大きい特徴量ほど loss を急激に減少させているので,そのような特徴量こそが重要なんだ」という発想ですね.

PCAによる特徴量作成

ここで紹介したい内容はいくつかあるんですが,かなり数理的な内容になってしまうので,とりあえず大雑把にどういう枠組みで手法が作られているのかの説明だけしておきます. まず,共分散行列の固有値を $\lambda_i, \, i=1,\cdots, d$ とします. このとき,固有値に対して以下を仮定するモデル (spiked model) が Johnstone (2001) によって提案されました.

$\lambda_i > 1,\, i=1,\dots, m $ は $d$ に依存しない定数で, $\lambda_i = 1,\, i=m+1, \dots, d$.

私が知っている範囲だと,spiked model を仮定してHDLSS data の幾何学的特徴に着目して固有値を推定し,PCAを行うという方法が提案されています (Yata and Aoshima, 2010; Aoshima and Yata 2011). 現在のところ,このモデルや,これをベースとした power spiked model (Yata and Aoshima, 2013) が仮定されて理論構築がされている手法が多いのではないかと思います.

固有値問題の解き方で工夫する方法

これは結構有名で,固有値問題や一般化固有値問題を解く時に対角成分だけを使ったりします. Bickel and Levina (2008) あたりが有名でしょうか. 他にも,いくつか仮定を入れて固有値固有ベクトルの漸近表現を陽に書けたりします. が,この内容は某ジャーナルで現在査読中のため詳細は伏せます.

まとめ

今回の記事では,HDLSS データに関することを適当にピックアップして紹介しました. 世の中DNNがもてはやされていますが(実際に凄いですけどねw),実務ではいつもデータがあるとは限らない,むしろデータが潤沢にある方が稀かと思います. そのようなとき,本記事で紹介した内容をふと思い出して論文に当たっていただけると幸いです.

なお,本記事の内容は殆どが学部生の頃の知識かつ,特にきちんとしたサーベイもせずに書き上げてしまったので,間違えているところがあるかもしれません. その時は御指摘願います.

  • 頂いたコメントを下に一部の記述を修正しました.

本文中で引用した論文リスト

  • Aoshima, M., & Yata, K. (2011). Effective methodologies for statistical inference on microarray studies. In Prostate Cancer-From Bench to Bedside. InTech.
  • Bickel, P. J., & Levina, E. (2008). Regularized estimation of large covariance matrices. The Annals of Statistics, 199-227.
  • Liu, B., Wei, Y., Zhang, Y., & Yang, Q. (2017, August). Deep neural networks for high dimension, low sample size data. In Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17 (pp. 2287-2293).
  • Yamada, M., Jitkrittum, W., Sigal, L., Xing, E. P., & Sugiyama, M. (2012). High-dimensional feature selection by feature-wise non-linear lasso. arXiv preprint. arXiv preprint arXiv:1202.0515.
  • Yamada, M., Tang, J., Lugo-Martinez, J., Hodzic, E., Shrestha, R., Ouyang, H., ... & Saha, A. (2018). Ultra high-dimensional nonlinear feature selection for big biological data. IEEE Transactions on Knowledge and Data Engineering.
  • Yata, K., & Aoshima, M. (2010). Effective PCA for high-dimension, low-sample-size data with singular value decomposition of cross data matrix. Journal of multivariate analysis, 101(9), 2060-2077.

【NAISTでの生活】就活編―IT/ICTベンチャーから日系大手まで―

どうも,接点QBです. 今回はNAISTでの生活シリーズ第2回目(3回目は無さそうな気がする)として,就活に関する記事を書きます. 就活を控えたNAIST生・NAISTの受験を考えている学部生・IT/ICT企業を目指している人達の参考になればと思います.

(注)企業様の敬称は省略します.また,多少の表記ゆれがあるかもしれません.

(注) 本記事は2018年1月に執筆したものです.

お前誰よ?

就活ということで,履歴書やエントリーシートに書きそうな内容を列挙しておきます. 多分この辺の情報を書いておかないと参考にならないので.

どこ受けたの?

実際に選考を受けたのは,以下の企業様です:

  • Albert:データソリューションカンパニー.データサイエンティストを多く抱えている会社.
  • Brain Pad:こちらもデータ分析コンサルみたいな会社.Albertとの違いは後述.
  • サイバーエージェントアメブロ・AbemaTV,FRESH!を始め,色々なコンテンツを展開している有名ITメガベンチャー
  • DeNA:キュレーションサイトの事が新聞に載る前に受けました.
  • NTTドコモ:携帯電話会社という印象が強いと思いますが,意外と色々やってます.
  • NTTコミュニケーションズ:通信回線の保守運用しかやってないという印象かもしれませんが,こちらもデータ分析をやっている部署があります.
  • ワークスアプリケーションズ:BtoBの会社なので知ってる人は少ないかもしれません.Hueという人工知能ERPが最近のウリのようです.
  • Yahoo! JAPAN:こちらも有名メガベンチャー.インフラエンジニアからデータサイエンティストまで,色々な人がいらっしゃいます.

大学院に入った当初は金融系も受けようと思っていましたが,Fintech系の話をウォッチしているうちに「もう日本の銀行は駄目だわ」と思って受けるのをやめました. が,ここではその話には深く立ち入らないことにします(長くなっちゃうので).

それで,実際にどこから(内)内定を貰ったかというと,DeNA以外の全社です. ということで,以下では(内)内定をゲットした会社に関する話を(公開できる範囲で)書いていきます.

M1の12月まで

学内での説明会

NAISTでは色々な企業の方に来ていただいて,学内で企業説明会みたいなものが行われます. もちろん,経団連に所属している企業様は3月までは説明会も出来ませんので,この時期は経団連に属していない企業様のみの説明会になります. 上に上げた企業様だと,Albert,NTTドコモNTTコミュニケーションズワークスアプリケーションズ以外はNAISTで説明会があったと思います. 大体は サポーターズキャリアセレクト といった「就活のサポートしますよ」企業様がイベントを設定してくださいます.

ブレインパッドはサポーターズのイベント経由で選考フローに乗りました.

逆求人イベント

さて,僕が一押しするのが「逆求人イベント」です. 逆求人を知らない方もいらっしゃるかもしれないので,一応説明しておきます. 逆求人イベントでは,普通の説明会と違って学生がブースを持ちます. 学生が個人ごとに設定したブースに座っていて,企業の方が学生のブースにいらっしゃって,お話をします.

僕は12月の中頃に 逆求人ナビ が主催する「データ逆求人」というイベントに参加しました. このイベントなくして僕の就活に成功は無かったと言っても良いぐらい素晴らしいイベントでした. 「データ逆求人」はデータ分析業務に就きたい人(俗に言うデータサイエンティスト)を対象とした逆求人イベントだったと記憶しています. 当日は時間が6つコマに区切られていて,各コマに1つの企業様と話せるというものでした. 幸運にも6コマ全て埋まったので朝から夕方まで喋りっぱなしで喉が結構きつかったです.

M1の1月以降

選考が本格化します. ここからは企業ごとに書いていこうと思います.

Albert

1月にアカリクの逆求人イベントに参加してお話しました. アカリク逆求人の前までは全く知らない企業だったんですが(最初は自転車の会社かと思いましたw),イベントでの説明を聞いていると「割りと興味あるな」となって,幸運にもお話できました. イベント中に「もう殆ど内定出たと思っといて」と言われて「マジかよwww」と思いましたが,結構マジでした. 選考フローは面接2回だけ(一回目は現場のデータサイエンティストの方と「何やりたい?」という感じの話で,二回目は雇用条件に関するお話)でした. つまり,殆ど選考らしい選考はありませんでした(逆求人で色々話したおかげ?).

Albert はかなりブレインパッドと事業内容が被っていますが,あえて違いを上げるとしたら

  • Albertの方がブレインパッドよりもアカデミックな印象を受けた
  • (あくまで個人的な印象ですが)Albertの方が先端技術をウォッチしている

という点でしょうか.

給与面に関しては月額346,761円 (30時間分のみなし残業代込み),近くに住むと家賃補助3万円/月,通勤手当は月3万円を上限に定期代で支給,その他も諸々の手当が付く感じです. ボーナスは年2回あるみたいです(何ヶ月分出るかは知りません). 基本的に年俸が500万をきることは無さそうな感じでした.

ブレインパッド

Albertよりもビジネス寄りな印象を受けました. 人事の方は良い人でしたし業務内容も良かったんですが,少し待遇が悪すぎましたね. 給与が月額28万円(30時間分のみなし残業代込み),家賃補助なし,ボーナスも決算賞与(つまり,出るかわからない)という感じでした.

選考は,Skypeでの面接が数回と,白金台にある会社での最終面接が1回でした.

サイバーエージェント

こちらも1月のアカリク逆求人経由で選考フローに乗りました. Skypeでの選考では,アドテクスタジオのAILab所属の方や,秋葉原ラボの方とお話させて頂きました. サイバーエージェントというとキラキラ系女子のイメージがある方が多いと思いますが,実際に中を見てみると「ベンチャーなのに意外ときっちりしている」と感じました. アドテクというかなり強固な収益基盤を持っていて,新規事業も展開し,かつ数年に一度きちんと事業を整理しているという,イケイケベンチャー感からは想像出来ないぐらいキチンとした経営戦略をとっている企業です.

選考フローはwebテストSkypeでの面接数回→本社(渋谷)での役員面接一回でした. webテストの詳細はさすがに書けませんが,僕は全然出来ませんでしたw(分野が違いすぎた) でも,webテストを受ける時にプロフィールシートも一緒に出すので,こちらをキチンと書いておけば何とかなる(会社側のニーズにマッチする内容なら)と思います. また,Skypeでの面接は「誰と話したい?」とこちらの要望をかなり聞いていただけました.

結局内定は辞退したんですが,正直な所かなり迷いました.

NTTドコモNTTコミュニケーションズ

秘密ですwww

他のITベンチャーに比べて古くからある「The 日系大手」な感じの企業なので,あまり書くと怒られるかと思いましてw

一つだけアドバイスするなら,最終の面接はスーツで行きましょうw

ワークスアプリケーションズ

逆求人ナビのイベント経由. ここが一番年俸が高かったです.最低でも600万円は保証すると言われました. 職種としては「AIスペシャリスト」という職で内定をいただきました. 人事の方がとても親身に相談に乗ってくださり(他の会社に関しての話など),ネットで言われている「ブラック」という雰囲気はあまり感じませんでした. また,現場のデータサイエンティストの方ともお話しましたが,技術レベルは高い(高い人としか話していない可能性あり)と感じました.

選考フローは,(ここまで免除)→プログラミング課題→技術面接→最終面接 でした. プログラミング課題は,「データサイエンス系の研究でプログラム書いてるならコレぐらいは出来るだろうな」という感じのレベルでした. 逆にあれぐらい書けないと入ってからキツイのでは?という印象でしたね. 技術面接は結構ゴリゴリ聞かれましたが,圧迫的な感じは一切なく,純粋に技術レベルを推し量ろうとしていると感じました.

Yahoo! JAPAN

逆求人ナビのイベント経由. イベントの時に「一人だけめっちゃきれいな人いる!やばい!!」と思っていた人がヤフーの人事の人でした(Twitterで有名な例のあの人).

メガベンチャーですが先進的な取り組み(社内システム含む)に積極的で,働きやすさ・面白さという2点では突出しているかなと思いました(面白さはサイバーもかなり良いと思いますが). ただ,新卒で入ると如何せん給料が(上記の会社と比べて)低いんですよね. 住宅手当も付きませんし. 中途採用で入るならかなり有力な選択肢になるかなと感じました.

選考フローは,Skype面談やらwebテスト→東京本社での面接→大阪での最終面接 でした. webテストはかなり基礎的な内容で,(コンピュータサイエンスの分野で普段から研究してるなら)特に対策とかしなくても大丈夫かと思います. 東京での面接は技術系の方との面接で,ここが一番重要だと思います. 僕の場合だと,スライドを使って研究紹介プレゼンをしました. 専門の方がいらっしゃるので,結構しっかり内容について質問が飛んできます. がしかし,こちらも普段からきちんと研究発表をしていれば何とかなるかと思います.

まとめ

結構色々な企業を受けましたが,特に何か就活用の対策をしたということはありませんでした―SPIの対策とかしたこと無い―. 個人的な意見としては,理系なら「就活の対策に時間を割くよりも,研究に時間を割いて自分の研究についてキチンと話せるようにしておく」ことの方が大事ではないかと思います. もちろん,これは一概には言えないのかもしれませんが,少なくとも僕はそう感じました. せっかく修士まで来ている(しかも,我々の研究に掛かる費用の多くは国民の血税から賄われている)のですから,その専門性を活かした職に就くべきだと思います.

また,本記事では具体的な給与額を含めて,できるだけ詳しく就活内容に関して書きました. こういう記事だとあまり給与を明示しているものが多くないのですが,

雇用される側からすると給与は切実かつ重要な要素なので,学生側がアクセスしやすいようにきちんと明示すべき

という考えのもとで出来るだけ明記しました. これから就活をする人の参考になれば幸いです.

最後になりましたが,就活を通して多くの企業の方(特に人事の方)にお世話になりました. この場を借りて御礼申し上げます. もし本記事の内容に問題がありましたら,ご連絡下さい.

【NAISTでの生活】ビジネスコンテスト・シリコンバレー編

どうも,接点QBです.最近修士論文の執筆に追われています. そんな中ですが,ふと思い立って大学院での生活をトピックごとに振り返る記事を書いてみようと思いました. なお,本当は紹介しようと思っていたNIPSの論文がまだ読み込めてないから,暇つぶしに記事を書こうとか思ったわけでは無いですよ.

今回は,タイトルに有る通りビジネスコンテストとそれに関連してシリコンバレーに行った時の話を書こうと思います. これからNAISTへの進学を考えている方や,「シリコンバレーって今どんな感じ?」といった方々,「大学院生がシリコンバレーで何してきたの?」という方々に,僕のドヤ顔体験記に興味ある方(←日本に一人ぐらいいてもいい)がメイン対象です. 上記以外の方も興味があればぜひ読んでみて下さい. 写真が多いのですぐ読み終わると思います.

何でビジコン?

僕はNAIST情報科学研究科に所属している大学院生なわけですが,そんな自分が何故突然ビジコン(ビジネスコンテスト)の話を書き始めたかということをまず説明します. NAISTではGEIOT (Global Entrepreneurs in Internet of Things)というプログラムが設置されています. これは,IoT分野での起業家を育成しようというプログラムで,文科省のEDGEプログラムに採択されていました. 現在はEdgeプログラムが終了して,EDGE-nextというプログラムになっているそうです. それに伴い,crossXcrossというプログラムがNAISTで設置されたようです(こっちはあまり知らない). つまり,僕はこのGEIOTに参加していたわけですね. プログラムの殆どが大阪イノベーションハブで開催されていたので,毎週土曜日にグランフロント大阪に行ってドヤ顔する生活を送っていました.

立命館大学学生ベンチャーコンテストで優秀賞

脱落していくチームも多い中,僕はチームメンバーに恵まれてなんとか製品のプロトタイプを作ることが出来ました. 実際に某お店で製品の実証実験もさせてもらい,かなり上手く行ったチームだと自負していますw

チームで考えたビジネスプランで国内のビジコンに応募して,2016年の立命館大学学生ベンチャーコンテストで優秀賞をもらいました(これ,検索すると一発で本名バレしてしまうのでは…).

f:id:setten-QB:20180107174556j:plain

それで,国内のビジコンで優秀な成績を収めたということで,大学院からお金を出してもらったお金とビジコンの賞金でシリコンバレーでの研修に連れて行ってもらいました!

シリコンバレー研修

さて,本題のシリコンバレーです. 研修は大阪イノベーションハブが主催するシリコンバレー・アントレプレナー・ネットワーキングプログラム2017に参加するという形でした.

研修の中身

  • 現地の日本人起業家や現地で働いている人達の前で自分のビジネスプランのプレゼンをする
  • 起業家の人たちからレビュー
  • 起業家の人たちの講演

というのが基本的な研修形式でした.研修場所は色々で,起業家の方たちにホテルまで来ていただく時もあれば,研修参加者が講師の方が働いている場所まで行って見学+研修という形式になることも多かったです. 以下,有名所の写真と共に現地でのことを振り返ってみます.

f:id:setten-QB:20180107180439j:plain

はい,まずは皆さん大好き「Apple」ですね. 実は宿泊していたホテルのすぐ近くにAppleの本社がありました. 道を歩いているとどこもかしこもAppleの看板だらけでした.

f:id:setten-QB:20180107180949j:plain

次はFacebookです. この看板を見るとこの構図で記念写真を撮るのは万国共通の模様. ちなみに,この裏面はサン・マイクロシステムズのロゴが入っています. 時代を感じますねぇ…

f:id:setten-QB:20180107181330j:plain

次はGoogleです. Googleでは日本人エンジニアの方のお話を聞かせてもらいました. ちなみに,Googleはマウンテンビューというエリアに本社があるんですが,そこら辺一体の土地を殆ど買い占めてるそうです. 社内が広すぎるので,移動するようの自転車(Googleカラー)がそこら辺に置いてあって,みんなそれに乗って移動するそうです(写真見たい人はTwitterFacebookで聞いてね).

f:id:setten-QB:20180107181903j:plain 次はスタンフォード大学です. 全然大学っぽくないですねw 個人的には,ここの医学部の先生から「リーダー像とは」みたいな講演をしていただいたのが一番面白かったです. ちなみに,医学部キャンパスの食堂でお昼を食べたんですが,noodleが不味すぎて衝撃でした. 皆さん,スタンフォードの学食ではnoodle以外を食べましょう.

f:id:setten-QB:20180107182304j:plain

最後はこちら,サンフランシスコの日本領事館. 当時,大阪市とサンフランシスコ市は姉妹都市だった(はた迷惑なお隣さんから来た人たちが頑張って某像を設置したおかげで,提携を解消すると大阪市長は言っていますが,どうなるのでしょうか?)ので,その繋がりもあり最終日は日本領事館でプレゼン+懇親会でした. アメリカでは色々食べましたが,ここで出てきた料理が一番美味しかったです. というか,日本と比べるとどこの国も食べ物が美味しくない. 別件で中国に行った時に比べるとアメリカはマシでしたが…

GEIOT勢は優秀

研修中のプレゼンは全て採点されて,最終日に順位発表があったのですが,上位3人は全てNAISTから参加した人でした. ちなみに,僕は2位でした( ・´ー・`)どや 自分で言うのも何ですが,NAISTのGEIOTから研修に参加した人たちは優秀な人ばかりで(そりゃ,みんな国内のビジコンで優秀な成績を収めた人だからね),心なしか引率の先生もドヤ顔でしたw 京大をはじめととする他の大学からの参加者や,社会人の人など色々な人が参加する中でも好成績を収められたことは,GEIOTというプログラムの素晴らしさを物語ってる気がします.

まとめ

NAISTでの生活を振り返るシリーズの第1段として,ビジコン・シリコンバレー編でした. 書ききれてないことや書けないこと(具体的なビジネスプランとか)もたくさんありますので,興味がある方は個人的に聞いてください. 出来る限りお答えします.

この活動をしている中で色々な起業の方・起業家の方・現場の方・省庁の方から意見を伺うことが出来ましたが,結構好評でした. なので,もしかすると将来的には僕達のプロダクトが皆さんに使われる日が来るかもしれません(←なお,現状は微妙な模様w)

【論文紹介】Interpretable Predictions of Tree-based Ensembles via Actionable Feature Tweaking【KDD 2017】

どうも,接点QB です.

KDD 2017の論文を漁っていたら,「Interpretable Predictions」というタイトルに目を惹かれて表題の論文を読んだので紹介したいと思います. 機械学習(教師あり)の手法は「予測・分類」を目的としていて,「説明性」に関してはあまり重視されない印象があります. しかし,実際の問題を解決しようとする際には説明性を求められるケースが多いと思います. たとえば,「この人は商品Aを購入しないグループに属する事はわかりました.じゃあ,どうすればこの人に商品Aを買わせることができますか?」といったケースがあります.

今回紹介する論文で提案されている手法は,上記のような例でも対応することが出来ます.

論文の概要

問題設定と前提条件は,

  • 2値分類
  • 使用するのは決定木ベースのアンサンブル分類器

です.そして,論文の目的は

  • negativeに分類された$\renewcommand{\vec}[1]{\boldsymbol{#1}}\vec{x}$を変換して,positiveに分類される$\boldsymbol{x}'$を作成する

ことです.

Notation

特徴量ベクトル
$\boldsymbol{x}\in \mathcal{X}\subset \mathbb{R}^n$.$\mathbb{E}[\vec{x}]=\vec{0},\, V[\vec{x}]=\vec{1}$とする.
ラベル
$\boldsymbol{x}$に対して$y\in \left\{-1, +1 \right\}=\mathcal{Y}$が対応する.
特徴量ベクトルとラベルの対応を表す写像
$f:\mathcal{X}\to \mathcal{Y}$.分類器はこの写像を推定したもの$\widehat{f}$.
アンサンブル学習器
弱学習器(決定木)$T_1,\cdots, T_K$のアンサンブル学習器を$\mathcal{T}$とする.

手法

まず注意してほしいのは,本稿で紹介する手法は「分類器を作成すること」ではなく「作成した分類器を使って特徴量ベクトルをより良いものに変換すること」です. なので,分類器$\mathcal{T}$は既に学習されているとします.

$\boldsymbol{x}$がtrue negativeすなわち$f(\boldsymbol{x})=\widehat{f}(\boldsymbol{x})=-1$を満たすとします. このとき,取り組むタスク \begin{equation} \boldsymbol{x}'\in \mathcal{X} \quad \text{such that} \quad \widehat{f}(\boldsymbol{x})=+1 \end{equation} を作成することです.さらに言えば,変換されたベクトルの中でも何かしらの意味で最良のベクトルを選択するという問題に落とし込むと,求めたいベクトルは一意に定まります:

\begin{align} \newcommand{\argmin}{\mathop{\rm arg\, min}\limits} \boldsymbol{x}' = \argmin_{x^{*}} \left\{ \delta \left( \boldsymbol{x}, \boldsymbol{x}'\right) \left| \widehat{f}(\boldsymbol{x})=-1 \text{ and } \widehat{f}\left( \boldsymbol{x}^{*}\right)=+1 \right. \right\}. \tag{1} \end{align}

(1)での$\delta$は,「何かしらの基準」を測るためのコスト関数です.ユークリッド距離でもコサイン類似度でも,目的によって使い分けるのが良いと思います.

Positive and Negative Paths

決定木の根から葉ノードへの経路は,一つの特徴量$x_i$に関する不等号条件の積み重ねによって形成されます. そして,葉ノードでnegativeかpositiveを分類します.

決定木$T_k$において$\vec{x}$をpositiveと分類する$j$番目の経路を$p_{k,j}^+$で表します. 簡単のため,$p_{k,j}$は次のように表せると仮定します: \begin{equation} p_{k,j} = \left\{ \left( x_1 \substack{\geq \\ \leq}\theta_1\right), \cdots ,\left( x_n \substack{\geq \\ \leq}\theta_n\right) \right\}. \end{equation} さらに,$P_{k}^{+}=\bigcup_{j\in T_k} p_{k,j}$とします.negativeの場合の$P_{k}^{-}$も同様に定義します.

さて,このパスを指定することで入力空間から特定の領域($p^+_{k,j}$を指定すればpositiveに分類される領域)を切り出す事ができます. たとえば,$\vec{x}=[x_1, x_2]$という特徴量ベクトルを分類する事を考えます. このとき,1つ目のノードの条件が$x_1\leq \theta_1$,2つ目のノードでの条件が$x_2\leq \theta_2$だとすると, これらの条件から$\mathcal{X}_{\rm posi}=\left\{ (x_1, x_2) \left| x_1\leq \theta_1,\, x_2\leq \theta_2\right. \right\}$という領域がpositiveと分類される$\vec{x}$を含む領域である事がわかります. 図1で緑になっている領域が$\mathcal{X}_{\rm posi}$に対応します.

f:id:setten-QB:20171021204449p:plain

図1

特徴量ベクトルの構成方法

[Tolomei et al., 2016] で提案されている手法は,$\vec{x}$から近くて,$\widehat{f}\left( \vec{x}'\right)=1$となるような$\vec{x}'$を発見する方法です.つまり(1)の$\vec{x}'$を見つけるという問題になるわけですね. その具体的な方法を説明していきます. まずは$p_{k,j}^{+}$の$\varepsilon$-satisfactory instance \begin{equation} \vec{x}^{+}_{j(\varepsilon)}[i] = \begin{cases} \theta_i - \varepsilon & x_i \leq \theta_{i} \\ \theta_i + \varepsilon & x_i > \theta_{i} \end{cases} \tag{2} \end{equation} を定義します. (2) において,$\vec{x}^{+}_{j(\varepsilon)}[i]$は$\vec{x}^{+}_{j(\varepsilon)}$の第$i$要素です. また,決定木$T_k$における$\vec{x}_{j(\varepsilon)}^{+}$の集合を$\Gamma_{k}=\bigcup_{j\in P_{k}^{+}}\vec{x}_{j(\varepsilon)}^{+}$で定義します. 注意してほしいのは,$\vec{x}_{j(\varepsilon)}^{+}\in \Gamma_{k}$は必ずしも$\widehat{f}\left( \vec{x}_{j(\varepsilon)}^{+}\right)=+1$とはならないことです. あくまで$\widehat{h}_k\left(\vec{x}_{j(\varepsilon)}^{+}\right)=+1$ということしか言えません($\widehat{h}$は決定木$k$での予測ラベル).

そして,$\Gamma=\bigcup_{k=1}^{K}\Gamma_k$(i.e. すべての決定木(弱学習器)での$\varepsilon$-satisfactory instance の集合)に属し $\widehat{f}\left(\vec{x}_{j(\varepsilon)}^{+}\right)=+1$を満たす$\vec{x}_{j(\varepsilon)}^{+}$の中から, $\delta\left(\vec{x}, \vec{x}_{j(\varepsilon)}^{+}\right)$を最小にするものを$\vec{x}'$とします. つまり, \begin{equation} \vec{x}'=\argmin_{\substack{\vec{x}_{j(\varepsilon)}^{+}\in \Gamma \\ \widehat{f}\left(\vec{x}_{j(\varepsilon)}^{+}\right)=+1}} \delta\left(\vec{x}, \vec{x}_{j(\varepsilon)}^{+}\right) \end{equation} によって,$\vec{x}'$を決定します.

$\vec{x}^{+}_{j(\varepsilon)}$を決定するために探索する範囲は,決定木の数$K$とpositiveパスの本数の$|P_{k}^+|$の積で書けるので $O(K|P_{k}^+|)$になります.

実装 ―Pythonでやってみる―

(注)2017/10/25に下記内容を変更しました.デバッグが出来ていないので,もしかしたらバグが含まれているかもしれません

擬似コードは論文に記載されているので,ここではPythonで本手法を使ってみたいとおもいます. とりあえずはirisデータでやってみます. これまでの説明は2クラスの分類問題で説明してきましたが,簡単に3クラス問題に拡張できるので,3クラス分類でやりたいと思います.

まずはデータを読み込んで,Random Forestで分類器を作ります.

import numpy as np
import pandas as pd
import scipy.stats
import copy
from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
import copy
import scipy.stats


iris = datasets.load_iris()

x_arr = iris['data']
mean_x = x_arr.mean(axis=0)
std_x = x_arr.std(axis=0)
x_arr = scipy.stats.zscore(x_arr)
y_arr = iris['target']

rfc = RandomForestClassifier()
rfc.fit(x_arr, y_arr)

rfcに各弱学習器(分類器)の情報が含まれていて取り出すことが出来ます. たとえば,rfc[0]には

DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=None,
            max_features='auto', max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, presort=False,
            random_state=2072047594, splitter='best')

が入っています.そして,この決定木に関する情報を取り出して目的のラベルに至るパスや,各ノードの閾値を得なければなりません. 決定木に関する情報は,オブジェクトの.tree_メンバが持っています.

ここで,まずはwc = rfr[0]の葉ノードから根ノードに至るパスを所得する方法を考えましょう. ここで使えそうなtree_メンバでchildren_leftchildren_rightがあります. children_leftの第$i$要素には,ノードIDが$i$の左側子ノードのIDが入っています. children_rightも同様で,こちらには右側子ノードのIDが入っています. もちろん,葉ノードには子ノードが存在しないので,ノード$i$が葉ノードになっているときは$-1$が入っています. 試しにwc.tree_.children_leftを実行してみると

array([ 1, -1,  3, -1,  5, -1,  7,  8, -1, 10, -1, -1, -1])

が帰ってきます.しがたって,子ノードが$-1$になっている要素$i$を取ってくれば,子ノードのIDが分かります.

children_left = wc.tree_.children_left  # information of left child node
# select the top-maxj prediction node IDs
leaf_nodes = np.where(children_left==-1)[0]

今回ほしいのはpositiveラベル(ここではラベル$3$にします)を出力する葉ノードへのパスだけなので, leaf_nodesでラベル$3$を出力するものだけを選択します.

class_labels = [0, 1, 2]
aim_label = 2

leaf_values = wc.tree_.value[leaf_nodes].reshape(len(leaf_nodes), len(class_labels))
leaf_nodes = np.where(leaf_values[:, aim_label] != 0)[0]

子ノードIDがわかったので,適当に子ノードのIDを一つ固定して話を進めます; leaf_node = leaf_nodes[0]. では,この葉ノードから根ノードまでのパスを求めます.

children_right = wc.tree_.children_right
feature = wc.tree_.feature
threshold = wc.tree_.threshold

leaf_node = leaf_nodes[0]  # 葉ノードを固定

child_node = leaf_node
parent_node = -100  # initialize
parents_left = []  # 左側親ノード
parents_right = []  # 右側親ノード
while (parent_node != 0):
    if (np.where(children_left == child_node)[0].shape == (0, )):
        parent_left = -1  # 左側親ノードが存在しない場合は-1
        parent_right = np.where(
            children_right == child_node)[0][0]
        parent_node = parent_right
    elif (np.where(children_right == child_node)[0].shape == (0, )):
        parent_right = -1  # 右側親ノードが存在しない場合は-1
        parent_left = np.where(children_left == child_node)[0][0]
        parent_node = parent_left
    parents_left.append(parent_left)
    parents_right.append(parent_right)
    """ 次のステップへの処理 """
    child_node = parent_node

上記プログラムですが,まずはchild_node = leaf_nodeによって子ノードIDを初期化しています. また,parent_node = -100で親ノードIDを初期化しています. parent_leftは少しややこしくて,「ノード$i$を左側子ノードにもつノードIDが,第$i$要素に格納」されています(これを左側親ノードと呼ぶことにします); parent_right = np.where(children_right == child_node)[0][0]parent_rightも同様で,「ノード$i$を右側子ノードにもつノードIDが,第$i$要素に格納」されています. parent_left, parent_right共に,条件を満たすノードが存在しなければ$-1$が格納されます;parent_left=-1. 決定木ではparent_leftの第$i$要素とparent_rightの第$i$要素に非$-1$である値が同時に入ることはありません (i.e. 特定のノードへ2つ以上のノードから矢印が伸びることはない). なので,親ノードはwhile文の中のif文によって一意に決まります;parent_node=parent_right or parent_node=parent_left. この親ノードを辿っていく処理をparent_nodeが$0$になるまで繰り返すことによって(根ノードのIDは必ず$0$), 葉ノードから根ノードまでのパスを求める事ができます.これらの処理を全ての葉ノードに対して行います.

""" search the path to the selected leaf node """
paths = {}
for leaf_node in leaf_nodes:
    """ correspond leaf node to left and right parents """
    child_node = leaf_node
    parent_node = -100  # initialize
    parents_left = []  # 左側親ノード
    parents_right = []  # 右側親ノード
    while (parent_node != 0):
        if (np.where(children_left == child_node)[0].shape == (0, )):
            parent_left = -1  # 左側親ノードが存在しない場合は-1
            parent_right = np.where(
                children_right == child_node)[0][0]
            parent_node = parent_right
        elif (np.where(children_right == child_node)[0].shape == (0, )):
            parent_right = -1  # 右側親ノードが存在しない場合は-1
            parent_left = np.where(children_left == child_node)[0][0]
            parent_node = parent_left
        parents_left.append(parent_left)
        parents_right.append(parent_right)
        """ 次のステップへの処理 """
        child_node = parent_node
    # nodes dictionary containing left parents and right parents
    paths[leaf_node] = (parents_left, parents_right)

これで,pathsはkeyが葉ノードIDで,valueが(左側親ノードID,右側親ノードID)になっているディクショナリになります:

{1: ([0], [-1]),
 3: ([2, -1], [-1, 0]),
 5: ([4, -1, -1], [-1, 2, 0]),
 8: ([7, 6, -1, -1, -1], [-1, -1, 4, 2, 0]),
 10: ([9, -1, 6, -1, -1, -1], [-1, 7, -1, 4, 2, 0]),
 11: ([-1, -1, 6, -1, -1, -1], [9, 7, -1, 4, 2, 0]),
 12: ([-1, -1, -1, -1], [6, 4, 2, 0])}

では次に,パス内のノードでの分岐条件を記述する情報を取り出します. 必要な情報は,条件分岐に使われる

  • 特徴量のインデックス;features
  • 閾値thresholds
  • 不等号の向き;inequality_symbols

です.

path_info = {}
for i in paths:
    node_ids = []  # node ids used in the current node
    # inequality symbols used in the current node
    inequality_symbols = []
    thresholds = []  # thretholds used in the current node
    features = []  # features used in the current node
    parents_left, parents_right = paths[i]
    for idx in range(len(parents_left)):
        if (parents_left[idx] != -1):
            """ the child node is the left child of the parent """
            node_id = parents_left[idx]  # node id
            node_ids.append(node_id)
            inequality_symbols.append(0)
            thresholds.append(threshold[node_id])
            features.append(feature[node_id])
        elif (parents_right[idx] != -1):
            """ the child node is the right child of the parent """
            node_id = parents_right[idx]
            node_ids.append(node_id)
            inequality_symbols.append(1)
            thresholds.append(threshold[node_id])
            features.append(feature[node_id])
        path_info[i] = {'node_id': node_ids,
                        'inequality_symbol': inequality_symbols,
                        'threshold': thresholds,
                        'feature': features}

各パスに対する処理の中(for文の中)で,ノードが左側親ノードを持つ場合と右側親ノードをもつ場合に分けて考えています. idxの初期値は葉ノードで,その親ノードに関する情報からnode_ids, inequality_symbols, thresholds, featuresに格納されていきます. そして,最終的にディクショナリpath_infoに葉ノードをkeyとしたディクショナリが格納されます.

一旦,ここまでの処理を関数にします.

def search_path(estimator, class_labels, aim_label):
    """
    return path index list containing [{leaf node id, inequality symbol, threshold, feature index}].
    estimator: decision tree
    maxj: the number of selected leaf nodes
    """
    """ select leaf nodes whose outcome is aim_label """
    children_left = estimator.tree_.children_left  # information of left child node
    children_right = estimator.tree_.children_right
    feature = estimator.tree_.feature
    threshold = estimator.tree_.threshold
    # leaf nodes ID
    leaf_nodes = np.where(children_left == -1)[0]
    # outcomes of leaf nodes
    leaf_values = estimator.tree_.value[leaf_nodes].reshape(len(leaf_nodes), len(class_labels))
    # select the leaf nodes whose outcome is aim_label
    leaf_nodes = np.where(leaf_values[:, aim_label] != 0)[0]
    """ search the path to the selected leaf node """
    paths = {}
    for leaf_node in leaf_nodes:
        """ correspond leaf node to left and right parents """
        child_node = leaf_node
        parent_node = -100  # initialize
        parents_left = []  # 左側親ノード
        parents_right = []  # 右側親ノード
        while (parent_node != 0):
            if (np.where(children_left == child_node)[0].shape == (0, )):
                parent_left = -1  # 左側親ノードが存在しない場合は-1
                parent_right = np.where(
                    children_right == child_node)[0][0]
                parent_node = parent_right
            elif (np.where(children_right == child_node)[0].shape == (0, )):
                parent_right = -1  # 右側親ノードが存在しない場合は-1
                parent_left = np.where(children_left == child_node)[0][0]
                parent_node = parent_left
            parents_left.append(parent_left)
            parents_right.append(parent_right)
            """ for next step """
            child_node = parent_node
        # nodes dictionary containing left parents and right parents
        paths[leaf_node] = (parents_left, parents_right)
        
    path_info = {}
    for i in paths:
        node_ids = []  # node ids used in the current node
        # inequality symbols used in the current node
        inequality_symbols = []
        thresholds = []  # thretholds used in the current node
        features = []  # features used in the current node
        parents_left, parents_right = paths[i]
        for idx in range(len(parents_left)):
            if (parents_left[idx] != -1):
                """ the child node is the left child of the parent """
                node_id = parents_left[idx]  # node id
                node_ids.append(node_id)
                inequality_symbols.append(0)
                thresholds.append(threshold[node_id])
                features.append(feature[node_id])
            elif (parents_right[idx] != -1):
                """ the child node is the right child of the parent """
                node_id = parents_right[idx]
                node_ids.append(node_id)
                inequality_symbols.append(1)
                thresholds.append(threshold[node_id])
                features.append(feature[node_id])
            path_info[i] = {'node_id': node_ids,
                            'inequality_symbol': inequality_symbols,
                            'threshold': thresholds,
                            'feature': features}
    return path_info

次に,$\varepsilon$-satisfactory instance を計算する部分を作ります. これはそんなに難しくないので一気に書いて関数としてまとめておきます.

def esatisfactory_instance(x, epsilon, path_info):
    """
    return the epsilon satisfactory instance of x.
    """
    esatisfactory = copy.deepcopy(x)
    for i in range(len(path_info['feature'])):
        # feature index
        feature_idx = path_info['feature'][i]
        # threshold used in the current node
        threshold_value = path_info['threshold'][i]
        # inequality symbol
        inequality_symbol = path_info['inequality_symbol'][i]
        if inequality_symbol == 0:
            esatisfactory[feature_idx] = threshold_value - epsilon
        elif inequality_symbol == 1:
            esatisfactory[feature_idx] = threshold_value + epsilon
        else:
            print('something wrong')
    return esatisfactory

最後に,上で定めた2つの関数を使って提案手法を実装します.

def feature_tweaking(ensemble_classifier, x, class_labels, aim_label, epsilon, cost_func):
    """
    This function return the active feature tweaking vector.
    x: feature vector
    class_labels: list containing the all class labels
    aim_label: the label which we want to transform the label of x to
    """
    """ initialize """
    x_out = copy.deepcopy(x)  # initialize output
    delta_mini = 10**3  # initialize cost
    for estimator in ensemble_classifier:
        if (ensemble_classifier.predict(x.reshape(1, -1)) == estimator.predict(x.reshape(1, -1))
            and estimator.predict(x.reshape(1, -1) != aim_label)):
            paths_info = search_path(estimator, class_labels, aim_label)
            for key in paths_info:
                """ generate epsilon-satisfactory instance """
                path_info = paths_info[key]
                es_instance = esatisfactory_instance(x, epsilon, path_info)
                if estimator.predict(es_instance.reshape(1, -1)) == aim_label:
                    if cost_func(x, es_instance) < delta_mini:
                        x_out = es_instance
                        delta_mini = cost_func(x, es_instance)
            else:
                continue
    return x_out

まとめ

よく「ブラックボックス」と言われる機械学習ですが,今回紹介した手法を使うと (どうすればpositiveなラベル方向へ改善出来るのかという事がわかるという意味で)解釈性を上げる事ができます. また,Decision Treeのアンサンブル手法であるRandom ForestやGradient Boostingは非常に手軽に試せる手法で, かつ比較的高精度な分類器を用意に作ることが出来ます.そのような分類器において解釈性の向上が得られるということは, 実際のデータマイニング応用にとても役立つと思われます.

【麻雀】配牌から和了を判定する【天鳳牌譜解析】

みなさんこんにちは.接点QBです.今回は天鳳の牌譜を解析して,配牌から和了出来るか否かを判定する分類器を作ってみたいと思います.

僕は学部の頃から麻雀をしていて,数理統計学の研究室で機械学習の研究をしていました.そして,今はコンピュータ・サイエンス領域の研究室に在籍しています. そうです.天鳳の牌譜解析にうってつけのバックグラウンドですねw. しかしながら,牌譜解析のプログラムを組むのが意外と面倒で,手を付けていませんでした. 麻雀関連のプログラムってあまりコードが公開されていなくて,公開されていてもCとJavaがほとんどなので,ちょっと僕には敷居が高かったわけです(PythonとRしか使ったことないので). まあ,せっかく就活も終わったので,時間がある時に牌譜解析やってみようかなーと思いまして,手始めにタイトルのようなことをやってみました. 決して,天鳳で八段からの降段が間近だから打つのを日和って時間が出来たからという理由ではありません!

今回の分類器を作ろうと思った理由ですが,学部時代にサークルの知人が「オレは配牌を見て和了れそうにない時は,最初から安牌を大量に抱える!」と言っていまして, 僕としては「上がれるかどうかなんてツモ次第で分からんやろ.みすみす和了の機会を逃すようなことは良くないのでは?」と考えていまして, 「機械学習で配牌から和了出来るかを判定出来るのか?」と気になったからです. まあ,その知人の脳内での処理が今回作成した分類器に匹敵または勝るという保証は全く有りませが…

天鳳の牌譜解析って結構めんどくさくて,「車輪の再発明」を避けるために出来るだけコードも載せていこうと思います. ただし,個々で書くコードが最適でない可能性が結構高いと思うので,その点はご了承下さいm( )m

牌譜ファイルの形式

牌譜の所得

まず天鳳の牌譜なんですが,日本プロ麻雀協会所属近藤千雄さんがご自身のブログで公開して下さっています(麻雀にどんよくです:天鳳牌譜解析をはじめたい人へ - livedoor Blog(ブログ)).

で,ここから適当な年の牌譜をダウンロードしてきて解凍しますと,txtファイルがあると思います.それが天鳳の牌譜です. 僕が解析するときはUTF8に変換してから使っています.

牌譜の形式

牌譜の形式を説明します.とりあえず2015年の牌譜の先頭20行ほどを見てみましょうか.

===== 天鳳 L0000 鳳東喰赤速 開始 2015/12/31 http://tenhou.net/0/?log=2015123123gm-00e1-0000-f2529b5b&tw=0 =====
  持点25000 [1]雷神魔理沙/七段/女 R2082 [2]天照/七段/男 R2199 [3]心の旋律/七段/男 R2176 [4]<>/七段/男 R2069
  東1局 0本場(リーチ0)  雷神魔理沙 -6000 天照 -3000 心の旋律 13000 <> -3000
    跳満ツモ 立直1 門前清自摸和1 平和1 ドラ2 赤ドラ1 裏ドラ0
    [1東]2m6m9m7p3s8s東南西西北発発
    [2南]3m5m7m7m2p4p5p6p7p6s南北中
    [3西]1m3m4p6p8p9p9p3s3s4s7s7s9s
    [4北]4m7m9m2p2p4p3s4s7s南西白発
    [表ドラ]2s [裏ドラ]5m
    * 1G2s 1d9m 2G発 2d中 3G中 3D中 4G6s 4d西 1G東 1d2m 2G1p 2d発 1N発発 1d北
    * 2G1p 2d北 3G8s 3d7s 4G9s 4D9s 1G7s 1d7p 2G9m 2d南 3G1m 3D1m 4G6m 4d9m
    * 1G4m 1D4m 2C3m5m 2d1p 3G1s 3D1s 4G西 4D西 1N西西 1d6m 2G1p 2D1p 3G8m 3D8m
    * 4G4m 4d発 1G3m 1D3m 2G9p 2d9m 3G1s 3d3m 4G5s 4d白 1G5s 1D5s 2G2m 2D2m
    * 3G4s 3d1m 4G6s 4d南 1G7m 1D7m 2N7m7m 2d1p 3G5p 3d8p 4G1p 4D1p 1G6p 1D6p
    * 2G6s 2d9p 3G2s 3R 3d1s 4G2p 4d7s 1G7p 1d南 2G5M 2D5M 3G中 3D中 4G8p
    * 4D8p 1G3p 1d7s 2G4m 2D4m 3G5S 3A

  東2局 0本場(リーチ0)  心の旋律 4900 <> -3900
    30符3900点3飜ロン 対々和2 断幺九1
    [1北]3m4m6m7m8m2p6p7s8s南西白発

大体何が書いてあるかは分かりますよね? この牌譜から必要な情報を抽出するプログラムを組みます.

これからやること

配牌をone-hot-vectorに変換して,上がれたか否か(0 or 1)をラベル付けします. one-hot-vectorって何?という人のために説明しておきますと,たとえば手牌が「123456m234p78s北北」だとすると,one-hot-vectorは $[1,1,1,1,1,1,0,1,1,1,0,0,\cdots,01,1,0,0,\cdots, 1,1,0,\cdots, 0]$となります.これは下の表に対応しています.

1m 2m 3m 4m 5m 6m 7m 1p 2p 3p 4p 5p 6s 7s 8s
1 1 1 1 1 1 0 0 1 1 1 0 0 1 1 0 2 0

このベクトルに対して,和了出来ていたら0,出来ていなかったら1をラベル付けします.

配牌と和了情報を取り出す

ここからは実際に牌譜を処理するコードを書きます. まずは正規表現を書いておきます.

end_pattern1 = re.compile("  ---- 試合結果 ----")
end_pattern2 = re.compile(" *[1-4]位")
end_pattern3 = re.compile('----- 終了 -----')

pattern_start = re.compile('=====*')  # ゲームスタート行の先頭文字列パターン
pattern_player = re.compile(' *持点')  # プレイヤー情報行の先頭文字列パターン
pattern_kyoku = re.compile(' *[東南西][1-4]局')  # 局情報行の先頭文字列パターン
pattern_tehai = re.compile(' *\[[1-4][東南西北]\]')  # 手牌情報行の先頭文字列パターン
pattern_dahai = re.compile(' *\*')  # 打牌情報行の先頭文字列パターン

次に,配牌と局ごとの和了者を取り出します.

game_id = 0
kyoku_id = 0
kyoku_lines = []  # 局情報
tehai_lines = []  # 手配
dahai_lines = []  # 打牌行格納用リスト
winner = []  # 和了者

file = '../data/houton2015_utf8.txt'
f = open(file, errors='ignore')
# line = f.readline()
for line in f:
    if end_pattern1.match(line) or end_pattern2.match(line) or end_pattern3.match(line):
        ''' 試合結果行の処理 '''
        continue
    elif pattern_start.match(line) != None:
        ''' 試合スタート行の処理 '''
        game_id += 1
    ryukyoku_flag = 0  # 流局フラグ初期化
    while line != '\n':
        '''
        1局分の処理
        '''
        if pattern_kyoku.match(line) != None:
            ''' 局スタート行の処理 '''
            kyoku_id += 1  # whileループを抜けるまで固定される
            tmp = line.strip()  # 行頭のスペース削除
            kyoku_lines.append(tmp)

        elif "流局" in line:
            ''' 流局フラグ '''
            ryukyoku_flag = 1

        elif pattern_tehai.match(line) != None:
            ''' 手牌情報行の処理 '''
            tmp = line.strip()  # 行頭のスペース削除
            tmp = re.sub('\[[1-4][東西南北]\]', '', tmp)
            # tehai_lines.append([kyoku_id, tmp])
            tehai_lines.append(tmp)

        elif pattern_dahai.match(line) != None:
            '''
            打牌情報行の処理
            スペース削除作業をしてリストに追加
            '''
            if ryukyoku_flag == 1:
                '''
                流局フラグが立っていたら、打牌情報から和了者を特定しなくてよい
                '''
                result = 'Ryukyoku'
            else:
                tmp = line.strip()  # 行頭のスペース削除
                tmp = tmp.replace("* ", "")  # 行頭のアスタリスクとスペース削除
                # print(tmp)
                if "A" in tmp:
                    result = tmp[-2]  # while を抜けるまで固定
        line = f.readline()
    winner.append([kyoku_id, result])

これでwinneerに局idと和了結果(流局の場合はRyukyoku,誰かが和了した場合は和了プレイヤーの番号)が格納され, tehai_linesには全ての手牌が格納されています. これらのリストから

| P1の配牌 | P2の配牌 | P3の配牌 | P4の配牌| 和了者 |

というデータフレームを作成します.

""" 和了者処理 """
winner_df = pd.DataFrame(winner, columns=['kyoku_id', 'winner'])
winner_df = winner_df.drop_duplicates('kyoku_id')
winner_df = winner_df.set_index('kyoku_id')
winner_df.columns = ['winner']
""" 手牌処理 """
tehai_df = pd.DataFrame([tehai_lines[i:i+4] for i in range(0, len(tehai_lines), 4)])
tehai_df.index = winner_df.index
tehai_df.columns = ['p1', 'p2', 'p3', 'p4']

result = pd.concat([tehai_df, winner_df], axis=1)

これでresultに手牌と和了者のデータフレームを格納出来ました. 次は,各列に入っている手牌をone-hot-vecrtorに変換し,分類器に学習させるためのデータを作成します. 具体的には(配牌のone-hot-vector, 和了判定)という形式のデータを作成します.

def convert_multi(result_part):
    data_y = []
    data_x = result_part[['p1', 'p2', 'p3', 'p4']].applymap(lambda x: mj.convert_tehai(x))
    for i in result_part['winner']:
        if i=='Ryukyoku':
            y_tmp = list(np.zeros(4))
            data_y.append(y_tmp)
        else:
            y_tmp = mj.one_hot_list(int(i)-4, 4)
            data_y.append(y_tmp)
    x_arr = np.vstack(
        [np.array(
            [list(j["p1"].values), list(j["p2"].values),
             list(j["p3"].values), list(j["p4"].values)]) 
         for i,j in data_x.iterrows()])
    y_arr = np.hstack(data_y)
    return (x_arr, y_arr)

processn = 20
step = int(len(result.index)/processn)
dfs = [result[i: i+step] for i in range(0, len(result.index), step)]
p = Pool(processn)  # make process for multi-processing
df = p.map(convert_multi, dfs)
p.close()
x_arr = np.vstack([i for i,j in df])
y_arr = np.hstack([j for i,j in df])

これで,x_arrに配牌のone-hot-vectorが格納され,y_arrに対応する和了判定が格納されした. ちなみに,上記コードの中にあるprocessnは自分の環境に合った数値に設定してください(並行処理のためのプロセス数です).

それでは,いよいよ分類器を作成します.

"""
split dataset
"""
np.random.seed(10)
sampler = np.random.permutation(len(y_arr))
test_size = 100000
train_index = sampler[2*test_size:]
vali_index = sampler[test_size: 2*test_size]
test_index = sampler[:test_size]
''' training set '''
x_train = x_arr[train_index, :]
y_train = y_arr[train_index]
''' validation set '''
x_vali = x_arr[vali_index, :]
y_vali = y_arr[vali_index]
''' test set '''
x_test = x_arr[test_index, :]
y_test = y_arr[test_index]

"""
Random Forest
"""
rfc = RandomForestClassifier(n_jobs=processn)
rfc.fit(x_train, y_train)
print(rfc.score(x_test, y_test))

結果

2015年の鳳東のデータを使った所,スコアは0.777でした. ハイパーパラメータのチューニングは特に何もしていないので,ちゃんとチューニングすればもう少し精度が良くなると思います. また,分類器を他のものに変えても精度が向上するかもしれません. 今回はRandom Forestを使いましたが,往々にしてGradient Boosting Treeの方が性能が良いことが多いです.

さて,confusion matrixを見てみると次のようになりました.

0 1
0 77323 1192
1 21155 330

やはりFalse negativeが多い(本来は和了できているのに,出来ないと予測しているケースが多い)ですねぇ… データの割合的に仕方ない気もしますが,どうなんでしょう. また,True negativeが多い事も目を引きます.和了れない配牌に対して「和了れない」と正しく判定出来るケースが多く, そういった意味では,和了れないと予測された場合に安牌を抱えておくのは割と悪くないのかなぁ*1とも思います.

*1:このような戦略を取った時の局収支がどうなるかを考えないといけないので,今回の解析結果だけからは何とも言えない

数学を捨てない統計学入門:2.確率変数

前回は抽象的な確率空間$(\Omega, \mathscr{B}, P)$を考えました. しかし,これはまだ実社会での「統計」に直結しそうにありません. そこで,今回は抽象的な確率空間と現実世界を結びつけるための「確率変数」と「分布関数」を導入します.

Borel Algebra

確率空間を現実の空間に結びつけるとき,$\Omega$は$\mathbb{R}$に対応します. そして,\sigma-algebra $\mathscr{B}$に対応するのがボレル集合族$\mathbb{B}$です.

Definition 全ての半開区間$(a,b]$を含む最小の\sigma-algebraを$\mathbb{R}$上のボレル集合族といい$\mathbb{B}$で表す. また,$\mathbb{B}$の元をボレル集合という.

さて,この定義に出てくる「全ての半開区間$(a,b]$を含む最小の\sigma-algebra」は本当に存在するのか? という疑問が(数学をやってる人なら)当然生じます. この疑問に答える命題を証明することが出来ます.

Lemma 標本空間$\Omega$の任意の部分集合族$\mathscr{F}$に対して, $\mathscr{F}$を含む最小の\sigma-algebraが存在する.

上で定義した$\mathbb{B}$は$(a,b]$の形の区間以外にも色々な区間を含みます. たとえば,$\{a\}, (a,b), [a, b], [a, b), (a, \infty), [a, \infty), (-\infty, a]$もボレル集合族に含まれます(つまりボレル集合).

確率変数

何やらややこしげなボレル集合を定義しましたが,これは確率変数(random variable; r.v.)を定義するためです.

Definition \begin{gather} X:\Omega\to\mathbb{R}\text{が}(\Omega, \mathbb{B})\text{上の確率変数} \stackrel{\textrm{def}}{\Longleftrightarrow} \forall B\in \mathbb{B},\, \left\{ \omega | X(\omega)\in B \right\}=X^{-1}(B)\in \mathscr{B} \end{gather}

確率変数は上記のようにして定義されます. 僕はよく確率変数の説明をする時に「取りうる値に対して確率が定まっている変数」という風に説明するのですが, 厳密な定義は上記のようになります. ただ,厳密な定義を知っていて何か役に立ったのか?と聞かれると正直微妙ですw. 僕自身は今は機械学習を専門にしていますが,役に立った記憶はありません. ただし,確率解析をやっていた数学科の後輩は必要不可欠な感じでした.

確率変数の定義は \begin{gather} \forall B\in \mathbb{B},\, \left\{ \omega | X(\omega)\leq a \right\}=X^{-1}(-\infty, a])\in \mathscr{B} \end{gather} と同値になります. これを使って確率変数のイメージを説明すると,次の図のようになります.

f:id:setten-QB:20170805174413p:plain

「実数空間$\mathbb{R}$で$a$以下である点の集合を$X$で$\Omega$に引き戻すとボレル集合になっている」 という図です. ちなみに,可測空間$(\mathbb{R},\mathbb{B})$において$f:\mathbb{R}\to\mathbb{B}$が \begin{gather} \forall a\in \mathbb{R},\,\left\{ x | f(x)\leq a \right\} \in \mathbb{B} \end{gather} を満たすとき,$f$を$(\mathbb{R},\mathbb{B})$上の可測関数といいます. したがって,確率変数は$(\Omega, \mathbb{B})$から$\mathbb{R}$への可測関数なわけですね.

分布関数と確率変数

ここまでで$(\Omega, \mathscr{B})$の実数空間版として$(\mathbb{R}, \mathbb{B})$が導入されました. 残りの確率測度$P$に対応するものとして分布関数$F$を導入します. 以下では,$P(X^{-1}(B))=P(\{ \omega | X(\omega)\in B\})$のことを$P(X\in B)$と表すことにします.

Definition $X$を$(\Omega, \mathscr{B}, P)$上の確率変数とする. このとき, \begin{gather} F:\mathbb{R}\to\mathbb{R}\text{が}X\text{の確率分布} \stackrel{\textrm{def}}{\Longleftrightarrow} F(x)=P(X\leq x) =P\left( \{ \omega | X(\omega)\leq x \} \right) , \, \forall x\in \mathbb{R}. \end{gather}

このように分布関数を定義したとき, \begin{gather} P(a<X\leq b)=P\left( X^{-1}((a,b])\right)=F(b)-F(a),\, \forall x\in \mathbb{R} \end{gather} が成り立ちます.

さて,ここまでの議論で確率空間での$P$に対応するものとして分布関数$F$を導入できました. こうやって$(\mathbb{R}, \mathscr{B}, F)$を導入することで以下のような利点があります.

  • 実数空間で「確率」を定義しようと思ったら,$(a,b]$に対して$F$を定義すれば十分
  • $X$は本当は関数なのですが,$P(X\in B)$と書くことで「$X$が$B$に含まれる確率」と$X$が変数であるかのように扱うことが出来る
  • 多少大雑把な表現になりますが,$(\mathbb{R}, \mathscr{B}, F)$を導入することで,もとの抽象的な確率空間に立ち戻って議論をする必要もなくなる

また,「確率変数$X$がどのような値をどのような確率でとるのか」を,その値と確率を同時に考えて$X$の確率分布あるいは単に分布といいます. もちろん分布関数$F$はそのような情報を与えてくれる関数で,$F(x)$は確率変数$X$が$x$以下の値を取る確率を表します. ちなみに,分布関数について以下のような性質があります.

  1. $F$は単調非減少
  2. $F$は右連続
  3. $F(\infty)=\lim_{t\to\infty}F(t)=1$
  4. $F(-\infty)=\lim_{t\to-\infty}F(t)=0$

確率変数は取りうる値の集合によって,離散型と連続型に分けらます(ルベーグ積分を知っていれば,特に分ける必要はないです). 取りうる値の集合を$E$とすると,$|E|=\aleph_0$なら$X$は離散型,$|E|=\aleph$なら$X$は連続型です. つまり,$X$が自然数・整数・有理数を値に取るなら離散型,無理数・実数に値を取るなら連続型ということです.

離散型確率分布

$X$が離散型確率変数であるとき, $E=\{x_1, \cdots, x_n,\cdots \},\, p(x_i)=P(X=x_i)$とすると \begin{gather} \sum_{i=1}^{\infty}p(x_i)=1,\, F(x)=\sum_{x_i\leq x}p(x_i) \end{gather} が成り立ちます.この$p(\cdot)$のことを$X$の確率質量関数 (p.m.f.) と呼びます.

二項分布

$X$のp.m.f.が \begin{gather} p(x)=\binom{n}{x}p^{x}(1-p)^{n-x},\, x=0,1,\cdots, n; 0<p<1 \end{gather} であるとき,この分布を二項分布といい$\textrm{Bin}(n, p)$で表します. 特に,$n=1$の時の2項分布$\textrm{B}(1,p)$をベルヌーイ分布といいます.

ポアソン分布

$X$のp.m.f.が \begin{gather} p(x)=P(X=x)=\frac{\lambda^x}{x!}\exp(-\lambda),\, \lambda>0, \, x=0,1,\cdots \end{gather} であるとき,この分布をポアソン分布といい$\textrm{Poisson}(n,p)$で表します. 実際には,一定の時間間隔内での機器の故障回数や交通事故数など,稀に起こる事象の発生回数がポアソン分布に従う事が知られています.

連続型確率分布

Definition $X$を$(\Omega, \mathscr{B}, P)$上の確率変数,$F$を$X$の分布関数とすると, \begin{gather} X\text{が連続型確率変数}\stackrel{\textrm{def}}{\Longleftrightarrow} \exists f:\mathbb{R}\to\mathbb{R}^{+}(\text{非負値可測関数})\quad \textrm{s.t.} \quad F(x)=\int_{-\infty}^{x}f(u)\, du,\, \forall x\in \mathbb{R}. \end{gather}

上記定義にあるような$F$を絶対連続な分布関数と呼び,$f$を$X$の確率密度関数(p.m.f.)と言います. そして \begin{gather} f(x)\geq 0,\, \int_{-\infty}^{\infty}f(u)\,du=1,\, \frac{d}{dx}F(x)=f(x) \end{gather} という性質が成り立ちます.

一様分布

$X$のp.d.f.が \begin{gather} f(x;a,b)=\begin{cases} \frac{1}{b-a} & x\in (a,b) \\ 0 & \textrm{otherwise} \end{cases} \end{gather} であるとき,$X$は$(a,b)$上の一様分布に従う言い,$X\sim \textrm{U}(a,b)$と表します($\sim$は「従う」の意味).

指数分布

$X$のp.d.f.が \begin{gather} f(x;\theta)=\begin{cases} \theta \exp(-\theta x) & x\geq 0\\ 0 & x<0 \end{cases} \end{gather} であるとき,$X$は指数分布$\textrm{Ex}(\theta)$に従うといいます.

正規分布

$X$のp.d.f.が

\displaystyle
f\left(x;\mu, \sigma^2\right)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left[-\frac{(x-\mu)^2}{2\sigma^2}\right]

であるとき,$X$は正規分布[tex:\mathcal{N}(\mu, \sigma2)]に従うといいます. 特に$\mathcal{N}(0, 1)$を標準正規分布といい,標準正規分布の密度関数を$\phi$,分布関数を$\Phi$と表すことが多いです.

正規分布はとにかく重要な分布です. 何かしらの例として出てくるのは大体が正規分布ですし,統計学の理論は殆どが何かしらが正規分布に従うという仮定を立てて理論を組み立てて,それを拡張していきます. 線形回帰モデルでも誤差に標準正規分布仮定して理論を組み立てます. 正規分布の形を知ってる人は結構多いと思いますが,標準正規分布を載せておきます.

f:id:setten-QB:20170806000850p:plain

意外と平ぺったいですね.

変数変換された確率変数

確率変数を可測関数で変換してできた変数は,また確率変数になります. このことは取り敢えず認めることとして,変換された確率変数がどのような確率変数となるかを見てみましょう.

実数$a>0,b$と確率変数$X$に対して,$Y=aX+b$も確率変数になります. このとき$Y$の分布を求めてみます.$Y,X$の分布関数をそれぞれ$F_Y, F_X$と書くことにします. $y\in \mathbb{R}$として,定義に忠実に計算を進めましょう. \begin{gather} F_Y(y)=P(Y\leq y)=P(aX+b\leq y)=P(aX\leq y-b) \end{gather} となります.よって, \begin{gather} F_Y(y)=P\left( X\leq \frac{y=b}{a}\right) =F_X\left( \frac{y-b}{a}\right) \end{gather} と求まります. さらに,$X$が連続型であると仮定してp.d.f.を求めます. \begin{gather} f_Y(y)=\frac{d}{dy}F_(y) =\frac{d}{dy}F_X\left( \frac{y-b}{a}\right) =\frac{1}{a}f_Y\left( \frac{y-b}{a}\right) \end{gather} $a<0$の場合に分布関数は \begin{gather} F_Y(y)=P\left( X\geq \frac{y-b}{a}\right) =1-P\left( X< \frac{y-b}{a}\right) =1-P\left( X\leq \frac{y-b}{a}\right) =1-F_X\left( \frac{y-b}{a}\right) \end{gather} となるので,p.d.f.は \begin{gather} f_Y(y)=\frac{d}{dy}F_Y(y) =\frac{d}{dy}\left\{ 1-F_X\left( \frac{y-b}{a}\right)\right\} =-\frac{1}{a}f_X\left( \frac{y-b}{a}\right) \end{gather} となります.したがって,一般に実数$a,b$に対して$Y=aX+b$と変換された確率変数のp.d.f.は \begin{gather} f_Y(y)=\frac{1}{|a|}f_X\left( \frac{y-b}{a}\right) \end{gather} となります.

正規分布に従う確率変数の変換

X\sim \mathcal{N}(\mu, \sigma^2) のとき,Y=(X-\mu)/\sigmaの分布を求めてみましょう. $y\in\mathbb{R}$に対して,先述の計算から

\displaystyle
P(Y\leq y)= P\left( \frac{X-\mu}{\sigma}\leq y\right)
= P(X\leq \sigma y+\mu)
=\int_{-\infty}^{\sigma y+\mu} \frac{1}{\sqrt{2\pi \sigma^2}}\exp\left[
-\frac{(x-\mu)^2}{2\sigma^2}
\right]\, dx
=\int_{-\infty}^{y}\frac{1}{\sqrt{2\pi}}\exp\left( -\frac{t^2}{2}\right)\,dt
=\Phi(y)

となります.ただし,t=(x-\mu)/\sigmaです. このことから,正規分布に従う確率変数は上記のような変換をすれば,標準正規分布に従う事が分かります. この変換は結構大事で,実際にデータを分析する時の前処理などで行います. もちろん,$\mu$や\sigma$には推定値(推定値については推定のところで説明します)を放り込みます.

他にも,「正規分布に従う確率変数の変換によって得られた確率変数」が従う分布には名前がつけられているものが多いです. たとえば,

  • $Y=\exp(X)$が従う分布は対数正規分布
  • $Y=\sum_{i=1}^{n}X_i^2$が従う分布は自由度$n$のカイ2乗分布

などです.

まとめ

今回は確率変数を定義して,確率変数が従う「確率分布」を紹介しました. 次回はこの確率変数を多次元に拡張する話を書こうと思います.