めもめも

このブログに記載の内容は個人の見解であり、必ずしも所属組織の立場、戦略、意見を代表するものではありません。

ロジスティック回帰のROC曲線を描くコード

ロジスティック回帰とは?

ロジスティック回帰は線形判別法の一種で、Type=0, Type=1 の2種類のカテゴリーに分類されるデータがあるときに、「あるデータが Type=1 に属する確率」を推定します。

たとえば、(x, y) 平面上に Type=1, Type=0 のサンプルが10個づつあるとして、次のような結果が得られます。確率(probability)が高い順にソートしてあり、確率が0.5以上/以下の堺に線を入れてあります。

            x          y  type  probability
0   25.100600  15.215185     1     0.972372
1   25.716642  10.509214     1     0.957894
2   26.328260   9.368392     1     0.955292
3   20.965230   6.646373     1     0.906886
4   16.940683   7.185182     1     0.876261
5   15.384313   4.033190     0     0.814820
6   -4.292937  19.269243     1     0.776546
7   19.822191  -3.042080     1     0.760062
8   12.615442  -4.103175     0     0.591393
9   15.223528  -6.948369     1     0.577949
--------
10  10.978713  -7.110957     1     0.475399
11  26.402198 -21.755863     0     0.450887
12  -7.920640  -3.024727     0     0.192805
13 -19.024978   6.163114     0     0.181953
14   1.152526 -14.663486     1     0.140983
15   6.027717 -20.457424     0     0.123408
16 -16.228075  -3.017577     0     0.099239
17 -28.237444   2.175471     0     0.058384
18  -6.533585 -19.910959     0     0.044175
19  -4.426988 -22.366264     0     0.041661

この結果は、次のように図示されます。

●が Type=1 で、×が Type=0 を表します。青い直線は Type=1 の確率が 0.5 の境界を表しており、この線より右上に進むと確率が上がって行き、左下に進むと確率が下がっていきます。このように単純な直線で分類するので「線形判別法」と呼ばれます。

このグラフは既存のサンプルデータ(トレーニングセット)を元に境界線を引きましたが、この後、新たな(どちらのタイプか分からない)未知のデータが来たときは、この境界線からの距離で Type=1 の確率が計算されます。ただし、図から分かるように、「確率 0.5 以上なら Type=1」という判断は必ずしも正しいとは限りません。この例では、20%のデータは誤って判定されています。

判定基準の決定

それでは、現実の問題にロジスティック回帰を適用する場合、どの程度の確率を境目に「Type=1」と判定するのがよいのでしょうか? これは、適用する問題の「中身」に依存します。

たとえば、これが、「重大な病気の一次検査」で、「Type=1」はその病気にかかっているものとします。医師の立場としては、わずかでも Type=1 の可能性があれば、念の為、精密検査に回したくなるでしょう。つまり、「確率 X 以上なら Type=1」と判断する X は、なるべく小さな値にしたくなります。(あまり小さくし過ぎると、「精密検査を受けたら実は病気ではなかった」という人が大量発生して、それはそれで困りますが。)

一方、これが、「馬券予想」で、「Type=1」は「当たり馬券」だとします。あなたが予想屋として、「当たり馬券」の情報を販売するならどうするでしょうか? 当然ながら、なるべく確度の高い情報を売るために、「確率 X 以上なら Type=1」と判断する X は、なるべく大きな値にしたくなります。(かといって、あまり大きくし過ぎると、「本当は当たりだったはず」の情報を売りそこねてしまいますが。)

このように、「確率 X 以上なら Type=1」と判断する X は、その大小によって、メリット/デメリットがあり、適用する問題に応じて決定する必要があります。一般的には、このようなメリット/デメリットの関係は、「真陽性(True Positive:TP)」と「偽陽性(False Positive:FP)」の割合で表すことができます。

  • 真陽性率:TP rate = 「Type=1の中で正しくType=1と判定した数」÷「Type=1の総数」
  • 偽陽性率:FP rate = 「Type=0と中で誤ってType=1と判定した数」÷「Type=0の総数」

医者の立場としては、病気の患者を見逃さないために、「TP rate = 1」を目指して X をなるべく小さくしますが、その分だけ「FP rate」も大きくなります。

予想屋の立場としては、はずれ情報で信頼を失わないために、「FP rate = 0」を目指して X をなるべく大きくしますが、その分だけ「TP rate」は小さくなります。

そこで、さまざまな X の値について、サンプルデータから計算した、TP rateとFP rateを縦横の軸にとったグラフを描いてみます。実際の操作としては、先に得られた「確率の高い順にソートしたサンプル」に対して、「上から何個目までをType=1と判定するか」という個数を1つずつ増やしながら、それぞれの場合の TP rate, FP rateを計算して、その点をグラフに打っていきます。このようなグラフを「ROC曲線」と呼びます。

ランダムなデータで数値シュミレーションするコードがあるので、これを実行してみましょう。

# -*- coding: utf-8 -*-
#
# ロジスティック回帰のROC曲線
#
# 2015/04/24 ver1.0
#

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pandas import Series, DataFrame
from numpy.random import rand, multivariate_normal
pd.options.mode.chained_assignment = None

#------------#
# Parameters #
#------------#
Variances = [50,150] # 両クラス共通の分散(2種類の分散で計算を実施)


