RandomForestClassifierの使い方

公式ドキュメント

パラメータ

  • DecisionTreeのアンサンブル学習なので多くはDecisionTreeと同じ。こちらを参照

    特有のパラメータ

  • n_estimators
  • bootstrap
  • oob_score
  • n_jobs
  • verbose
  • warm_start(調査中)

パラメータを変えて様子をみる。

n_estimators

  • 他のパラメータはデフォルトのまま。つまり、汎化能力は低いdecision treeの組み合わになる。
  • AUCの推移 f:id:nsb248:20170224185532p:plain
  • estimatorの数が1のときに比べ、数個のestimatorを加えるだけでかなり精度が上がっている。

bootstrap

  • ツリー構築時に学習データからbootstrapをするかのフラグ。bootstrapをすることで性能が向上する。相関が減るからだっけ?(確認中)
  • AUCの推移
    • 赤がTrue、青がFalse f:id:nsb248:20170224192550p:plain

verbose

  • 途中のツリー構築処理のログを出力してくれる。
  (prop.get_family(), self.defaultFamily[fontext]))
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed:    0.0s finished
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed:    0.0s finished
[Parallel(n_jobs=1)]: Done  10 out of  10 | elapsed:    0.0s finished

DecisionTreeClassifierの限界

これまで

限界を探る。

関係のないデータを混ぜる。

  • x2というデータを入れる。
  • targetはx2に全く依存していない。

    結果

  • accuracy: 0.810
  • std: 0.100
  • 分岐ツリー f:id:nsb248:20170224174722p:plain
  • 重要度
{'x2': 0.24604190914667379, 'x1': 0.45724873019357137, 'x0': 0.29670936065975489}
  • かなり精度が悪化した。
  • 関係ないはずのx2の重要度が0.25もある。

データを45度回転させる。

結果

  • accuracy: 0.835
  • std: 0.037
  • 分類 f:id:nsb248:20170224180953p:plain
  • 分岐ツリー f:id:nsb248:20170224181001p:plain
  • 多少ではあるが精度が悪化している。
  • 分類結果を見ると境界が直線であってほしいにも関わらず、階段状になってしまう。理論上仕方ないが。

DecisionTreeClassifierのパラメータ調整

設定

  • DecisionTreeClassifierの使い方はこちら
  • 下記のデータの分類をDecisionTreeClassifierを使って行う
    • 約10%ほどノイズが入っている。
    • x0: (0.5, 1.5), x1:(0.5, 1.5)で分岐すれば精度90%で分類するごとができる。 f:id:nsb248:20170224165542p:plain
  • 何も制約を入れずに全データを使って学習させると、下記の通りオーバーフィッティングしている。

とりあえずパラメータを変えてみる

  • クロスバリデーションを使ってモデルの評価を行う。

    デフォルトパラメータの場合(ツリーへの制約がない状態)

  • accuracy: 0.800
  • std: 0.039
  • 分類 f:id:nsb248:20170224165616p:plain
  • 分岐ツリー f:id:nsb248:20170224165639p:plain
  • かなり複雑になっている。
  • 90%くらいの精度になってほしいが、80%くらい。

max_depthを調整する。

max_depth=3

  • accuracy: 0.777
  • std: 0.035
  • 分類 f:id:nsb248:20170224170115p:plain
  • 分岐ツリー f:id:nsb248:20170224170128p:plain

max_depth=4

  • accuracy: 0.889
  • std: 0.038
  • 分類 f:id:nsb248:20170224170242p:plain
  • 分岐ツリー f:id:nsb248:20170224170248p:plain

  • max_depth=4でかなりいい感じの結果になった。

GridSearchを使って最適なパラメータを探す。

  • accuracy: 0.882
  • std: 0.044
  • 分類 f:id:nsb248:20170224171535p:plain
  • 分岐ツリー f:id:nsb248:20170224171543p:plain

