kanayamaのブログ

twitter: @tkanayama_

ポケモンを題材に因果推論を実践してみる

@tkanayama_です。最近「計量経済学*1」と「効果検証入門 *2」を読んだので、せっかくなので実際に手を動かすことによって理解の整理をしたいと思いました。

www.yuhikaku.co.jp

gihyo.jp

そこで今回は、人工データを用いて「ボールの性能と捕獲確率」の関係性を効果検証してみました(人工データの生成方法は記事の末尾に記述しました)。

問題設定

今は昔、モンスターボールしか存在せず、スーパーボールが世の中で出回り始めたばかりの頃、オーキド博士が「スーパーボールは本当にモンスターボールより捕まえやすいのか?」という仮説を検証しようとしています。

そこでオーキド博士は世界中のトレーナーたちからデータを収集し、下記のような100件のデータを収集することができました。

データID 使ったボール*3 捕まえたポケモン ポケモンを捕まえるのに使ったボールの個数
1 モンスターボール コラッタ 3個
2 スーパーボール イワーク 2個
・・・ ・・・ ・・・
100 スーパーボール フシギダネ 4個

(実際のデータはここに置きました)

オーキド博士はさっそくこれらのデータを用いて、モンスターボール・スーパーボールそれぞれに対し、捕まえるのに使ったボールの個数の平均値を算出しました。

ボール ポケモンを捕まえるのに使ったボールの個数の平均値
モンスターボール 5.22個
スーパーボール 6.74個

予想に反して、スーパーボールのほうがモンスターボールよりも個数が多くなってしまいました。この結果だけ見ると、「スーパーボールはモンスターボールよりポケモンを捕まえにくい」と言えそうです。

有意差検定

トレーナーたちから収集したデータが偶然偏っていたために、スーパーボールの性能が低く出てしまったとオーキド博士は考えました。そこで、今回の結果がどれほど信頼できる値なのか、統計的仮説検定を使って検証することにしました。

モンスターボールを使ってポケモンを捕まえるのに使うボールの個数」を表す確率変数を  Y^{モンスター}、「スーパーボールを使ってポケモンを捕まえるのに使うボールの個数」を表す確率変数を  Y^{スーパー}とおき、  Y^{モンスター} Y^{スーパー}がそれぞれ平均  \mu_{モンスター} および  \mu_{スーパー}・分散  \sigma *4正規分布に独立に従うと仮定します。

帰無仮説を「 \mu_{モンスター} = \mu_{スーパー}」・対立仮説を「 \mu_{モンスター} \lt \mu_{スーパー}」として、有意水準5%で片側検定を行います。 モンスターボールおよびスーパーボールのデータ数をそれぞれ  m, nとすると、母分散の不偏標本分散  s_{モンスター}^2, s_{スーパー}^2はそれぞれ

 s_{モンスター}^2 = \frac{1}{m-1} \Sigma(Y^{モンスター}_i  - \bar{Y}^{モンスター})^2

 s_{スーパー}^2 = \frac{1}{n-1} \Sigma(Y^{スーパー}_i - \bar{Y}^{スーパー})^2

で表されます。ただし、 \bar{Y}^{モンスター}および \bar{Y}^{スーパー}はそれぞれ  Y^{モンスター} Y^{スーパー}の標本平均です。

このとき、

 \frac{( \bar{Y}^{モンスター} -  \bar{Y}^{スーパー}) - (\mu_{モンスター} -  \mu_{スーパー})}{\sqrt{s^2(\frac{1}{m} + \frac{1}{n})}}

は、自由度  m+n-2 のt分布に従います。ただし、

 s^2 = \frac{1}{m+n-2}({(m-1) s_{モンスター}^2 + (n-1) s_{スーパー}^2} )

とします。

帰無仮説のもとでは  \mu_{モンスター} -  \mu_{スーパー} = 0 が成り立つので、結局

 \frac{\bar{Y}^{モンスター} -  \bar{Y}^{スーパー}}{\sqrt{s^2(\frac{1}{m} + \frac{1}{n})}} \tag{1}

と「自由度  m+n-2 のt分布における95%限界値」を比較すれば良いことになります。実際のデータを用いてそれぞれ計算すると、

 (式(1)の値) = -1.17

 (95\%限界値) = -1.66