# データセット {x_n,y_n,type_n} を用意
def prepare_dataset(variance):
    n1 = 80
    n2 = 200
    mu1 = [9,9]
    mu2 = [-3,-3]
    cov1 = np.array([[variance,0],[0,variance]])
    cov2 = np.array([[variance,0],[0,variance]])

    df1 = DataFrame(multivariate_normal(mu1,cov1,n1),columns=['x','y'])
    df1['type'] = 1
    df2 = DataFrame(multivariate_normal(mu2,cov2,n2),columns=['x','y'])
    df2['type'] = 0
    df = pd.concat([df1,df2],ignore_index=True)
    df = df.reindex(np.random.permutation(df.index)).reset_index()
    return df[['x','y','type']]

# ロジスティック回帰を実施
def run_simulation(variance, subplot):
    training_set = prepare_dataset(variance)
    training_set1 = training_set[training_set['type']==1]
    training_set2 = training_set[training_set['type']==0]
    ymin, ymax = training_set.y.min()-5, training_set.y.max()+10
    xmin, xmax = training_set.x.min()-5, training_set.x.max()+10
    subplot.set_ylim([ymin-1, ymax+1])
    subplot.set_xlim([xmin-1, xmax+1])
    # 分類データを表示
    subplot.scatter(training_set1.x, training_set1.y, marker='o')
    subplot.scatter(training_set2.x, training_set2.y, marker='x')

    training_set['bias'] = 1
    w = np.array([[0],[0.1],[0.1]])
    phi = training_set[['bias','x','y']]
    phi = phi.as_matrix(columns=['bias','x','y'])
    t = training_set[['type']]
    t = t.as_matrix()

    # 最大100回のIterationを実施
    for i in range(100):
        # IRLS法によるパラメータの修正
        y = np.array([])
        for line in phi:
            a = np.dot(line, w)
            y = np.append(y, [1.0/(1.0+np.exp(-a))])
        r = np.diag(y*(1-y)) 
        y = y[np.newaxis,:].T
        tmp1 = np.linalg.inv(np.dot(np.dot(phi.T, r),phi))
        tmp2 = np.dot(phi.T, (y-t))
        w_new = w - np.dot(tmp1, tmp2)
        # パラメータの変化が 0.1% 未満になったら終了
        if np.dot((w_new-w).T, (w_new-w)) < 0.001 * np.dot(w.T, w):
            w = w_new
            break
        w = w_new

    # 分類誤差の計算と確率付きデータの用意
    d0, dx, dy = w[0], w[1], w[2]
    err = 0
    training_set['probability'] = 0.0
    for index, line in training_set.iterrows():
        a = np.dot(np.array([1, line.x, line.y]), w)
        p = 1.0/(1.0+np.exp(-a))
        training_set.ix[index, 'probability'] = p
        if (p-0.5)*(line.type*2-1) < 0:
            err += 1
    err_rate = err * 100 / len(training_set)
    # 境界線(P=0.5)を表示
    linex = np.arange(xmin-5, xmax+5)
    liney = - linex * dx / dy - d0 / dy
    label = "ERR %.2f%%" % err_rate
    subplot.plot(linex,liney ,label=label, color='blue')
    subplot.legend(loc=1)
    result = training_set.sort_index(by=['probability'], ascending=[False])
    # 確率付きデータを返却
    return result.reset_index()

# ROC曲線の表示
def draw_roc(result, subplot):
    positives = len(result[result['type']==1])
    negatives = len(result[result['type']==0])
    tp = [0.0] * len(result)
    fp = [0.0] * len(result)
    for index, line in result.iterrows():
        for c in np.arange(0, len(result)):
            if index < c:
                if line.type == 1:
                    tp[c] += 1
                else:
                    fp[c] += 1
    tp_rate = np.array(tp) / positives
    fp_rate = np.array(fp) / negatives

    subplot.set_ylim([0, 1])
    subplot.set_xlim([0, 1])
    subplot.set_xlabel("False positive rate")
    subplot.set_ylabel("True positive rate")
    subplot.plot(fp_rate, tp_rate)

# Main
if __name__ == '__main__':
    fig = plt.figure()
    for c, variance in enumerate(Variances):
        subplots = fig.add_subplot(2,2,c+1)
        result = run_simulation(variance, subplots)
        subplots = fig.add_subplot(2,2,c+2+1)
        draw_roc(result, subplots)
    fig.show()

実行結果は次のようになります。2種類のデータを用いた結果が左右にあります。左は、もともとのデータがうまく分かれていて正しい判定が(ある程度は)容易な例で、右は、データが入り混じっていて正しい判定が困難な例です。

どちらのグラフも「TP rateを上げるとFP rateも上がる」「FP rateを下げるとTP rateも下がる」というメリット/デメリットの交換条件を示しています。

理想の判定は、「TP rate=1、FP rate=0」、すなわち上記グラフの左上の点を実現することですので、一般論としては、なるべく左上に近い点を実現する X を選択するのがよいことになります。左の例は、正しい判定が容易なので左上に近い点があります。右の例は正しい判定が困難なので左上にはあまり近づけていないことが分かります。

ちなみにグラフの右上は、「とにかく全部 Type=1 と言っとけ」という極端な例です。この場合、TP rateだけを見ると「TP rate=1」ですばらしい判定ロジックのように思えますが、実際には何も考えていないことに注意してください。「おれは当たり馬券を全部予想した(TP rate=1)!」という宣伝文句にだまされてはいけません。

同様にグラフの左下は、「とにかく全部 Type=0 と言っとけ」という例です。さらに、左下と右上を結ぶ直線(TP rate = FP rate = Y)上の点は、「根拠はないけど確率 Y で Type=1 と言ってみる(サイコロの出た目で当たり馬券を予想するようなもの)」という判定法にあたります。したがって、あなたが考えだした独自の判定法がこの直線よりも右下に来たとすると、これは、「適当な予想よりもさらに精度の悪い判定法」ということになります。