まとめ

  • ツリーの大きさを制約するようなパラメータを上手く設定することで、汎化性能が上がることを確認。

DecisionTreeClassifierの使い方

公式ドキュメント

パラメータを変えて様子をみる。

サンプルデータ

decision treeで分類しやすいように格子状のデータを作成する。 f:id:nsb248:20170224141608p:plain

パラメータを変えて実験

デフォルトのパラメータのまま

  • 分類結果 f:id:nsb248:20170224141857p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224141902p:plain:w200
  • 重要度
{'x1': 0.69674185463659166, 'x0': 0.30325814536340839}
  • ちゃんとすべて分類できている。線形で分類できるからな。
  • feature importance で本当はx0とx1の重要度は同じはず、ただ、最初の分岐がx0だったので、x0の方が高くなったいる。使用する乱数によって変わる。

criterion=‘entropy’

  • 分類結果 f:id:nsb248:20170224144138p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224144147p:plain:w200
  • 今回のケースだと何もかわらず。。
  • どう結果に影響するかを他のケースを作って検証しよう(また今度)

splitter=‘random’

  • 分類結果 f:id:nsb248:20170224144549p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224144600p:plain:w200
  • 重要度
{'x1': 0.51479236812570128, 'x0': 0.48520763187429866}
  • 境界をランダムに決めているので、分岐ツリーが巨大になっている。
  • 重要度は近い値に。

max_features=1

  • 今、特徴数は2なので、1とした場合をみる。2の場合は良い方で分岐するが、1だと乱数次第。
  • 分類結果 f:id:nsb248:20170224145125p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224145138p:plain:w200
  • 重要度
{'x0': 0.052631578947368363, 'x1': 0.94736842105263164}
  • 複雑な分岐になる。
  • 重要度も乱数次第なので、当てにならない。

max_depth=3

  • ツリーの深さの最大を3とする。3に達したらそこで終了
  • 分類結果 f:id:nsb248:20170224150236p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224150239p:plain:w200

min_samples_split=80

  • nodeに残っているサンプル数が80個未満になったら分岐せずにそこで終了。
  • 分類結果 f:id:nsb248:20170224150033p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224150035p:plain:w200

min_samples_leaf=11

  • リーフに残るサンプル数の最低数を指定
  • 分類結果 f:id:nsb248:20170224151055p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224151058p:plain:w200

min_weight_fraction_leaf = 6 / 400

  • class_weightを指定していないので、min_samples_leaf=6と同じ。
  • 分類結果 f:id:nsb248:20170224151323p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224151325p:plain:w200

max_leaf_nodes=10

  • リーフの最大個数を指定。今回は10個のリーフができている。
  • 分類結果 f:id:nsb248:20170224151457p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224151501p:plain:w200

min_impurity_split=0.48

  • criterionが指定した値以下になったら、そこで分岐終了
  • 分類結果 f:id:nsb248:20170224151817p:plain:w200
  • 分岐ツリー f:id:nsb248:20170224151819p:plain:w200

まとめ

  • 今回は各パラメータの意味を確認できた。
  • キレイに分類できるようにデータを作っている&予測をしていないので、当然、パラメータ調整しても意味がない。

ソースコードの一部

import numpy as np
import pandas as pd
from sklearn import tree
import nsb
import pydotplus


def grid_data():
    n = 20
    x0, x1 = np.meshgrid(np.linspace(0, 2, n), np.linspace(0, 2, n))

    def is_even(x):
        return np.round(x) % 2 == 0

    y = np.logical_xor(is_even(x0), np.logical_not(is_even(x1))).astype(np.int64)
    return pd.DataFrame({'x0': x0.ravel(), 'x1': x1.ravel(), 'y': y.ravel()})


