ともにゃん的データ分析ブログ

勉強したことの備忘録とかね

離散選択理論入門(その1)

はじめに

最近、大学院時代に読んでいた

Discrete Choice Methods with Simulation

Discrete Choice Methods with Simulation

を読み直しています。やっぱり面白いです。でも忘れてることも多いです。

ということで簡単にまとめるついでに、離散選択理論入門というシリーズで記事を書いていこうと思います。


離散選択理論とは

人々は日々意思決定をしています。意思決定にも、決定するものの種類によって色々考えられますが、身近な例としては物を買うという行為も意思決定のひとつです。
物を買うという行為は、いくつかの選択肢の中(商品群)の中から買うのもを選びます。このとき、選択肢は離散的な、つまりカテゴリカルなものとなっています。このように、離散的な選択肢の中からあるものを選択するという行為を離散選択と呼びます。そして、なぜそれを選択したのか?それを分析する道具として離散選択理論があります。そこで用いられる代表的なモデルとしては、ロジスティック回帰モデルがよく知られています。

効用と選択の原理

人々は何に基づいて選ぶという行為をしているでしょうか。直感で?好みで?色々と要因が考えられますが、その要因を定量化したものとして、効用関数というものを考えてみます。あることの効用とは、人がそのあることをした結果得られる満足感のことで、特にここではそれを定量的に表したもののことを指します。例えばある人が牛乳をスーパーに行って、明治の牛乳と雪印の牛乳を発見したとします。このとき、もし明治の牛乳を選んだのならば、その人にとっては明治の牛乳を選んだほうが雪印の牛乳を選んだときよりも満足感が高いからそうしたのだと考えます。これを効用を使った表現すると


\begin{eqnarray*}
\mbox{明治の牛乳を選んだときの効用} > \mbox{雪印の牛乳を選んだときの効用}
\end{eqnarray*}

であると言えます。


離散選択理論では、人々の選択の背後にはこのような効用が存在しており、それぞれの選択肢を選んだときの効用を考えながら、その中で最も効用が高い選択肢を選ぶという原理を採用しています。これは、僕たちが実際に脳内で効用を計算して、それによって選択という行為をしていると言っているのではなく、あくまでそういう"モデル"を考えるということです。ただこういったモデル化が、人々の選択という行動の本質をうまく抽出することの助けとなります。


さてこの効用ですが、様々な要素の関数となっていると考えられます。例えば、上記の牛乳の例であれば、"価格"という要因の関数になっていると考えられます。明治の牛乳を選択したとき、それは雪印の牛乳と比べて価格が安かったから選んだのかもしれません。つまり、価格が安いがために、明治の牛乳を選んだときの効用が雪印の牛乳を選んだときの効用を上回ったのかもしれません。他にも上記の例における選択という行為に影響を与える要因としては、牛乳の味や賞味期限の長さといったことも考えられるでしょう。したがって、効用というのはこういった選択肢に関連する様々な要因の関数になっていると考えられます。したがって、今後は効用のことを効用関数と呼んだりします。また効用についての記号を導入し、 n さんの選択肢 i に対する効用を U_{ni} という記号で表したりすることにします。上記の牛乳の例ではそれぞれの選択肢について


\begin{eqnarray*}
i &=& 1:\ \mbox{明治の牛乳を選んだとき}\\
i &=& 2:\ \mbox{雪印の牛乳を選んだとき}
\end{eqnarray*}

とすると


\begin{eqnarray*}
U_{n1} > U_{n2}
\end{eqnarray*}
と表せます。

効用の定式化

効用についての説明しました。しかし、これは選択という意思決定者の脳内のモデル化であって、選択という行為を分析したい人が意思決定者本人でない限り、その人の効用完全に知ることはできません。"外から"人々の選択という行為分析するためにはどうすればいいでしょうか。

上記の牛乳の例に戻りましょう。明治の牛乳を購入したのは、その人にとってその味が好きだから選んだのかもしれません。このとき、分析者はその人の味の好みを知りませんので、明治の牛乳の方が好ましい度、つまり効用がどういったものなのか知ることはできません。

では価格に注目するとどうでしょうか。明治の牛乳を購入したのは単に雪印の牛乳よりも価格が安かったからかもしれません。このとき価格は分析者も分かりますので、外から人の選択行為を分析する要因として使えそうです。

そこで、効用を外から観測できる部分と、観測できない部分に分割することを考えます。 n さんの選択肢 i に対する外から観測できる効用を V_{ni} 、観測できない効用を \epsilon_{ni} とすると、分割された効用は以下のように表せます:


