Pythonで形態素解析器Sudachiを使う (SudachiPy)

仕事でSudachiをお試しする機会がありましたので、基本的な使い方を備忘録として整理しておきます。

Sudachi

Sudachiはワークスアプリケーションズさんで開発されたOSSで、日本語に特化した形態素解析器です。分かち書き・品詞付け・正規化を行います。

github.com

大きな特徴は3つです。特に3番目は実務上もどかしい問題です。例えば"焼き肉"と"焼肉"、"醤油"と"しょう油"といった表現の違いを、Sudachiで正規化することで同じ単語であると認識できるようになります。

  • UniDicとNeologdがベースとした豊富な語彙
  • 分割単位を3段階 (辞書単位から固有表現単位まで) から選択可能
  • 送り仮名や字種の違いなども正規化

SudachiPy

SudachiをPythonから使うためのライブラリとしてSudachiPyも開発されています (同じくワークスアプリケーションズさん主体のOSS) 。今回もこれを使います。

github.com

インストール

SudachiPyのインストール -> 辞書のダウンロードと配置 の順番で行います (なおPythonは3.6です) 。

最初にSudachiPyをインストールします。このときインストールしたパスを確認しておきます。

$ pip install -e git+git://github.com/WorksApplications/SudachiPy@develop#egg=SudachiPy
$ pip list | grep sudachipy
SudachiPy                          0.1.0     /Users/ohke/src/sudachipy

次にSudachiのレポジトリから辞書をダウンロードして、SudachiPyのインストールディレクトリに配置します。

  • 上で確認しましたインストールディレクトリ直下のresources配下に"system.dic"に配置します
    • /Users/ohke/src/sudachipy/resources/system.dic
  • 辞書にはcoreとfullの2つがあり、fullの方には雑多な固有表現まで含まれます
    • 今回はcoreを使いますが、fullも同じ手順でOKです
$ wget https://github.com/WorksApplications/Sudachi/releases/download/v0.1.1/sudachi-0.1.1-dictionary-core.zip
$ unzip sudachi-0.1.1-dictionary-core.zip
$ mv system_core.dic /Users/ohke/src/sudachipy/resources/system.dic
$ ls /Users/ohke/src/sudachipy/resources/
char.def        rewrite.def     sudachi.json    system.dic      unk.def

ここまでで準備OKです。

コマンドラインから使う

まずコマンドラインから使ってみます。sudachipyで起動します。

  • -mで分割モード (A, B, C) を指定できます
    • Aは"スカイ"と"ツリー"に名詞として分割されましたが、Cは"スカイツリー"という1つの固有名詞として分割されていることがわかります
$ sudachipy -m A -a
友人・我孫子と東京スカイツリーへ走って向かった。
友人    名詞,普通名詞,一般,*,*,*        友人    友人    ユウジン        0
・      補助記号,一般,*,*,*,*   ・      ・              0
我孫子  名詞,固有名詞,地名,一般,*,*     我孫子  我孫子  アビコ  0
と      助詞,格助詞,*,*,*,*     と      と      ト      0
スカイ  名詞,普通名詞,一般,*,*,*        スカイ  スカイ          0
ツリー  名詞,普通名詞,一般,*,*,*        ツリー  ツリー          0
へ      助詞,格助詞,*,*,*,*     へ      へ      ヘ      0
走っ    動詞,一般,*,*,五段-ラ行,連用形-促音便   走る    走る    ハシッ  0
て      助詞,接続助詞,*,*,*,*   て      て      テ      0
向かっ  動詞,一般,*,*,五段-ワア行,連用形-促音便 向かう  向かう  ムカッ  0
た      助動詞,*,*,*,助動詞-タ,終止形-一般      た      た      タ      0
。      補助記号,句点,*,*,*,*   。      。              0
EOS