if __name__ == '__main__':
    df = grid_data()
    clf = tree.DecisionTreeClassifier(min_impurity_split=0.48)
    clf.fit(df[['x0', 'x1']], df['y'])
    nsb.plot.scatter_with_boundary(df['x0'], df['x1'], df['y'], clf, 'fig.png')
    dot_data = tree.export_graphviz(clf, out_file=None, feature_names=['x0', 'x1'])
    graph = pydotplus.graph_from_dot_data(dot_data)
    graph.write_png('tree.png')

    print(dict(zip(['x0', 'x1'], clf.feature_importances_)))

To Do, Doing, Done

To Do

Doing

Done

因果推論3 操作変数法

統計的因果推論の勉強

問題設定

  • 目的
    • バウチャー制度の効果を評価したい。
  • 背景
    • くじ引きで当たると、「私立高校へ進学したら授業料が半額になる権利」がもらえる。
    • くじ引きは希望者全員が引ける。当たりは40%とする。
    • 私立高校、公立高校ともに入学試験はないとする。
    • しかし、私立高校の方が質の高い授業を行っているとする。
    • 卒業時に、私立・公立高校で統一テストを実施している。
  • 状況整理
    • 全部で4つのパターンがある
      • くじ引き当たり、私立へ進学する
      • くじ引き当たり、公立へ進学する
      • くじハズれ、私立へ進学する
      • くじハズれ、公立へ進学する
    • バウチャー制度の評価ポイント
      • くじ引きが当たれば、私立へ進学し、外れたら公立へ進学する生徒の卒業時の成績。

観察データ作成(ケース1)

  • 全体の10%当たりは結果に関係なく必ず公立へ進学する。
  • 生徒の20%は裕福で、上記の10%に入った生徒以外はくじ引きの結果に関係なく私立へ進学する。
  • 生徒の成績の平均は、y = 600 + 80 * (私立) + 150 * (裕福)
    • 観測できない影響を大きくするために、150 * (裕福)とする。
  • 標準偏差は50
  • 生徒数は1000
  • 生成したデータはこんな感じ。
  • 生徒が裕福かどうかは不明とする。
 lottery  private       score
0        0        0  605.830138
1        0        0  595.145346
2        0        0  660.159643
3        1        1  676.014558
4        0        0  612.031560
5        0        0  592.204991
6        0        1  704.909145
7        0        0  618.559375
8        0        0  626.251403
9        0        0  536.343088
       lottery      private        score
count  1000.00000  1000.000000  1000.000000
mean      0.41700     0.485000   649.887460
std       0.49331     0.500025    72.562223
min       0.00000     0.000000   476.666647
25%       0.00000     0.000000   595.276742
50%       0.00000     0.000000   644.046202
75%       1.00000     1.000000   704.860884
max       1.00000     1.000000   856.068556

簡単な分析(ケース1)

  • 単純に私立と公立で平均をとってみる。
    • データ生成時のモデル上、差は80程度のはず。
    • 他の要素の影響が含まれている。
106.589641791

操作変数法

  • 割当
    • くじ引きが当たる: z=1
    • くじ引きが外れる: z=0
  • 操作変数
    • z=1のとき、私立へ進学する: d_1=1
    • z=1のとき、公立へ進学する: d_1=0
    • z=0のとき、私立へ進学する: d_0=1
    • z=0のとき、公立へ進学する: d_0=0
    • 結果に影響しないので、d_1 \geq d_0仮定する(単調性)
  • 私立へ進学するかどうか
    • d=z d_1 + (1 - z) d_0
  • テストの成績
    • y=d y_1 + (1 - d) y_2
  • 独立性
    • yz|dは独立
    • d_zzは独立
  • 推定したいもの
    • くじ引きが当たれば私立へ進学するが、外れれば公立へ進学する生徒の成績の差
      • E[y_1 - y_0|d_1=1, d_0=0]
  • 計算する上での問題
    • 条件を満たす生徒が誰なのかがわからない。(当たっても私立へ進学しない生徒の存在など。)
    • 当然、y_1y_0のどちらかしか観測できていない。
  • どう計算するか(詳細な計算は本を見てください)。
    • 準備
      • E[y|z=1-E[y|z=0]=E[y_1-y_0|d_1-d_0=0]p(d_1-d_0=1)]
      • p(d_1=1)=p(d_1=1, d_0=1) + p(d_1=1, d_0=0)
      • p(d_0=1)=p(d_0=1,d_1=1)
      • E(d|z=1)=p(d_1=1)
      • E(d|z=0)=p(d_0=1)
    • 以上より
      • E(y_1 - y_0|d_1=1, d_0=0)=\frac{E(y|z=1)-E(y|z=0)}{p(d_1-d_0=1)}=\frac{E(y|z=1)-E(y|z=0)}{E(d|z=1)-E(d|z=0)}
    • 最終的にはデータから計算できるものばかりで表される!!
  • 結果は期待通り80に近い。