\begin{eqnarray*}
U_{ni} = V_{ni} + \epsilon_{ni}
\end{eqnarray*}

ここで観測できる効用 V_{ni} は、上記の例でいうと、牛乳の価格の関数になっていると考えられます。価格等の要因に関わらない n さんにとっての牛乳 i そのもの好みを \alpha_{ni} 、牛乳 i の価格を p_{i} とすれば、例えば V_{ni} は次のような次のような形を考えることができます:


\begin{eqnarray*}
V_{ni} = \alpha_{ni} + \beta_{n}p_{i}
\end{eqnarray*}

ここで \beta_{n} は、牛乳 i の価格が n さんの牛乳 i への効用にどれだけ影響を与えるかの係数(重み)、言い換えると価格にどれほど反応するかを表していると考えられます。もちろん価格以外にも、もし購入日から賞味期限までの日数などがデータとして手に入っていれば、"賞味期限への係数 x 賞味期限の長さ"といった形を上記の式に組み入れることも可能です。また、上記の V_{ni} はデータと係数の線形和となっていますが、非線形な関係を仮定しても良いでしょう。そして \epsilon_{ni}n さんの味の好みといったような、観測できない要因の効果が全て詰め込まれた、分析者にとっては未知の値です。

このように定式化すると、 n さんの牛乳の購買情報(データ)がいっぱい手に入れば、 n さんの牛乳 i への好み \alpha_{ni} や、 n さんが価格にどれだけ反応するかを表す \beta_{n} といった値が購買情報から推定できそうです。実際に推定できるのですが、その方法は後ほどということで、ひとまずは


\begin{eqnarray*}
U_{ni} = V_{ni} + \epsilon_{ni}
\end{eqnarray*}

という効用の観測できる部分(分析者が得られるデータから推定できそうな部分) V_{ni} と、そもそもデータが手に入らなくて外からだと分からない部分 \epsilon_{ni} に効用を分割する、という一般的な定式化の下で選択確率という概念を導入したいと思います。

選択確率

効用関数を定式化しました。ここで牛乳の例に戻ると、ある人にとって明治の牛乳のほうが雪印の牛乳のほうが、"ある時"は好ましい、つまり U_{n1} > U_{n2} だったのでした。しかし、効用関数を U_{ni} = V_{ni} + \epsilon_{ni} と定式化したため、未観測の値である \epsilon_{ni} の値によっては雪印のほうが好ましい場合もあるでしょう。たとえば、目の前の人が雪印の牛乳を手にとっていたので、なんとなく明治の牛乳よりも雪印の牛乳のほうが"良いように見えて"、そのときは雪印の牛乳のほうが効用が高くなったということもあるでしょう。こういった"目の前の人が雪印の牛乳を手にとっていた"というが与える効用への影響は(ほぼ間違いなく)データとして観測することができないので、効用の未観測の部分である \epsilon_{ni} に押し込まれることになります。つまり、 \epsilon_{ni} はその時々で値を変えるが分析者にとっては不明の確率変数であると考えることができます。

いつも明治の牛乳のほうが雪印の牛乳のほうが好ましいわけではなくて、雪印の牛乳のほうが好ましいときもあるということであれば、確率で議論をするのが良さそうです。つまり、どのくらいの確率で] n] さんは明治の牛乳を好むのか(雪印の牛乳を好まないのか)、という議論をしたほうが良さそうです。そこで出てくる概念が、選択確率というものです。ある選択肢を選ぶ確率のことです。

選択確率を導出しましょう。選択確率は、商品 i よりも商品 j の効用が高いときに商品] i] を選択するという原理を確率で表現したものです。つまり、商品 i よりも商品 j の効用が高い確率が選択確率です。 n さんが商品 i を選ぶ確率、選択確率を P_{ni} とすると、以下のように表現できます:


\begin{eqnarray*}
P_{ni} = \mbox{Prob}(U_{ni} > U_{nj}\ \forall j \neq i)
\end{eqnarray*}

ここで \forall j \neq i は、すべての商品 i でない商品 j に対して、という意味です。この式を更に書き下していきます:


\begin{eqnarray*}
P_{ni} &=& \mbox{Prob}(U_{ni} > U_{nj}\ \forall j \neq i)\\
&=& \mbox{Prob}(V_{ni} + \epsilon_{ni} > V_{nj} + \epsilon_{nj}\ \forall j \neq i)\\
&=& \mbox{Prob}(V_{ni} - V_{nj} > \epsilon_{nj} - \epsilon_{ni}\ \forall j \neq i)\\
&=&  \mbox{Prob}(\epsilon_{nj} - \epsilon_{ni} < V_{ni} - V_{nj}\ \forall j \neq i)
\end{eqnarray*}