$ sudachipy -m C -a
友人・我孫子とスカイツリーへ走って向かった。
友人    名詞,普通名詞,一般,*,*,*        友人    友人    ユウジン        0
・      補助記号,一般,*,*,*,*   ・      ・              0
我孫子  名詞,固有名詞,地名,一般,*,*     我孫子  我孫子  アビコ  0
と      助詞,格助詞,*,*,*,*     と      と      ト      0
スカイツリー    名詞,固有名詞,一般,*,*,*        スカイツリー    スカイツリー            0
へ      助詞,格助詞,*,*,*,*     へ      へ      ヘ      0
走っ    動詞,一般,*,*,五段-ラ行,連用形-促音便   走る    走る    ハシッ  0
て      助詞,接続助詞,*,*,*,*   て      て      テ      0
向かっ  動詞,一般,*,*,五段-ワア行,連用形-促音便 向かう  向かう  ムカッ  0
た      助動詞,*,*,*,助動詞-タ,終止形-一般      た      た      タ      0
。      補助記号,句点,*,*,*,*   。      。              0
EOS

Pythonコードから使う

PythonコードからSudachiPyを使って形態素解析します。

最初にconfig.SETTINGFILEの設定ファイルをJSONでロードします。

  • このファイルがSudachiPyのデフォルトの設定になります
    • このJSONを変更することで、デフォルトの設定からカスタマイズできます
import json
from sudachipy import config
from sudachipy import dictionary

print(config.SETTINGFILE)
# /Users/ohke/src/sudachipy/sudachipy/../resources/sudachi.json

with open(config.SETTINGFILE, 'r', encoding='utf-8') as f:
    settings = json.load(f)
    
print(settings)
# {'systemDict': 'system.dic',
#  'characterDefinitionFile': 'char.def',
#  'inputTextPlugin': [{'class': 'com.worksap.nlp.sudachi.DefaultInputTextPlugin'}],
#  'oovProviderPlugin': [{'class': 'com.worksap.nlp.sudachi.MeCabOovProviderPlugin',
#    'charDef': 'char.def',
#    'unkDef': 'unk.def'},
#   {'class': 'com.worksap.nlp.sudachi.SimpleOovProviderPlugin',
#    'oovPOS': ['補助記号', '一般', '*', '*', '*', '*'],
#    'leftId': 5968,
#    'rightId': 5968,
#    'cost': 3857}],
#  'pathRewritePlugin': [{'class': 'com.worksap.nlp.sudachi.JoinNumericPlugin',
#    'joinKanjiNumeric': True},
#   {'class': 'com.worksap.nlp.sudachi.JoinKatakanaOovPlugin',
#    'oovPOS': ['名詞', '普通名詞', '一般', '*', '*', '*']}]}

そして上の設定を反映したトークナイザを生成し、tokenizeメソッドを呼び出すことで、形態素解析できます。

  • コマンドラインと同様に、tokenizer.Tokenizer.SplitModeの1つを引数に渡すことで分割モードを指定できます (ここではCを使ってます)
  • 結果はMorphemeのリストとなっており、表層形 (surface()) 、品詞 (part_of_speech()) 、読み (reading_form()) 、正規化した表現 (normalized_form()) を取得できます
from sudachipy import tokenizer

tokenizer_obj = dictionary.Dictionary(settings).create()
print(type(tokenizer_obj))
# <class 'sudachipy.tokenizer.Tokenizer'>

text = '友人・我孫子とスカイツリーでスパゲティを食った。'

tokens = tokenizer_obj.tokenize(tokenizer.Tokenizer.SplitMode.C, text)
print(type(tokens))
# <class 'sudachipy.morphemelist.MorphemeList'>

for t in tokens:
    print(t.surface(), t.part_of_speech(), t.reading_form(), t.normalized_form())

形態素解析後の結果を見ると以下のようになります。品詞はリストになっていること、"スパゲティ" -> "スパゲッティ" と正規化されていることに注目してください。

