【論文紹介】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]$という特徴量ベクトルを分類する事を考えます． このとき，１つ目のノードの条件が$x_1\leq \theta_1$，２つ目のノードでの条件が$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}$に対応します．

図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_leftとchildren_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':
        '''
        １局分の処理
        '''
        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:このような戦略を取った時の局収支がどうなるかを考えないといけないので，今回の解析結果だけからは何とも言えない

数学を捨てない統計学入門：２．確率変数

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

Borel Algebra

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

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

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

Lemma 標本空間$\Omega$の任意の部分集合族$\mathscr{F}$に対して， $\mathscr{F}$を含む最小の-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} と同値になります． これを使って確率変数のイメージを説明すると，次の図のようになります．

「実数空間$\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$の時の２項分布$\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.が

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

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

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

変数変換された確率変数

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

実数$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} となります．

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

のとき，の分布を求めてみましょう． $y\in\mathbb{R}$に対して，先述の計算から

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

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

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

などです．

まとめ

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

数学を捨てない統計学入門 １．確率

私は学部が数学系で大学院が情報系なのですが，情報系に移ってきた時に「一部の人達は統計学の数学的側面をあまり理解してないな」と感じました． ということで，何回かに分けて「数学を捨てない統計学入門」と題して記事を書こうかと思います． 内容としては，学部レベルの数理統計学で扱う内容です． 「数理統計学」というタイトルが付いた本を見てもらえれば，大体同じ内容が載っているのでそちらを見ていただいても良いと思います．

確率の導入

「数学を捨てない」と題しているからには，中高で習うような統計的確率ではなくKolmogorovが定めた公理的確率を確率の定義として採用したいと思います． そのために，まずは-algebraを定義します．

-algebra

確率は標本空間（集合）の部分集合族（部分集合の集合）で定義されます． ただし，その部分集合族はどんな集合族でも良いわけではありません． 適当な条件を満たす，の部分集合族の元に対してのみ確率が定義されます． その部分集合族というのが-algebraで，以下の条件を満たすの部分集合族$\mathscr{A}$です．

• 
• 
• 

さて，この-algebra上に測度を用いて確率を定義していくわけですが， 「測度って何だよ？速度の誤変換か？」という人のために測度の大雑把な説明をしておきたいと思います．

測度

測度を一言で説明すると，「集合の大きさを測るための尺度」です． たとえば，「で構成される集合の大きさは？」と聞かれれば， 「一辺の長さが1の正方形なので，」と答えるでしょう． これは，特定の尺度（ここでは正方形の面積）を通して，図形に対して実数値を対応させていることに他なりません． この集合と実数値の対応を測度といいます．

確率測度

さて，上記の-algebraと測度を用いて，やっと確率を定義できます． 以下の3条件を満たすを，可測空間上の確率測度といいいます．

\begin{gather} \forall A\in \mathscr{A},\, 0\leq P(A)\leq 1 \tag{P1} \\ P(\Omega)=1 \tag{P2} \\ {\displaystyle A_i\in \mathscr{A}, A_i\cap A_j=\emptyset\,(i\neq j)\Longrightarrow P\left( \bigcup_{i=1}^{\infty}A_i\right) =\sum_{i=1}^{\infty}P(A_i)} \tag{P3} \\ \end{gather}

1つ目の条件は，皆さんご存知の「確率は1を超えない」というやつですね． 2つ目は，全事象の確率は1という事を表していて， 3つ目は，排反な事象の和事象の確率はそれぞれの事象の確率の和になるという事を表しています．

ちなみに，の3つの組を確率空間と言います． この確率空間上で，以下のことが成り立ちます．

\begin{align} & P(\emptyset)=0 \\ & A_i \in \mathscr{A}\,(i=1,\cdots, n), A_i\cap A_j=\emptyset\,(i\neq j) \Longrightarrow P\left( \bigcup_{i=1}^{n}A_i\right)=\sum_{i=1}^{n}P(A_i)\\ & A\in \mathscr{A}\Longrightarrow P(A)+P(A^{c})=1 \\ & A,B\in \mathscr{A}\Longrightarrow P(A\cup B)=P(A)+P(B)-P(A\cup B) \\ & A_i \in \mathscr{B}\,(i=1,\cdots, n)\Longrightarrow P\left( \bigcup_{i=1}^{n}A_i\right) \leq \sum_{i=1}^{n}P(A_i) \\ & A_1\subset A_2\subset \cdots \subset,\, A_i\in \mathscr{A}\Longrightarrow P\left( \bigcup_{n=1}^{\infty}A_n\right)=\lim_{n\to \infty}P(A_n) \\ & A_1\supset A_2\supset \cdots \supset,\, A_i\in \mathscr{A}\Longrightarrow P\left( \bigcap_{n=1}^{\infty}A_n\right)=\lim_{n\to \infty}P(A_n) \\ \end{align}

いずれも「確率なんだから成り立って当然でしょ？」と思われる命題ですが， 数学では「定義→定理」という流れを大切にします． 先に定めた確率測度の定義のみから，上記の命題たちが成り立つことがわかります．

条件付き確率

確率空間$(\Omega, \mathscr{A}, P)$において，$B\in \mathscr{A}$が$P(B)>0$を満たすとする． このとき，$\forall A\in \mathscr{A}$, \begin{equation} P(A|B)=\frac{P(A\cap B)}{P(B)} \tag{1.1} \end{equation} を，事象$B$が与えられたときの$A$の条件付き確率といいます． そして，$P(\cdot\,|B)$は確率測度の定義の3条件(P1), (P2), (P3)を満たしていることが確認できます． 実際に証明してみましょう．

まずは(P1)の成立を確認します． $P(\cdot)$は確率測度なので，$\forall A\in\mathscr{A},\, P(A\cap B)\geq 0$です． また，仮定から$P(B)>0$なので，$P(A|B)\geq 0$となります． 次に(P2)の成立を確認します． $$ P(\Omega |B)=\frac{P(\Omega \cap B)}{P(B)} =\frac{P(B)}{P(B)} = 1 $$ となり，(P2)を満たすことがわかります． 最後に(P3)の成立を確認します．$A_i\in \mathscr{A}, A_i\cap A_j=\emptyset\,(i\neq j)$とすると \begin{align} P\left( \left. \bigcup_{i=1}^{\infty} A_i \right| B\right) =& \frac{P\left( \bigcup_{i=1}^{\infty} A_i \cap B\right) }{P(B)}\\ =& \frac{P(\bigcup_{i=1}^{\infty} (A_i\cap B))}{P(B)} \\ =& \frac{\sum_{i=1}^{\infty}P(A_i\cap B)}{P(B)}\\ =& \sum_{i=1}^{\infty}\frac{P(A_i\cap B)}{P(B)}\\ =& \sum_{i=1}^{\infty}P(A_i|B) \end{align} となり，(P3)を満たすことが分かります． 以上のことから，$P(\cdot\, |\,B)$は$(\Omega, \mathscr{A}, P)$上の確率測度であることが示されました．

事象の独立性

$A, B\in\mathscr{A}$に対して \begin{equation} P(A\cap B)=P(A)P(B)\tag{1.2} \end{equation} が成り立つとき，$A$と$B$は独立であると言います．

高校までは，独立なら(1.2)が成り立つと教わったと思いますが， 本来は逆です．(1.2)が成り立つときに独立という概念が定義されます．

今回のまとめ

• 確率は特別な集合族に属する集合の大きさを測るための測度

次回は確率変数や分布の話を書こうと思います．

テスト

$a+b=c$