となります。したがって有意差検定の結果、対立仮説は棄却できず、どちらの性能が良いとも言えないということになりました。

交絡因子の存在

この結果に違和感があったオーキド博士は、データを提供してくれたトレーナーたちに聞き取り調査を行いました。その結果、「捕まえやすそうなポケモンにはモンスターボールを使い、捕まえにくそうなポケモンにはスーパーボールを使っている」ということがわかりました。スーパーボールを節約しながら必要な場面でだけ使いたいトレーナーとしては、納得の行動原理です。

「捕まえやすさ」と「捕まえるのに使ったボールの数」を、ボールの種類別に色分けして散布図を書くと下のようになりました。

f:id:tepppei:20200504225226p:plain
ボールの種類と捕まえやすさ・使ったボールの個数の関係

横軸が「ポケモンの捕まえやすさ」を表す数値で、数値が大きいほど捕まえやすいことを表します。縦軸が、そのポケモンを捕まえるのに使ったボールの個数です。赤の点がモンスターボール・青の点がスーパーボールを表します。これを見ると、たしかに捕まえやすいポケモンほどモンスターボールが使われやすいことがわかります。

この場合、「ポケモンのつかまえやすさ」は、「どのボールを使うか」と「捕まえるのに使ったボールの個数」の両方に対して相関を持っています。このような変数を交絡因子と呼び、交絡因子を無視して効果検証を行うと実際の効果と異なる結果が出てしまいます。直感的にも、今回の問題設定ではスーパーボールが不利になってしまうことがわかると思います。そこでオーキド博士は線形重回帰を用いたモデル化を行い、「ポケモンのつかまえやすさ」の影響を取り除いたうえで、改めて「スーパーボールは本当にモンスターボールより捕まえやすいのか?」という仮説を検証することにしました。

線形重回帰によるモデル化

データ数をNとし、「ポケモンのつかまえやすさ」を表す確率変数を  X_1, ,,, X_N、「どのボールを使ったか(スーパーボールなら1、モンスターボールなら0)」を表す確率変数を  Z_1, ,,, Z_N、「ポケモンを捕まえるのに使ったボールの数」を表す確率変数を Y_1, ,,, Y_Nとします。ここで、「ポケモンのつかまえやすさ」はポケモンごとに一意に定まる値で、その数値はあらかじめわかっているものとします。あらためて、データを示すと下記のようになります。

データID 使ったボールの種類(Z) つかまえたポケモン ポケモンの捕獲しやすさ(X) ポケモンを捕まえるのに使ったボールの個数(Y)
1 0 コラッタ 255 3
2 1 イワーク 45 2
・・・ ・・・ ・・・
100 1 フシギダネ 45 4

 X, Y, Zには次のような関係が成り立っていると仮定します。

 Y = \textbf{X} \beta + \epsilon \tag{2}

ここで、  \textbf{X} *5 N \times 3の行列で、1列目は  (1, ,,, 1)^T・2列目は  (X_1, ,,, X_N)^T・3列目は (Z_1, ,,, Z_N)^Tです。また、 \betaは3次元ベクトル   (\beta_0, \beta_1, \beta_2) とします。また、 \epsilonは誤差項と呼ばれるN次元の確率変数で、 \epsilon  \textbf{X}で条件づけた分布は、平均0・各要素の分散が \sigma^2で各要素が無相関のN次元正規分布に従うと仮定します。このとき、式(2)の両辺に対してXの条件付き期待値をとると、 E[\epsilon| \textbf{X}] = 0 を用いて

 E[Y| \textbf{X}] = \textbf{X} \beta \tag{3}

となることがわかります。つまり、式(2)のモデル化は、Xの条件のもとでのYの期待値を線形モデルとしてモデル化しているという見方ができます。

さて、この仮定のもとで何らかの方法で  \beta X, Yを用いて推定できたとします。この  \betaの推定値を用いて、式(3)からYの期待値を下記のように計算できます。

 E[Y|X, Z=1] = \beta_0 + \beta_1 X + \beta_2 \times 1

 E[Y|X, Z=0] = \beta_0 + \beta_1 X + \beta_2 \times 0

したがって、 E[Y|X, Z=1]  -  E[Y|X, Z=0] = \beta_2が得られ、 \beta_2は「X(ポケモンの捕まえやすさ)の条件を揃えたうえでZ(使うボール)だけ変えたとき、捕まるまでの回数(Y)の平均がどれだけ変化するか?」を表しているということがわかります。