81.2776422538

今後の課題

  • 結果の検定をどうするのか?
  • 入学試験がある場合はどうすればよいか?(選択バイアスが出るよね?)

使ったプログラムコード

import numpy as np
import pandas as pd


def sampling():
    n = 1000
    score_base = 600
    score_private = 80
    score_rich = 150
    score_sigma = 50
    rich_ratio = 0.2
    public_ratio = 0.1
    winning_ticket_ratio = 0.4

    data = pd.DataFrame({
        'rich': np.random.binomial(1, rich_ratio, n),
        'lottery': np.random.binomial(1, winning_ticket_ratio, n)
    })
    data['private'] = data.apply(lambda d: 1 if d['rich'] == 1 or d['lottery'] == 1 else 0, axis=1)
    data['private'] = data['private'] * np.random.binomial(1, 1 - public_ratio, n)
    data['score'] = score_base + score_private * data['private'] + score_rich * data['rich'] \
        + np.random.normal(0, score_sigma, n)

    return data.drop(['rich'], axis=1)


def simple_analysis(data):
    d0 = data[data['private'] == 0]
    d1 = data[data.apply(lambda d: d['private'] == 1 and d['lottery'] == 1, axis=1)]
    print(np.average(d1['score']) - np.average(d0['score']))


def instrumental_variable_method(data):
    data0 = data[data['lottery'] == 0]
    data1 = data[data['lottery'] == 1]
    y0 = np.average(data0['score'])
    y1 = np.average(data1['score'])
    d0 = np.average(data0['private'])
    d1 = np.average(data1['private'])
    print((y1 - y0) / (d1 - d0))