友人 ['名詞', '普通名詞', '一般', '*', '*', '*'] ユウジン 友人
・ ['補助記号', '一般', '*', '*', '*', '*']  ・
我孫子 ['名詞', '固有名詞', '地名', '一般', '*', '*'] アビコ 我孫子
と ['助詞', '格助詞', '*', '*', '*', '*'] ト と
スカイツリー ['名詞', '固有名詞', '一般', '*', '*', '*']  スカイツリー
で ['助詞', '格助詞', '*', '*', '*', '*'] デ で
スパゲティ ['名詞', '普通名詞', '一般', '*', '*', '*']  スパゲッティ
を ['助詞', '格助詞', '*', '*', '*', '*'] ヲ を
食っ ['動詞', '一般', '*', '*', '五段-ワア行', '連用形-促音便'] クッ 食う
た ['助動詞', '*', '*', '*', '助動詞-タ', '終止形-一般'] タ た
。 ['補助記号', '句点', '*', '*', '*', '*']  。

まとめ

SudachiをPythonから使って形態素解析を行いました。

Python: lifetimesで顧客のトランザクション数を予測する

前回紹介したBG/NBDなど、顧客の未来のライフタイムバリュー (ここではトランザクション数) を予測するモデルを、簡単に利用できるPythonパッケージ lifetimes を紹介します。

lifetimes

lifetimesは、顧客の未来のトランザクションを計算するためためのパッケージです。
BG/NBD以外にも、ベースとなっているPareto/NBDや、トランザクション数ではなく金額を直接的に予測するGamma-Gammaモデル (提案論文) なども提供されています。

github.com

lifetimesを使った予測

$ pip install lifetimes

データセット

今回は、UCI Machine Learning Repositoryで提供されているECサイトの取引データを使って予測を行っていきます。

2010/12/1〜2011/12/9に発生した取引 (=トランザクション) 約54万件 (キャンセル含む) が収録されており、顧客ID・取引日時・金額などの情報を含んでいます。

UCI Machine Learning Repository: Online Retail Data Set

最初に、このデータをPandasのDataFrameにロードします。

import pandas as pd
import urllib.request

url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00352/Online%20Retail.xlsx'
path = './OnlineRetail.xlsx'

# xlsxのダウンロード
urllib.request.urlretrieve(url, path)

# DataFrameにロード
df = pd.read_excel(path)
df.head()

で、ちょこっと加工して、顧客ID、購入日時、合計金額のレコードに集約します。

  • キャンセルのトランザクション (InvoiceNoが'C'で始まる) はここで除去してしまいます
    • キャンセルされた取引がどの購入取引と紐付いているのかがわからないため
  • 簡単化のため、購入日時 (InvoiceDate) が同じ場合は同じトランザクションとしてます
# キャンセル分は、削除
df = df.loc[df['SubTotalPrice'] > 0.0]

# 合計金額を計算
df['SubTotalPrice'] = df['Quantity'] * df['UnitPrice']

# 顧客ID, 購入日時, 金額 のレコードを生成
transaction_df = df.groupby(['CustomerID', 'InvoiceDate'])['SubTotalPrice'].sum().reset_index()

将来のトランザクション数の予測 (BG/NBDモデル)

トランザクション数をBG/NBDモデルを使って予測していきます。

復習となりますが、BG/NBDでは、顧客ごとに3つのデータを使います。詳しくは前回の投稿をご覧ください。

  • T: 最初にトランザクションしてから現在までの期間 (日数)
  • frequency: 期間T中のトランザクション数
  • reacency: 期間T中で最後にトランザクションしてから現在までの期間 (日数)

トランザクションデータから上の3つの項目を持つデータを作り、BG/NBDモデルを学習させます。

  • PandasのDataFrameに対応したインタフェースとなってます
  • トランザクションごとにレコードとなっているデータでも、lifetimes.utils.summary_data_from_transaction_dataメソッドを使うと、上3つのカラムを持つ顧客ID単位のテーブルに集約してくれます (自前でGroupbyする必要がありません)
    • observation_period_endは2011/11/30として、ちょうど1年分となるように調整してます
# transaction_dfから、frequency, recency, Tを抽出
from lifetimes.utils import summary_data_from_transaction_data

summary_df = summary_data_from_transaction_data(
    transaction_df, 'CustomerID', 'InvoiceDate', 
    observation_period_end='2011-11-30'
)
summary_df.head()