ここで、 \epsilon_{n} = (\epsilon_{1},\ldots,\epsilon_{J})^{'} という J 次元ベクトルで、 J 個ある選択肢それぞれに対する効用の未観測部分をベクトル化したものです。また I(\cdot) は指示関数で、カッコ内の条件が満たされていれば 1 、そうでなければ 0 を返す関数です。

以上が選択肢 i の選択確率です。最後の式を見ると、選択確率は未観測な要因 \epsilon_{j},\ \epsilon_{i} の差よりも、観測できる効用 V_{ni},\ V_{nj} の差の方が大きい確率と言い換えることができます。

選択確率が定義できました。しかし、上記のままでは具体的に確率を計算することはできなそうです。そこで。効用関数の未観測部分 \epsilon_{ni} に様々な確率部分布を仮定することで、上記の選択確率を具体的に計算できる形にします。 \epsilon_{ni} にどんな確率分布を仮定すると、上記の確率はどのように計算できて、その確率はどういった性質を持っていて、そこから何が分かるのか、これを調べるのが離散選択理論です。

ここで、 \epsilon_{n} = (\epsilon_{1},\ldots,\epsilon_{J})^{'} という J 次元ベクトルで、 J 個ある選択肢それぞれに対する効用の未観測部分をベクトル化したものとします。効用関数の未観測部分 \epsilon_{ni} に確率分布 f(\epsilon_{n}) を仮定したとすると、上記の選択確率は次のように書けます:


\begin{eqnarray*}
P_{ni} &=&  \mbox{Prob}(\epsilon_{nj} - \epsilon_{ni} < V_{ni} - V_{nj}\ \forall j \neq i)\\
&=& \int_{\epsilon}I(\epsilon_{nj} - \epsilon_{ni} < V_{ni} - V_{nj}\ \forall j \neq i)f(\epsilon_{n})d\epsilon_{n}
\end{eqnarray*}

例えば \epsilon_{ni} にガンベル分布を仮定すると、上記の積分が解析的に解けて、選択確率はみなさんがよく知っている(多項)ロジスティック回帰における確率になります。様々な f((\epsilon_{n}) を仮定することで、個性的な性質をもった様々な選択確率を考えることができ、その性質と人々が選択するという行為を照らし合わせながら、どの確率を用いることで現実の現象をうまく表現できるのか、といったことが検討できます。

次回以降の記事では、 f(\epsilon_{n})パラメトリックな分布を色々と仮定することで導出される様々な選択確率を紹介していきたいと思います。

なぜベイズ予測分布はあれなのか?

お久しぶりです

久しぶりの投稿となります。いつの間にかXGBoostについての記事

kefism.hatenablog.com

が好評を得ていて、ブログのアクセス数も随分と増えていました!非常に嬉しいことです☆(ゝω・)v

今回は、ベイズ関連の教科書を読むと必ず紹介されているアレ(予測分布)がなぜそういう形をしているのかを説明します☆(ゝω・)v

参考文献は

です☆(ゝω・)v

不等式の準備

一般的に、任意の実数 z に対して以下が成り立ちます:

 \log z \leq z -1

したがって、Q(x)P(x) をそれぞれ密度関数としたときに、上記 zP(x)/Q(x) を代入することによって


\begin{eqnarray}
&& \log \frac{P(x)}{Q(x)} \leq \frac{P(x)}{Q(x)} - 1 \\
&\iff& Q(x)\log \frac{P(x)}{Q(x)} \leq P(x) - Q(x)
\end{eqnarray}

を得ることができます。2行目は両辺に  Q(x)を掛けました。

さらに上記を両辺 x について積分をして整理すると


\begin{eqnarray}
&&\int Q(x)\log \frac{P(x)}{Q(x)} dx \leq \int P(x) dx - \int Q(x) dx \\
&\iff& \int Q(x)\log \frac{P(x)}{Q(x)} dx \leq 0 \\
&\iff& \int Q(x)\log P(x) dx - \int Q(x)\log Q(x) dx \leq 0 \\
&\iff& \int Q(x)\log P(x) dx \leq \int Q(x)\log Q(x) dx
\end{eqnarray}

となります。2行目の右辺は、密度関数の積分は1になることを利用しました。3行目の左辺は分数の対数の性質を利用しました。

上記の不等式をみると、密度関数 Q(x) における2つの対数尤度関数 \log P(x)\log Q(x) の期待値を比較したときに、必ず自分自身  Q(x) の対数尤度関数の期待値の方が任意の密度関数  P(x) の対数尤度関数 \log P(x) よりも大きくなることを意味しています。

汎化損失

 x を現在得ているデータ、 q(x) をデータ  x の真の分布、p^{\ast}(x) を我々が構築した予測分布とします。

SlideShareに上げているスライド

www.slideshare.net

でも説明していますが、汎化損失 G とは

G = -\int q(x)\log p^{\ast}(x)dx

で計算される数値のことをいいます。

この汎化損失は、推定された予測分布  p^{\ast}(x) の対数尤度を"真のデータの分布" q(x) で期待値を取ったものにマイナスをかけた形をしています。したがって"真のデータの分布" q(x)がもし分かっている場合に、汎化損失は値が小さければ小さいほど推定された予測分布がどのくらい望ましいかを表す数値であるとみることができます。


さて、はじめに準備した不等式によると、密度関数  q(x) に対する対数尤度の期待値は、自分自身の対数尤度  \log q(w) の場合が最も高くなるのでした。つまり、汎化誤差  G

 -\int q(x)\log q(x)dx \leq -\int q(x)\log p^{\ast}(x)dx = G

となります。データから推定された予測分布よりも、"真の分布"を予測分布としたときに、その予測分布が最も望ましいのは当然のことでしょう。


しかし、真のデータ分布は分かりません。パラメータの事後分布が構築できたときに、何を予測分布とするのが良いのでしょうか?それがどのベイズの教科書を読んでも出てくる、おなじみの予測分布の式

 p(y|x) = \int p(y|\theta)p(\theta|x)d\theta

です。ここで  y は予測値、\theta はパラメータ、p(y|\theta)\theta をパラメータとしてもつ y の分布、p(\theta|x) はデータ x の下で得られたパラメータ \theta の事後分布です(以降同様)。


 p(y|x) = \int p(y|\theta)p(\theta|x)d\theta が良い予測分布である理由

 p(y|x) = \int p(y|\theta)p(\theta|x)d\thetaが良い予測分布であるということですが、それを導出しましょう。

まず、どんな  p(y|x) が良いのかの指標となる、 p(y|x) の対数尤度を以下のように用意します:

 \int p(y|\theta) \log p(y|x) dy

この対数尤度を最大にする  p(y|x) が良い予測分布であるとします。

ここで、データ x は様々な値を取りうるので、上記の対数尤度を x の分布 p(x|\theta) で期待値を取ります:

 \int p(x|\theta) \int p(y|\theta) \log p(y|x) dy dx

さらに、ここでかなり強い仮定、"真のパラメータ  \theta の事前分布 p(\theta) が分かっている"を置きます。このときにこの真の事前分布で上記の式の期待値を取り、式変形を行っていきます:

 
\begin{eqnarray}
&& \int p(\theta) \int p(x|\theta) \int p(y|\theta) \log p(y|x) dy dx d\theta \\ 
&=& \int \int \int p(x|\theta) p(\theta) p(y|\theta) \log p(y|x) dy dx d\theta \\
&=& \int \int \int p(x, \theta) p(y|\theta) \log p(y|x) dy dx d\theta \\
&=& \int \int \int p(\theta|x) p(x) p(y|\theta) \log p(y|x) dy dx d\theta \\
&=& \int p(x) \int \log p(y|x) \int p(y|\theta)p(\theta|x) d\theta dy dx\\
\end{eqnarray}

1つ目の等号は被積分関数をすべて積分記号の中に入れました。1つ目の等号から2つ目の等号、2つ目の等号から3つ目の等号は、よく知られた条件付き確率の式 p(x|\theta) p(\theta) = p(x, \theta) = p(\theta|x)p(x) を使用しました。4つ目の等号は積分の順序を入れ替えて整理しました。

上記の2つ目以降の積分の部分

 \int \log p(y|x) \int p(y|\theta)p(\theta|x) d\theta dy

を確認してください。これは最初に用意した不等号で比較されている値

 \int Q(x)\log P(x) dx \leq \int Q(x)\log Q(x) dx

とよく似た形をしています。したがって

 p(y|x) = \int p(y|\theta)p(\theta|x) d\theta

としてやれば、他の任意の予測分布 p^{'}(y|x) に対して、最初に用意した不等号より

 \int p(y|x) \log p^{'}(y|x) dy \leq \int \log p(y|x) \int p(y|\theta)p(\theta|x) d\theta dy = \int p(y|x) \log p(y|x) dy

となるわけです。以上の不等式の両辺に-1をかけて不等号を入れ替えて、小さければ小さい程よい値という、最初に紹介した汎化損失の形で不等式を表現できます:

  - \int p(y|x) \log p(y|x) dy \leq - \int p(y|x) \log p^{'}(y|x) dy


ところで、 p(x) についての期待値計算を無視していましたが、この不等式は両辺を  p(x) で期待値をとっても不等号の向きは変わりませんのでここでは無視をして大丈夫です。

以上により、事後分布  p(\theta|x) が求まった下で汎化損失を最も小さくする予測分布は

 p(y|x) = \int p(y|\theta)p(\theta|x) d\theta

ということになるのです。


あの強い仮定はそのままで大丈夫?

"真のパラメータの事前分布 p(\theta) が分かっている"という仮定を置いて、この分布についての期待値をとりながら上記の結論を導きました。これはとても強い仮定です。なぜなら、誰も真のパラメータの事前分布なんて知らないからです。

どうやったって分かりっこないのだから、そこは目を瞑るしかありません。なので自分が設定した事前分布を信じて事後分布を構築し、事後分布を構築したならば

 p(y|x) = \int p(y|\theta)p(\theta|x) d\theta

によって予測を行うしかないのです。

ただ、これが汎化損失の意味で"手元にある事後分布"における最良の予測分布なのであって、この予測分布そのものが良い予測分布なのかは分かりません。あくまで、事後分布  p(\theta|x) を得たらならば、あれこれ色々な予測分布を考えるのではなく、上記の予測分布で予測をすればいいと言っているだけです。

そこで良い予測分布はどれなのか?の疑問に答えるのが、渡辺澄夫先生が開発されたWAICです。WAICは"真のデータの分布"と構築した予測分布とのある意味での"近さ"を表している値です。したがって、いくつか事後分布を構築してWAICによってモデル選択を行い、上記の予測分布によって予測するのが良いのではないでしょうか。

Pythonのpandasで大きなデータを扱うときにメモリ効率を上げる方法メモ

最近Pythonを使って大きなデータフレームを結合したり、そのデータフレームに対してメソッドを使って処理をしていますが、頻繁に

"Memory Error"

に遭遇しています。そんな中で色々探し回って見つけたTips的なものをメモとして残しておきます。

  • pandas.DataFrame.merge をする時に copy=False のオプションをつける

そもそもこのcopyってなんだかよくわからないのですが、使用メモリの挙動を見ている限りは copy=True とするとデータフレームを結合するときにそのcopyを用意して、それを参照している感じがします。copy=Falseで劇的にメモリ効率が良くなった気がします。

ちなみにpandas.DataFrame.mergeよりもpandas.DataFrame.joinの方が結合の効率が良いという報告もあります
参照:
python - Improve Pandas Merge performance - Stack Overflow

  • del 変数名 で変数を消して gc.collect() でメモリを開放

参照:
Pythonで少なくメモリを使用する方法 - のんびりしているエンジニアの日記

  • メソッドを呼び出すときは inplace=True オプションをつける

参照:
pandas でメモリに乗らない 大容量ファイルを上手に扱う - StatsFragments


特に中間テーブルなどを結合した最終的な結果データフレームだけを使いたい時などは、途中過程で生成されるデータフレームは結局使用されないと思うので、上記のオプションをつけても問題無いと思います。そうでなくても問題ないのかもしれません。

pandasのDataFrameのplotメソッドを使って描画するときにプロットを並べて表示する方法

pandasのDataFrameのplotメソッドを使って描画するときに、プロットを並べて表示する方法のメモ

import matplotlib.pyplot as plt
import pandas as pd

fig, axes = plt.subplots(nrows=1, ncols=2, squeeze=False)

df1.plot(kind="bar", x="x", y="y", ax=axes[0,0], figsize=(15, 5))
df2.plot(kind="bar", x="x", y="y", ax=axes[0,1], figsize=(15, 5))

XGBoostの概要

XGBoostの凄さに最近気がついたので、もうちょっと詳しく知りたいと思って以下の論文を読みました。

XGBoost: A Scalable Tree Boosting System

せっかくなので、簡単にまとめてみたいと思います。。。と思っていたら結構な量になってしいました。
何か間違い等がありましたらコメントをしていただくか、@kefism へ連絡をしてくださると嬉しいです。

続きを読む