if __name__ == '__main__':
    np.random.seed(123)
    data = sampling()
    print(data.head(10))
    print(data.describe())
    simple_analysis(data)
    instrumental_variable_method(data)
    print('hello')```

因果推論2 傾向スコア:母集団に偏りがあるときの施策の効果を測る

統計的因果推論の勉強

問題

  • AとBの2つの学校があり、Aには何か数学の点数が上がるような施策をし、Bには何もしない。
  • AとBには、親の学歴の分布が異なる。
  • 2つのグループの設定が異なるので、単純には比較できない。
  • 親の学歴とかを除いて、施策が効果があったのかを比較したい。

単純に無作為サンプリングできる場合。

  • 上記のように異なる特徴をもつ集団があったとしても、AとBの生徒をグチャ混ぜにして、無細工に選んで施策を行った場合、単純に比較できる。
  • この場合、施策以外の特徴量を揃えることができるから。
  • 例えば、数学の点数は施策に関係なく、親の学歴に依存している場合を考える。 f:id:nsb248:20170207200206p:plain
  • 赤が施策あり。青が施策なし。
  • 縦軸が数学の点数、横軸が親の学歴。
  • 赤の平均が49.56。青の平均が49.67。
  • 期待通り、今回の施策には効果がないことがわかる。
    n = 1000
    parent_school_career = np.random.normal(3, 1, n)
    treated = np.random.binomial(1, 0.5, n)
    math_score = np.random.normal(50, 10, n) + 10 * (parent_school_career - 3)
    data = pd.DataFrame({'math': math_score, 'parent': parent_school_career, 'treated': treated})

    # plot all sampling data
    plt.scatter(data['parent'], data['math'], c=data['treated'])
    plt.show()

    idx1 = data['treated'] == 1
    idx0 = data['treated'] == 0
    math0 = data[idx0]['math']
    math1 = data[idx1]['math']
    print('z=0:', np.average(math0), np.std(math0))
    print('z=1:', np.average(math1), np.std(math1))

母集団に偏りがある場合

  • Aの親の方がBに比べ、学歴が高いとする。
  • Aに施策を行った。 f:id:nsb248:20170207200905p:plain
  • Aが赤。Bが青。
  • 赤の平均が54.74。青の平均が45.49。
  • 施策には効果がないのに、平均を単純に比較すると有意に差が出てしまう。
    n = 1000
    treated = np.random.binomial(1, 0.5, n)
    parent_school_career = np.random.normal(2.5, 1, n) + treated
    math_score = np.random.normal(50, 10, n) + 10 * (parent_school_career - 3)
    data = pd.DataFrame({'math': math_score, 'parent': parent_school_career, 'treated': treated})

    # plot all sampling data
    plt.scatter(data['parent'], data['math'], c=data['treated'])
    plt.ylabel('math score')
    plt.xlabel("parent's school career")
    plt.show()

    idx1 = data['treated'] == 1
    idx0 = data['treated'] == 0
    math0 = data[idx0]['math']
    math1 = data[idx1]['math']
    print('z=0:', np.average(math0), np.std(math0))
    print('z=1:', np.average(math1), np.std(math1))

傾向スコアを用いて共変量(親の学歴)の影響をなくす。IPW推定を行う。

  • 赤の平均が49.81。青の平均が49.52。
  • 共変量の影響がなくなり、期待通り平均の差がなくなった。
def calculate_score(x, alpha):
    a = np.dot(x, np.array([alpha]).T).flatten()
    s = 1 / (1 + np.exp(-a))
    return s


def loglikelihood(x, z, alpha):
    s = calculate_score(x, alpha)
    p = z * np.log(s) + (1 - z) * np.log(1 - s)
    return np.sum(p)


def non_random_sampling_ipw():
    n = 1000
    treated = np.random.binomial(1, 0.5, n)
    parent_school_career = np.random.normal(2.5, 1, n) + treated
    math_score = np.random.normal(50, 10, n) + 10 * (parent_school_career - 3)
    data = pd.DataFrame({'math': math_score, 'parent': parent_school_career, 'treated': treated})
    data['const'] = np.ones(n)

    idx1 = data['treated'] == 1
    idx0 = data['treated'] == 0

    cons = (
        {'type': 'ineq', 'fun': lambda theta: 10 - theta},
        {'type': 'ineq', 'fun': lambda theta: theta + 10},
    )
    alpha = np.random.uniform(0, 1, 2)
    res = minimize(lambda a: -loglikelihood(data[['parent', 'const']], data['treated'].as_matrix().flatten(), a),
                   alpha,
                   constraints=cons,
                   method="SLSQP")

    score = calculate_score(data[['parent', 'const']], res.x)
    w0 = (1 - data['treated']) / (1 - score)
    w1 = data['treated'] / score
    n0 = np.sum(w0)
    n1 = np.sum(w1)
    y0 = np.sum(w0 * data['math']) / n0
    y1 = np.sum(w1 * data['math']) / n1
    print(y0)
    print(y1)

今後の課題

  • IPW推定を使ったときの分散はどう求める? どうやって検定を行う?
  • 傾向スコアのモデルが適切でないと最終結果がよくない。今回  y = \alpha x + \betaでなく y = \alpha x としていたときはかなり悪かった。
  • 潜在的な共変量があるときはどうすれば?