モデルを生成して学習します。

  • BG/NBDモデルはlifetimes.filters.BetaGeoFilter
    • 他のモデルもlifetimes.filters.BaseFitterを継承したモデルとして実装されてます
    • penalizer_coefは正則化の係数で、ここでは0に設定することで正則化無しのモデルとしてます
# BG/NBDモデルを学習する
from lifetimes import BetaGeoFitter

bgf_model = BetaGeoFitter(penalizer_coef=0.0)

bgf_model.fit(summary_df['frequency'], summary_df['recency'], summary_df['T'])
# <lifetimes.BetaGeoFitter: fitted with 4297 subjects, a: 0.00, alpha: 70.93, b: 3.13, r: 0.83>

学習の結果、各パラメータは以下の値で推定されています。

パラメータ
r 0.83
alpha 70.96
a 0.00
b 2.00

モデルの可視化をしていきます。

1日後の購入数はplot_frequency_recency_matrixメソッドで可視化できます。

import matplotlib
from lifetimes.plotting import plot_frequency_recency_matrix

plot_frequency_recency_matrix(bgf_model)

f:id:ohke:20190302104443p:plain

生存確率はplot_probability_alive_matrixです。recencyが小さいほど・frequencyが大きいほど、確率が高く計算されることがわかります。

from lifetimes.plotting import plot_probability_alive_matrix

plot_probability_alive_matrix(bgf_model)

f:id:ohke:20190302104456p:plain

このモデルを使って予測していきます。BetaGeoFilterのconditional_expected_number_of_purchases_up_to_timeメソッドで、未来の期間t (日) でのトランザクション数を予測できます。

  • ここではt=7 (12/1〜12/7の1週間) を指定してます
# 未来の7日間 (2011-12-01〜2011-12-07) のトランザクション数を顧客ごとに予測
t = 7
summary_df['predicted_purchases'] = bgf_model.conditional_expected_number_of_purchases_up_to_time(
    t, summary_df['frequency'], summary_df['recency'], summary_df['T']
)

# 予測トランザクション数が多いトップ5
summary_df.sort_values(by='predicted_purchases', ascending=False).head()

実際に2011/12/1〜2011/12/7に発生したトランザクション数のTOP5と比較すると、5人中2人が含まれており、傾向レベルでは予測できていそうです。

# 実際のトランザクション数TOP5を表示
transaction_df.loc[
    (transaction_df['InvoiceDate'] >= '2011-12-01') &
    (transaction_df['InvoiceDate'] < '2011-12-07')
].groupby('CustomerID').count().sort_values(by='InvoiceDate', ascending=False)['InvoiceDate'].head()
# CustomerID
# 15856.0    6
# 14911.0    5
# 12748.0    5
# 12569.0    4
# 13089.0    3

まとめ

lifetimesパッケージを使うことで、顧客のトランザクション数を予測できることを確認しました。

ライフタイムバリューを予測する (BG/NBDモデル)

今回は、顧客のライフタイムバリューを予測する方法の1つとして、BG/NBDモデルを紹介します。

モチベーション

顧客のライフタイムバリューを予測できると、その顧客に対してどれだけ投資して良いかがわかります。ここで言う投資は、クーポンやポイントなどのインセンティブ、DMや電話などの営業努力のことです。

また、離脱 (churn) の予測にも転用できます。過去に売上に大きく貢献してきた人が、未来のライフタイムバリューはそれと比べると極めて小さい、または、ゼロとなることが予測されるケースです。
新規顧客の獲得は既存顧客の維持よりも遥かに大変なことはよく知られてます。また、過去には大きく売上に貢献していたので、引き止められれば優良顧客まで戻る見込みは高いはずです。

BG/NBDモデル

このライフタイムバリューを予測するモデルの1つに、BG/NBD (Beta-Geometric / Negative Binomial Distribution) が2005年に提案されています (PDF) 。

repository.upenn.edu

