<jupyter_start><jupyter_text>XGBoost による顧客離反分析 (Churn Analysis)------ コンテンツ1. [背景](1.背景)1. [セットアップ](2.セットアップ)1. [データ](3.データ)1. [学習](4.学習)1. [ホスティング](5.ホスティング)  1. [評価](5-1.評価)  1. [推論エラーのコスト](5-2.推論エラーのコスト)  1. [最適な閾値を探す](5-3.最適な閾値を探す)1. [エンドポイントの削除](6.エンドポイントの削除)--- 1.背景_このノートブックで実施する内容は、[AWS blog post](https://aws.amazon.com/blogs/ai/predicting-customer-churn-with-amazon-machine-learning/)にも記載されています。_どのようなビジネスであっても、顧客を失うことは大きな損害です。もし、満足していない顧客を早期に見つけることができれば、そのような顧客をキープするためのインセンティブを提供できる可能性があるでしょう。このノートブックでは、満足していない顧客を自動で認識するために機械学習 (Machine Learning, ML) を利用する方法を説明します。このような顧客の離反分析は Customer Churn Prediction と呼ばれています。機械学習モデルは完璧な予測を行えないので、このノートブックでは予測のエラーが生じたときの相対的なコストを考慮して、機械学習を利用したときの成果を金額で評価します。ここでは、私達にとってなじみのある離反分析、携帯電話会社からの離反を取り上げます。携帯電話会社が、ある顧客が離反しそうと察知したら、その顧客にタイムリーにインセンティブを与えます。つまり、電話をアップグレードしたり、新しい機能を使えるようになったりして、引き続き携帯電話会社を使おうと思うかもしれません。インセンティブは、顧客が離反して再度獲得するまでにかかるコストよりもずっと小さいことが多いです。--- 2.セットアップまず、このノートブックインスタンスに付与されている IAM role を `get_execution_role()` から取得しましょう。後ほど、SageMaker の学習やホスティングを行いますが、そこで IAM role が必要になります。そこで、ノートブックインスタンスの IAM role を、学習やホスティングでも利用します。通常、role を取得するためにはAWS SDKを利用した数行のコードを書く必要があります。ここでは `get_execution_role()` のみで role を取得可能です。SageMaker Python SDK は、データサイエンティストが機械学習以外のコードを簡潔に済ませるために、このような関数を提供しています。<jupyter_code># bucket = '<your_s3_bucket_name_here>'
# prefix = 'sagemaker/DEMO-xgboost-churn'

# Define IAM role
import boto3
import re
from sagemaker import get_execution_role

role = get_execution_role()<jupyter_output><empty_output><jupyter_text>以降で利用するライブラリをここで読み込んでおきます。<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import io
import os
import sys
import time
import json
from IPython.display import display
from time import strftime, gmtime
import sagemaker

print('Current SageMaker Python SDK Version ={0}'.format(sagemaker.__version__))<jupyter_output><empty_output><jupyter_text>上記セルを実行して、SageMaker Python SDK Version が 1.xx.x の場合、以下のセルのコメントアウトを解除してから実行してください。実行が完了したら、上にあるメニューから [Kernel] -> [Restart] を選択してカーネルを再起動してください。再起動が完了したら、このノートブックの一番上のセルから再度実行してください。その場合、以下のセルを実行する必要はありません。<jupyter_code># !pip install -U --quiet "sagemaker==2.16.1"<jupyter_output><empty_output><jupyter_text>--- 3.データ携帯電話会社は、どの顧客が最終的に離反したか、または、サービスを使い続けたかの履歴データをもっています。この履歴データに対して学習を行うことで、携帯電話会社の顧客離反を予想するモデルを構築します。モデルの学習が終わった後、任意の顧客のデータ (モデルの学習で利用したものと同じ情報を利用します）をモデルに入力すると、モデルはその顧客が離反しそうかどうかを予測します。もちろん、モデルは誤って予測することも考えられるので、将来を予測することはやはり難しいですが、そのような誤りに対応する方法も紹介します。ここで利用するデータセットは一般的に利用可能で、書籍 [Discovering Knowledge in Data](https://www.amazon.com/dp/0470908742/) の中で <NAME> が言及しているものです。そのデータセットは、著者によって University of California Irvine Repository of Machine Learning Datasets に提供されています。ここでは、そのデーセットをダウンロードして読み込んでみます。Jupyter notebook では、冒頭に `!` を入力することで、シェルコマンドを実行することができます。AWS CLIを用いてS3からデータをダウンロードします。<jupyter_code>!aws s3 cp s3://sagemaker-sample-files/datasets/tabular/synthetic/churn.txt ./<jupyter_output><empty_output><jupyter_text>カレントディレクトリ下にダウンロードされた `churn.txt` を `pandas` を利用して読み込んでみます。 `pandas` は、表形式のデータを読み込んで、様々な加工ができるライブラリです。例えば、以下を実行すると表形式でのデータ表示が可能です。<jupyter_code>churn = pd.read_csv('./churn.txt')
pd.set_option('display.max_columns', 500)
churn<jupyter_output><empty_output><jupyter_text>データをみると 5,000 行のデータしかなく、現在の機械学習の状況から見ると、やや小さいデータセットです。各データのレコードは、ある米国の携帯電話会社の顧客のプロフィールを説明する21の属性からなります。その属性というのは、- `State`: 顧客が居住している米国州で、2文字の省略形で記載されます (OHとかNJのように)- `Account Length`: アカウントが利用可能になってからの経過日数- `Area Code`: 顧客の電話番号に対応する3桁のエリアコード- `Phone`: 残りの7桁の電話番号- `Int’l Plan`: 国際電話のプランに加入しているかどうか (yes/no)- `VMail Plan`: Voice mail の機能を利用しているかどうか (yes/no)- `VMail Message`: 1ヶ月の Voice mail のメッセージの平均長- `Day Mins`: 1日に通話した時間(分)の総和- `Day Calls`: 1日に通話した回数の総和- `Day Charge`: 日中の通話にかかった料金- `Eve Mins, Eve Calls, Eve Charge`: 夜間通話にかかった料金- `Night Mins`, `Night Calls`, `Night Charge`: 深夜通話にかかった料金- `Intl Mins`, `Intl Calls`, `Intl Charge`: 国際通話にかかった料金- `CustServ Calls`: カスタマーサービスに電話をかけた回数- `Churn?`: そのサービスから離反したかどうか (true/false)最後の属性 `Churn?` は目的変数として知られ、MLのモデルで予測する属性になります。目的変数は2値 (binary) なので、ここで作成するモデルは2値の予測を行います。これは2値分類といわれます。それではデータを詳しく見てみます。まずはカテゴリデータごとにデータの頻度をみてみます。カテゴリデータは、`State`, `Area code`, `Phone`, `Int’l Plan`, `VMail Plan`, `Churn?`で、カテゴリを表す文字列や数値がデータとして与えられているものです。`pandas`ではある程度自動で、カテゴリデータを認識し、`object`というタイプでデータを保存します。以下では、`object` 形式のデータをとりだして、カテゴリごとの頻度を表示します。また `describe()`を利用すると各属性の統計量を一度に見ることができます。<jupyter_code># Frequency tables for each categorical feature
for column in churn.select_dtypes(include=['object']).columns:
    display(pd.crosstab(index=churn[column], columns='% observations', normalize='columns'))

# Histograms for each numeric features
display(churn.describe())
%matplotlib inline
hist = churn.hist(bins=30, sharey=True, figsize=(10, 10))<jupyter_output><empty_output><jupyter_text>データを見てみると以下のことに気づくと思います。- `State` の各頻度はだいたい一様に分布しています。- `Phone` はすべて同じ数値になっていて手がかりになりそうにありません。この電話番号の最初の3桁はなにか意味がありそうですが、その割当に意味がないのであれば、使うのは止めるべきでしょう- 数値的な特徴量は都合の良い形で分布しており、多くは釣り鐘のようなガウス分布をしています。ただ、`VMail Message`は例外です。- `Area code` は数値データとみなされているようなので、非数値に変換しましょうさて、実際に`Phone`の列を削除して、`Area code`を非数値に変換します。<jupyter_code>churn = churn.drop('Phone', axis=1)
churn['Area Code'] = churn['Area Code'].astype(object)<jupyter_output><empty_output><jupyter_text>それでは次に各属性の値を、目的変数の True か False か、にわけて見てみます。<jupyter_code>for column in churn.select_dtypes(include=['object']).columns:
    if column != 'Churn?':
        display(pd.crosstab(index=churn[column], columns=churn['Churn?'], normalize='columns'))

for column in churn.select_dtypes(exclude=['object']).columns:
    print(column)
    hist = churn[[column, 'Churn?']].hist(by='Churn?', bins=30)
    plt.show()<jupyter_output><empty_output><jupyter_text>データ分析の結果から、離反する顧客について、以下のような傾向が考えられます。- 地理的にもほぼ一様に分散している- 国際通話を利用している- VoiceMailを利用していない- 通話時間で見ると長い通話時間と短い通話時間の人に分かれる- カスタマーサービスへの通話が多い (多くの問題を経験した顧客ほど離反するというのは理解できる)加えて、離反する顧客に関しては、`Day Mins` と `Day Charge` で似たような分布を示しています。しかし、話せば話すほど、通常課金されるので、驚くことではないです。もう少し深く調べてみましょう。`corr()` を利用すると相関係数を求めることができます。<jupyter_code>display(churn.corr())
pd.plotting.scatter_matrix(churn, figsize=(12, 12))
plt.show()<jupyter_output><empty_output><jupyter_text>いくつかの特徴は互いに100%の相関をもっています。このような特徴があるとき、機械学習のアルゴリズムによっては全くうまくいかないことがあり、そうでなくても結果が偏ったりしてしまうことがあります。これらの相関の強いペアは削除しましょう。Day Mins に対する Day Charge、Night Mins に対する Night Charge、Intl Mins に対する Intl Charge を削除します。<jupyter_code>churn = churn.drop(['Day Charge', 'Eve Charge', 'Night Charge', 'Intl Charge'], axis=1)<jupyter_output><empty_output><jupyter_text>ここまででデータセットの前処理は完了です。これから利用する機械学習のアルゴリズムを決めましょう。前述したように、数値の大小 (中間のような数値ではなく)で離反を予測するような変数を用意すると良さそうです。線形回帰のようなアルゴリズムでこれを行う場合は、複数の項（もしくはそれらをまとめた項）を属性として用意する必要があります。そのかわりに、これを勾配ブースティング木 (Gradient Boosted Tree)を利用しましょう。Amazon SageMaker は、マネージドで、分散学習が設定済みで、リアルタイム推論のためのホスティングも可能な XGBoost コンテナを用意しています。XGBoost は、特徴感の非線形な関係を考慮した勾配ブースティング木を利用しており、特徴間の複雑な関連性を扱うことができます。Amazon SageMaker の XGBoostは、csv または LibSVM 形式のデータを学習することができます。ここでは csv を利用します。csv は以下のようなデータである必要があります。- 1列目が予測対象のデータ- ヘッダ行はなしまずはじめに、カテゴリ変数を数値データに変換する必要があります。`get_dummies()` を利用すると数値データへの変換が可能です。そして、`Churn?_True`のデータを最初の列にもってきて、`Churn?_False.`, `Churn?_True.`のデータを削除した残りのデータをconcatenate (連結) します。<jupyter_code>model_data = pd.get_dummies(churn)
model_data = pd.concat([model_data['Churn?_True.'], model_data.drop(['Churn?_False.', 'Churn?_True.'], axis=1)], axis=1)<jupyter_output><empty_output><jupyter_text>ここで学習用、バリデーション用、テスト用データにわけましょう。これによって overfitting (学習用データには精度が良いが、実際に利用すると制度が悪い、といった状況) を回避しやすくなり、未知のテストデータに対する精度を確認することができます。<jupyter_code>train_data, validation_data, test_data = np.split(model_data.sample(frac=1, random_state=1729), [int(0.7 * len(model_data)), int(0.9 * len(model_data))])
train_data.to_csv('train.csv', header=False, index=False)
validation_data.to_csv('validation.csv', header=False, index=False)<jupyter_output><empty_output><jupyter_text>学習には学習用データとバリデーション用データのみが必要です。上で csv に出力したデータをS3にアップロードして学習に利用できるようにします。<jupyter_code>sagemaker_session = sagemaker.Session()
input_train = sagemaker_session.upload_data(path='train.csv', key_prefix='sagemaker/DEMO-xgboost-churn')
input_validation = sagemaker_session.upload_data(path='validation.csv', key_prefix='sagemaker/DEMO-xgboost-churn')<jupyter_output><empty_output><jupyter_text>`input_train` と `input_validation` にはアップロードしたファイルのS3パスが保存されています。これらは csv ファイルですが、Amazon SageMaker が用意している XGBoost のコンテナは、ファイルをデフォルトで libsvm 形式と認識してしまうため、このままだとエラーが発生します。`TrainingInput`という関数を利用して、`content_type='text/csv'`を明示的に指定することで、csv 形式と認識させることができます。<jupyter_code># from sagemaker.session import s3_input
from sagemaker.inputs import TrainingInput

content_type='text/csv'
s3_input_train = TrainingInput(input_train, content_type=content_type)
s3_input_validation = TrainingInput(input_validation, content_type=content_type)<jupyter_output><empty_output><jupyter_text>--- 4.学習それでは学習を始めましょう。まず、XGBoost のコンテナの場所を取得します。コンテナ自体は SageMaker 側で用意されているので、場所を指定すれば利用可能です。<jupyter_code>container = sagemaker.image_uris.retrieve("xgboost", boto3.Session().region_name, "1.2-1")<jupyter_output><empty_output><jupyter_text>学習のためにハイパーパラメータを指定したり、学習のインスタンスの数やタイプを指定することができます。XGBoost における主要なハイパーパラメータは以下のとおりです。- `max_depth` アルゴリズムが構築する木の深さをコントロールします。深い木はより学習データに適合しますが、計算も多く必要で、overfiting になる可能性があります。たくさんの浅い木を利用するか、少数の深い木を利用するか、モデルの性能という面ではトレードオフがあります。- `subsample` 学習データのサンプリングをコントロールします。これは overfitting のリスクを減らしますが、小さすぎるとモデルのデータが不足してしまいます。-  `num_round` ブースティングを行う回数をコントロールします。以前のイテレーションで学習したときの残差を、以降のモデルにどこまで利用するかどうかを決定します。多くの回数を指定すると学習データに適合しますが、計算も多く必要で、overfiting になる可能性があります。- `eta` 各ブースティングの影響の大きさを表します。大きい値は保守的なブースティングを行います。- `gamma` ツリーの成長の度合いをコントロールします。大きい値はより保守的なモデルを生成します。XGBoostのhyperparameterに関する詳細は [github](https://github.com/dmlc/xgboost/blob/master/doc/parameter.rst) もチェックしてください。<jupyter_code>sess = sagemaker.Session()

hyperparameters = {"max_depth":"5",
                        "eta":"0.2",
                        "gamma":"4",
                        "min_child_weight":"6",
                        "subsample":"0.8",
                        "objective":"binary:logistic",
                        "num_round":"100"}

xgb = sagemaker.estimator.Estimator(container,
                                    role, 
                                    hyperparameters=hyperparameters,
                                    instance_count=1, 
                                    instance_type='ml.m4.xlarge',
                                    sagemaker_session=sess)

xgb.fit({'train': s3_input_train, 'validation': s3_input_validation})<jupyter_output><empty_output><jupyter_text>--- 5.ホスティング学習が終われば、`deploy()`を実行することで、エンドポイントを作成してモデルをデプロイできます。<jupyter_code>xgb_predictor = xgb.deploy(initial_instance_count = 1, instance_type = 'ml.m4.xlarge')<jupyter_output><empty_output><jupyter_text>5-1.評価現在、エンドポイントをホストしている状態で、これを利用して簡単に予測を行うことができます。予測は http の POST の request を送るだけです。ここではデータを `numpy` の `array` の形式で送って、予測を得られるようにしたいと思います。しかし、endpoint は `numpy` の `array` を受け取ることはできません。このために、`csv_serializer` を利用して、csv 形式に変換して送ることができます。<jupyter_code>xgb_predictor.serializer = sagemaker.serializers.CSVSerializer()<jupyter_output><empty_output><jupyter_text>作成済みのテストデータを受け取ると、これをデフォルト500行ずつのデータにわけて、エンドポイントに送信する `predict` という関数を用意します。あとは `predict` を実行して予測結果を受け取ります。<jupyter_code>def predict(data, rows=500):
    split_array = np.array_split(data, int(data.shape[0] / float(rows) + 1))
    predictions = ''
    for array in split_array:
        predictions = ','.join([predictions, xgb_predictor.predict(array).decode('utf-8')])

    return np.fromstring(predictions[1:], sep=',')

dtest = test_data.values
predictions = []
predictions.append(predict(dtest[:, 1:]))
predictions = np.array(predictions).squeeze()<jupyter_output><empty_output><jupyter_text>機械学習の性能を比較評価する方法はいくつかありますが、単純に、予測値と実際の値を比較しましょう。今回は、顧客が離反する `1` と離反しない `0` を予測しますので、この混同行列を作成します。<jupyter_code>pd.crosstab(index=test_data.iloc[:, 0], columns=np.round(predictions), rownames=['actual'], colnames=['predictions'])<jupyter_output><empty_output><jupyter_text>_注意点, アルゴリズムにはランダムな要素があるので結果は必ずしも一致しません._48人の離反者がいて、それらの39名 (true positives) を正しく予測できました。そして、4名の顧客は離反すると予測しましたが、離反していません (false positives)。9名の顧客は離反しないと予測したにもかかわらず離反してしまいました (false negatives)。重要な点として、離反するかどうかを `np.round()` という関数で、しきい値0.5で判断しています。`xgboost` が出力する値は0から1までの連続値で、それらを離反する `1` と 離反しない `0` に分類します。しかし、その連続値 (離反する確率) が示すよりも、顧客の離反というのは損害の大きい問題です。つまり離反する確率が低い顧客も、しきい値を0.5から下げて、離反するとみなす必要があるかもしれません。もちろんこては、false positives （離反すると予測したけど離反しなかった）を増やすと思いますが、 true positives (離反すると予測して離反した) を増やし、false negatives (離反しないと予測して離反した）を減らせます。直感的な理解のため、予測結果の連続値をみてみましょう。<jupyter_code>plt.hist(predictions)
plt.show()<jupyter_output><empty_output><jupyter_text>連続値は0から1まで歪んでいますが、0.1から0.9までの間で、しきい値を調整するにはちょうど良さそうです。<jupyter_code>pd.crosstab(index=test_data.iloc[:, 0], columns=np.where(predictions > 0.3, 1, 0))<jupyter_output><empty_output><jupyter_text>例えば、しきい値を0.5から0.3まで減らしてみたとき、true positives は 1 つ、false positives は 3 つ増え、false negatives は 1 つ減りました。全体からみると小さな値ですが、全体の6-10%の顧客が、しきい値の変更で、予測結果が変わりました。ここで5名にインセンティブを与えることによって、インセンティブのコストが掛かりますが、3名の顧客を引き止めることができるかもしれません。つまり、最適な閾値を決めることは、実世界の問題を機械学習で解く上で重要なのです。これについてもう少し広く議論し、仮説的なソリューションを考えたいと思います。 5-2.推論エラーのコスト2値分類の問題においては、しきい値に注意しなければならないという、似たような状況に直面することが多いです。それ自体は問題ではありません。もし、出力の連続値が2クラスで完全に別れていれば、MLを使うことなく単純なルールで解くことができると考えられます。重要なこととして、MLモデルを正版環境に導入する際、モデルが false positives と false negatives に誤って入れたときのコストがあげられます。しきい値の選択は4つの指標に影響を与えます。4つの指標に対して、ビジネス上の相対的なコストを考える必要があるでしょう。 コストの割当携帯電話会社の離反の問題において、コストとはなんでしょうか？コストはビジネスでとるべきアクションに結びついています。いくつかの仮定をおいてみましょう。まず、true negatives のコストとして \$0 を割り当てます。満足しているお客様を正しく認識できていれば何も実施しません。false negatives が一番問題で、なぜなら、離反していく顧客を正しく予測できないからです。顧客を失えば、再獲得するまでに多くのコストを払う必要もあり、例えば逸失利益、広告コスト、管理コスト、販売管理コスト、電話の購入補助金などがあります。インターネットを簡単に検索してみると、そのようなコストは数百ドルとも言われ、ここでは `$500` としましょう。これが false negatives に対するコストです。最後に、離反していくと予測された顧客に `$100` のインセンティブを与えることを考えましょう。携帯電話会社がそういったインセンティブを提供するなら、2回くらいは離反の前に考え直すかもしれません。これは true positive と　false negative のコストになります。false positives の場合 (顧客は満足していて、モデルが誤って離反しそうと予測した場合)、 `$100` のインセンティブは捨てることになります。その `$100` を効率よく消費してしまうかもしれませんが、優良顧客へのロイヤリティを増やすという意味では悪くないかもしれません。  5-3.最適な閾値を探すfalse negatives が false positives よりもコストが高いことは説明しました。そこで、顧客の数ではなく、コストを最小化するように、しきい値を最適化することを考えましょう。コストの関数は以下のようなものになります。```txt$500 * FN(C) + $0 * TN(C) + $100 * FP(C) + $100 * TP(C)```FN(C) は false negative の割合で、しきい値Cの関数です。同様にTN, FP, TP も用意します。この関数の値が最小となるようなしきい値Cを探します。最も単純な方法は、候補となる閾値で何度もシミュレーションをすることです。以下では100個の値に対してループで計算を行います。<jupyter_code>cutoffs = np.arange(0.01, 1, 0.01)
costs = []

for c in cutoffs:
    _predictions = pd.Categorical(np.where(predictions > c, 1, 0), categories=[0, 1])
    matrix_a = np.array([[0, 100], [500, 100]])
    matrix_b = pd.crosstab(index=test_data.iloc[:, 0], columns=_predictions, dropna=False)
    costs.append(np.sum(np.sum(matrix_a * matrix_b)))

costs = np.array(costs)
plt.plot(cutoffs, costs)
plt.show()
print('Cost is minimized near a cutoff of:', cutoffs[np.argmin(costs)], 'for a cost of:', np.min(costs))<jupyter_output><empty_output><jupyter_text>6.エンドポイントの削除エンドポイントは起動したままだとコストがかかります。不要な場合は削除します。<jupyter_code>xgb_predictor.delete_endpoint()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11528185947101256, 'fraction_numerical': 0.015428800427464601, 'mean_word_length': 9.604107648725213, 'pattern_counts': {'":': 7, '<': 99, '<?xml version=': 0, '>': 102, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '243051', 'license': 'Apache-2.0', 'n_tokens_mistral': 10117, 'n_tokens_neox': 9405, 'n_words': 670, 'path': 'xgboost_customer_churn/xgboost_customer_churn.ipynb', 'repo_name': 'kazuhitogo/amazon-sagemaker-examples-jp', 'size': '14982'}