回帰係数の推定

  \textbf{X} \textbf{Y}を使って \betaを推定します。いろいろな推定方法がありえますが、最小二乗法

 \hat{\beta} = argmin_\beta ||Y -   \textbf{X} \beta|| ^ 2 \tag{4}

を用いて推定された \hat{\beta}は、最小分散線形不偏推定量(線形な不偏推定量の中で、分散が最小な推定量)になります。

最小二乗法の解は、式(4)を \betaについて偏微分した値が0になることを用いて計算すると、

 \hat{\beta} = (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T Y \tag{5}

となります。念のため  \hat{\beta}が不偏推定量であることを確認しておくと、

 E[\hat{\beta} | \textbf{X}] =  (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T E[Y|\textbf{X}] = (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T (\textbf{X} \beta + E[\epsilon| \textbf{X}]) = 0

となり、たしかに不偏であることが確認できます*6

ここで実際のデータに式(5)を適用すると、

 \hat{\beta} = (14.60, -0.05, -3.20)^T

となりました。今回興味があるのは  \hat{\beta_2}です。この値が-3.2ということは、「スーパーボールを用いてポケモンを捕まえると、モンスターボールを用いた場合に比べて、捕まえるまでに必要なボールの数が平均で3.2個少なくなる」ということを表しています。

回帰係数の仮説検定

最後に、仮説検定を行うことによって上記で得られた \hat{\beta_2}が統計的にどの程度信頼できるのかを確認します。帰無仮説は「 \beta_2 = 0(つまり、スーパーボールを使ってもつかまえやすくならない)」、対立仮説は「  \hat{\beta_2} \lt 0」とします。

まず、  \hat{\beta_2}が従う分布を調べます。式(5)より、

   \hat{\beta} = (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T Y = (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T (\textbf{X} \beta + \epsilon) = (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T \textbf{X} \beta + (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T \epsilon = \beta +  (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T \epsilon

が成り立ちます。仮定より \epsilonはXの条件下で平均0・分散 \sigma^2正規分布に従うため、その線形結合で表される  \hat{\beta}もXの条件下で正規分布に従います。平均は

  E[\hat{\beta} | \textbf{X}] =  \beta +  (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T  E[\epsilon | \textbf{X}] =  \beta

分散は

  Var[\hat{\beta} | \textbf{X}] =  Var[(\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T \epsilon] =  AVar[ \epsilon ]A^T = \sigma^2 AA^T =  \sigma^2  (\textbf{X}^T \textbf{X})^{-1}

です。なお、途中で A = (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T と置きました。したがって、 (\textbf{X}^T \textbf{X})^{-1} の2行目・2列目の成分を  q^{22} とすると、

 Var[\hat{\beta_2} | \textbf{X}] = \sigma^2 q_{22}

と計算できます。 \hat{\beta_2}正規分布に従うことと、その平均・分散が計算できたので、標準正規分布に従う確率変数

 \frac{\hat{\beta_2} - \beta_2}{\sigma \sqrt{q_{22}}} \tag{6}

を作ることができます。

しかし、真の分散  \sigma^2はわからないので、これを不偏推定量  S^2で置き換えます。ただし、

 S^2 = \frac{||Y - X\hat{\beta}||^2}{n - 3}

です。式(6)の \sigma^2 S^2で置き換えた式

 \frac{\hat{\beta_2} - \beta_2}{S \sqrt{q_{22}}} \tag{6}

は自由度N-3のt分布に従います。

これで、 \hat{\beta_2}が従う分布がわかりました。あとは、帰無仮説のもとでは \beta_2 = 0なので、

 \frac{\hat{\beta_2}}{S \sqrt{q^{22}}} \tag{7}

を計算してt分布の95%限界値と比較することにより、仮説検定を行います。

実際にデータを当てはめると、

 (式(7)の値) = -2.20

 (95\%限界値) = -1.66

であるので、「スーパーボールはモンスターボールより捕まえやすい」という仮説が95%有意水準で採択されることがわかりました。

補足など

残差の分布について

 e = Y -  \textbf{X} \beta = (I -  \textbf{X}( \textbf{X}^T \textbf{X})^{-1} \textbf{X}^T)Y

として、eとX, eとZの関係を散布図にすると下記のようになります。

f:id:tepppei:20200505110916p:plain
eとXの関係

f:id:tepppei:20200505110938p:plain
eとZの関係

この散布図は実現値の一つに過ぎないため、これだけ見ても正確なことは何も言えませんが、なんとなくXの値が大きくなるとeの分散が小さくなっているように見えます。線形重回帰でモデル化を行う際に、「誤差項の各要素の分散が \sigma^2であること」、つまり誤差項の分散が全て等しいことを仮定していましたが、この仮定が満たされていない可能性が考えられます。等分散性が満たされないと、例えば最小二乗法で推定された    \hat{\beta}は分散最小不偏推定量になりませんし、式(6)もt分布に従わないことになります。そのため、より正確な推定・検定を行うためにはeの分散も考慮したモデルを構築する必要があるかもしれません。

他の交絡因子について

今回は交絡因子として「ポケモンの捕まえやすさ」を取り上げましたが、他にも「状態異常」や「どれだけHPを削ったか」も交絡因子になりえます。これらの因子が捕まるまでの回数に影響を与えることはjudge_captureの実装から明らかなので、ボールの選択にも相関を持つ場合に交絡因子になります。今回のシミュレーションでは「状態異常なし」・「HPを削っていない」という条件で揃えたため交絡因子にはなりませんでしたが、実データを集めた場合にはこれらの因子にも考慮する必要がありそうです。

データの生成方法について

下記サイトを参考にしてシミュレーションデータを生成しました。

wiki.xn--rckteqa2e.com

実装はPythonで下記のように行いました。judge_captureポケモンが捕まったかどうかを判定し、count_thrown_ballsでそのポケモンを捕まえるのに何個のボールを使ったかを計算します。

from typing import Optional, Tuple

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


class Ball:
    def __init__(self, ball_type: str) -> None:
        self._random_variable_upper_bounds = dict(
            monster=256,
            super=201,
            hyper=151
        )
        self._ball_coefficient_mapper = dict(
            monster=12,
            super=8,
            hyper=12
        )
        self.ball_type = ball_type

    def get_random_variable(self) -> int:
        random_variable_upper_bound = self._random_variable_upper_bounds[self.ball_type]
        return np.random.randint(0, random_variable_upper_bound)

    def get_coefficient(self) -> int:
        return self._ball_coefficient_mapper[self.ball_type]


class Pokemon:
    def __init__(self, rarity: int, hp: int, hp_max: int, state: Optional[str] = None) -> None:
        self._state_mapper = dict(
            asleep=25,
            freeze=25,
            burn=12,
            paralysis=12,
            poison=12
        )
        self._state = state
        self.rarity = rarity
        self.hp = hp
        self.hp_max = hp_max

    def get_state_adjustment(self) -> int:
        return self._state_mapper[self._state] if self._state else 0


def _calculate_f_value(ball: Ball, pokemon: Pokemon) -> int:
    numerator = int(pokemon.hp_max * 255 / ball.get_coefficient())
    denominator = int(pokemon.hp / 4)
    f_value = numerator / denominator if denominator > 0 else numerator
    return min(255, f_value)


def judge_capture(ball: Ball, pokemon: Pokemon) -> bool:
    if ball.ball_type == 'master':
        return True
    first_random_variable = ball.get_random_variable() - pokemon.get_state_adjustment()
    if first_random_variable < 0:
        return True
    if first_random_variable > pokemon.rarity:
        return False
    second_random_variable = np.random.randint(0, 256)
    if second_random_variable <= _calculate_f_value(ball, pokemon):
        return True
    return False


def count_thrown_balls(ball: Ball, pokemon: Pokemon) -> int:

    def f(n: int) -> int:
        if judge_capture(ball=ball, pokemon=pokemon):
            return n
        return f(n + 1)

    return f(1)

参考文献

*1:計量経済学 (y21) -- 浅野 皙, 中村 二朗

*2:効果検証入門〜正しい比較のための因果推論/計量経済学の基礎 -- 安井 翔太

*3:1匹のポケモンを捕まえるのに1種類のボールしか使っていないとします。

*4:つまり、母分散が等しいことを仮定しています。

*5:わかりにくいですが、太字にして区別しています。

*6:最小分散線形不偏推定量であることを示すためには、他の線形不偏推定量よりも分散が小さいことを別途示す必要があります