先行していたPareto/NBD (提案論文) の簡易版として提案されたモデルです。

  • 顧客ごとの過去のアクション (ここではトランザクションという) を使う
    • 期間Tにおけるトランザクション数  x
      • Tは最初のトランザクションの時刻  t_0 を起点
    • 最後にトランザクションした時間  t_x
  • 予測したいことは3つ
    • 長さtの時間のトランザクション数の (全体の) 期待値  E[X(t) ]
    • 長さtの時間でトランザクション数xとなる (全体の) 確率  P[X(t)=x ]
    • (T, T+t] の時間における、ある顧客 (X=x, t, T) のトランザクション数の期待値  E(\Upsilon(t) | X=x, t, T)

ある顧客Aをイメージするとこんな感じです。期間T中にx回トランザクションが発生し、最後のトランザクションは時間  t_x でした。では、Tの後の時間tの間に何回トランザクションが発生するか、を予測したいのです。

前提

BG/NBDモデルは5つの前提を置いています。

  • 顧客は アクティブ (active) と 離脱 (inactive) の2つの状態を持つ
    • アクティブな顧客のみがトランザクションを行う

前提1. アクティブな顧客はポアソン分布に従ってトランザクションする


f(t_j | t_{j-1}, \lambda) = \lambda e^{t_j - t_{j-1}}, t_{j} > t_{j-1} \geq 0

前提2. λはガンマ分布に従う


f(\lambda | r, \alpha) = \frac{\alpha^{r} \lambda^{r-1} e^{-\lambda\alpha}}{\Gamma(r)}, \lambda>0

前提3. トランザクション後に確率pで離脱する

Pareto/NBDと異なる前提を置いています。Pareto/NBDでは、離脱は任意のタイミングで起きる前提に立っています。

P(inactive immediately after j-th transaction) = 
p(1-p)^{j-1}, j=1, 2, ...

前提4. pはベータ分布に従う

B(a, b)はベータ関数です。


f(p|a, b) = \frac{p^{a-1} (1-p)^{b-1}}{B(a, b)}, 0 \leq p \leq 1


B(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}

前提5. λとpは顧客間で独立

そのため前提1.〜4.で挙げたパラメータがわかれば、他の顧客を考慮すること無く、λとpを計算できます。

パラメータの推定

上の前提をまとめると、求める必要があるパラメータは  r, \alpha, a, b の4つとなります。尤度関数を定義してパラメータを求めます。

尤度関数Lは以下の式となります。 \delta_{x>0} は x>0のときに1、x=0で0となります。したがって、まだ1回もトランザクションしていない場合は、第1項だけで計算されます。


L(r \alpha, a, b | X=x, t_x, T) = \frac{B(a, b + x)}{B(a, b)} \frac{\Gamma(r + x)a^r}{\Gamma(r)(\alpha + T)^{r + x}} + \delta_{x>0} \frac{B(a+1, b+x-1)}{B(a, b)} \frac{\Gamma(r + x)\alpha^r}{\Gamma(r)(a + t_x)^{r+x}}

最大化する対数尤度関数LLは以下の式です。Nは顧客数です。


LL(r, \alpha, a, b) = \sum_{i=1}^{N} ln[L(r, \alpha, a, b | X_i = x_i, t_{x_i}, T_i) ]

予測

上で求まったパラメータを使って、3つの値を予測します。

トランザクション数がxとなる確率 (全体)


P(X(t)=x | r, \alpha, a, b) = \frac{B(a, b+x)}{B(a, b)} \frac{\Gamma(r+x)}{\Gamma(r)x!} (\frac{\alpha}{\alpha+t})^r (\frac{1}{\alpha+t})^x + \delta_{x>0} \frac{B(a+1, b+x-1)}{B(a, b)} (1-(\frac{\alpha}{\alpha+t})^r (\sum_{j=0}^{x-1} \frac{\Gamma(r+j)}{\Gamma(r)j!}(\frac{t}{\alpha+t})^j))

トランザクション数の期待値 (全体)


E(X(t) | r, \alpha, a, b) = \frac{a+b-1}{a-1}(1-(\frac{\alpha}{\alpha+t})^r {}_2F_1(r, b; a+b-1; \frac{t}{\alpha+t}))

 {}_2F_1() はガウス超幾何関数 (Gaussian hypergeometric function) です。このあたりは勉強不足でよくわからんですが、以下で計算できるようです。


_2F_1(r, b; a+b-1; \frac{t}{\alpha+t}) = \sum_{n=0}^{\infty} \frac{(r)_n (b)_n}{(a+b-1)_n} \frac{(\frac{t}{\alpha+t})^n}{n!}

 (x)_n はポッホハマー記号で  (x)_0 := 1 ,  (x)_n := \prod_{k=0}^{n-1}(x+k) で定義されます。

ある顧客のトランザクション数の期待値


E(\Upsilon(t) | X=x, t_x, T, r, \alpha, a, b) = \frac{a+b+x-1}{a-1} \frac{1-(\frac{\alpha + T}{\alpha + T + t})^{r+x}{}_2F_1(r + x, b + x; a+b+x-1; \frac{t}{\alpha + T + t})}{1 + \delta_{x>0}\frac{a}{b+x-1}(\frac{\alpha + T}{\alpha + t_x})^{r+x}}

複雑な式に見えますが、パラメータさえ定まってしまえば、簡単に計算できます。 以下のパラメータ (原文Figure 1から抜粋) を使ってPythonで計算してみます。

パラメータ
r 0.243
α 4.414
a 0.793
b 2.426
import math
from scipy.special import hyp2f1  # 超幾何級数2F1の実装

# パラメータ
r, alpha, a, b = 0.243, 4.414, 0.793, 2.426

# トランザクション数の期待値
def expect_transaction_count(T, x, t_x, t):
    term1 = (a + b + x - 1)/(a - 1)
    
    numerator = 1 - ((alpha + T) / (alpha + T + t))**(r + x) * hyp2f1(r + x, b + x, a+b+x-1, t / (alpha + T + t))
    denominator = 1 + (a / (b + x - 1))*((alpha + T) / (alpha + t_x))**(r + x)
    term2 = numerator / denominator
    
    return term1 * term2

3パターンの顧客でトランザクション数の期待値を比較してみます。いずれもT=t=30、つまり過去の30期間 (例えば日数や週数) のxとt_xを使って、未来の30期間のトランザクション数を予測します。

  • 顧客1と顧客2を比較すると、過去のトランザクション数xが多いほうが、未来のトランザクション数も増えると予想されることがわかります
  • 顧客1と顧客3を比較すると、最後のトランザクションt_xが遠いほうが、未来のトランザクション数も減ると予想されることもわかります
# 顧客1
T, x, t_x, t = 30, 4, 25, 30
expect_transaction_count(T, x, t_x, t)
# 2.3526752183407695

# 顧客2
T, x, t_x, t = 30, 8, 25, 30
expect_transaction_count(T, x, t_x, t)
# 4.388865813295875

# 顧客3
T, x, t_x, t = 30, 4, 15, 30
expect_transaction_count(T, x, t_x, t)
# 1.1367813390968273

所感

BG/NBDモデルについての、私の所感をまとめます。

  • 期間T内のトランザクション数と最終トランザクション日時しか見ていないため、周期性のある顧客については表現できていない
    • 周期性のある顧客は、tの設定によっては、予想されるトランザクション数と実際のトランザクション数が大きく異なることが考えられる
  •  t_1  t_{x-1} のトランザクション発生日時は捨ててしまって良いのだろうか?
    • 例えば過去30日間で見たとき、同じx=10だとしても、1〜9日目と30日目にトランザクションした人と、平均的にトランザクションした人は同じと扱って良いかは疑問が残る (後者の方が継続して利用しそう)
  • あくまでトランザクション数であるため、ビジネスサイドが求めるライフタイムバリューと揃えるためには何らかの計算が必要
    • その顧客の過去の平均金額と掛け算する、など

まとめ

顧客のライフタイムバリューを表現するBG/NBDモデルを紹介しました。