らくがき入門

環境と研究テーマが大幅に変わりました。だいたい何かに入門しています。

ユーザーベース協調フィルタリングを実装してみた

amazonNetflixで使われているレコメンドシステムに突然興味を持ったので、推薦システムの勉強をしたのでpythonで実装してみました。

勉強のアウトプットなので、間違っていたらご指摘いただけると幸いです。

推薦システムとは

amazonでモノを買ったりすると、「この商品をチェックした人はこんな商品もチェックしています」みたいな推薦をしてくれますよね。 まさにあれが推薦システムで、特定のユーザーに対してどのモノを推薦すべきかを決定するシステムのことを推薦システムといいます。 その中でも、すべてのユーザーの嗜好に応じて異なる推薦のリストを提示するシステムが個人化推薦といいます。

種類としては、

  • 協調型推薦
  • 内容ベース型推薦
  • 知識型ベース型
  • 複数のアプローチを組み合わせたハイブリッド型

があります。

今回は、協調型推薦の中のユーザーに基づいた推薦が対象です。

協調型推薦

協調型推薦は、 過去に同じ興味を共有したユーザーは将来的にも同じような興味をもつだろう 仮定しています。

なので、ユーザーAとユーザーBの購入履歴が似ていて、ユーザーBがまだ知らないモノをユーザーAが最近購入したとすると、 これをユーザーBに提案することは合理的だといえます。

ユーザーがお互いに示し合わせたわけではなく暗黙的に協調して、膨大なモノの集合から最も有望なモノをフィルタリングするため、 この技術を協調フィルタリングと呼んでいます。

協調型推薦を考える上で、以下の問いに答える必要がでてきます。

  • 類似したユーザーというけれど類似したユーザーを膨大なユーザーの中から探し出すのか?
  • そもそも類似ってなに?
  • 類似ってどうやって測る?
  • 新規のユーザーでまだシステムの利用履歴がないけどどう対処する?
  • ユーザーに関して利用できるデータが少ない場合はどうする?
  • 類似ユーザーを探す以外に特定のユーザーがあるモノを好きかどうかを予測できる手法は存在する?

ユーザーベース推薦

ユーザーベース推薦のシステムのアイディアは非常にシンプルで、

  • 対象ユーザーの評価データ入力して、過去の嗜好と似ている他のユーザー(=ピアユーザー)を特定する
  • 対象ユーザーがまだ見ていないモノに対する評価をピアユーザーの評価を用いて予測する

ユーザーベース推薦は以下の2つの仮説に基づいています。

  • ユーザーが過去に似た嗜好をもっているなら、その嗜好は将来においても似ている
  • ユーザーの好みは長い間一貫している

実際この辺りの仮定は少し無理があるので、amazonとかは別の手法を用いているみたいです。この辺りはまだ知識不足なので、よくわかってません笑

このあたりを基礎知識として、ユーザーベース推薦を実装してみました。

使用するデータセット

今回利用するのはMovieLensデータセットです。 https://grouplens.org/datasets/movielens/

MovieLensデータセットは推薦システムの開発やベンチマークとしてミネソタ大学の研究グループが公開してくれています。ありがたいですね。

今回はMovieLens 100K Datasetを使用します。

ベンチマーク(参考コード)

推薦システムの実装で使われていることが多いので、他の人のコードを眺めながらできるのは非常にありがたいです。 ちなみに以下のgithubのコードをベンチマークとして実装しています。

今回の実装はMovieLens100kに対するユーザーベース協調フィルタリングになるので、 precisionで19.69%を目標とすることになります。

実装

データの読み込み

データは手元にダウンロードして、使用しています。 ua.baseがtrain dataで、ua.testがtest dataとなっているみたいです。

import numpy as np
import pandas as pd

train = pd.read_csv('../data/input/ml-100k/ua.base', names=["user_id", "item_id", "rating", "timestamp"], sep='\t')
test = pd.read_csv('../data/input/ml-100k/ua.test', names=['user_id', 'item_id', 'rating', 'timestamp'], sep='\t')

データとしては以下のような内容です。 データの全容としては、943ユーザーによる1682アイテムに対して1~5の評価をしているレビューサイトのデータという感じです。

amazonをイメージしてもらえると理解の助けになると思うのですが、全員がすべてのアイテムのレビューをしているわけではないので、非常に疎なデータセットになっています。

user_id item_id rating timestamp
0 186 302 3 891717742
1 22 377 1 878887116
2 244 51 2 880606923
3 166 346 1 886397596
4 298 474 4 884182806
  • 類似したユーザーというけれど類似したユーザーを膨大なユーザーの中から探し出すのか?

→今回はすべてのユーザーと総当たりで類以度を計算して、上位何人かを類似したユーザーと定義します。

  • そもそも類似ってなに?

→類以度という概念を導入してユーザー間の類似を表現します。

  • 類似ってどうやって測る?

→類似度にも様々な表現方法があるのですが、今回はコサイン類以度を採用します。

  • 新規のユーザーでまだシステムの利用履歴がないけどどう対処する?
  • ユーザーに関して利用できるデータが少ない場合はどうする?
  • 類似ユーザーを探す以外に特定のユーザーがあるモノを好きかどうかを予測できる手法は存在する?

→上の3つに関しては今回は範囲外としてまた別のタイミングで扱うことにします。

とりあえず今回の記事で扱う話を整理したので、それぞれに注目して行きたいと思います。

コサイン類以度

今回類以度として採用したコサイン類以度は2つのベクトルを用いて以下のように定式化できます。  \displaystyle
  sim(\boldsymbol{a}, \boldsymbol{b}) = \frac{\boldsymbol{a} \cdot \boldsymbol{b}}{|\boldsymbol{a} |\times |\boldsymbol{b}|}

コサイン類以度は、0から1の値を取りベクトル同士の成す角度の近さを表現しています。

推薦システムの文脈でコサイン類以度を考える場合には、ユーザーが異なることを気をつける必要があります。例えば、ユーザーAは甘めに評価する一方で、ユーザーBは辛めに評価する傾向にあるといったところを考慮する必要があるということです。 これは、ユーザーの評価値の平均をそのユーザーの評価から引くことで解決でき、これを調整コサイン類以度といいます。

実装としては、以下のように実装しています。下のmean_adjustment=Trueとすると調整コサイン類以度になります。

def cosine_similarity(v1, v2, mean_adjustment=False):
    if mean_adjustment:
        v1 = v1 - np.mean(v1)
        v2 = v2 - np.mean(v2)
    return np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2))

下準備

ここではindexがuser_id、columnsがitem_idになるような評価値行列を作成しています。 ユーザーが今回評価をつけていないところは0としました。調べてみると、Rの推薦システムのライブラリrecommenderlabも欠損は0としているみたいなので、とりあえずこれで進めていきます。

train_rating_mat = pd.pivot_table(train, index='user_id', columns='item_id', values='rating')
train_rating_mat.fillna(0,  inplace=True)

メインの実装部

以下がメインの実装部になっています。 user1には参考にしたコードでテスト用に5人選抜していたので、そのまま流用しています。

手順としては、

  • user1に対して全ユーザーのコサイン類以度を算出
  • 算出されたコサイン類以度の上位10人を選抜
  • 上位10人の近接性とuser1の平均評価値を使ってuser1のアイテムに対する評価値を予測
  • 予測結果からトップ10を用いてレコメンドリストを作成
  • 実際にuser1が購入したリストと比較してprecisionで評価

これを実行するとprecisionが22%になり、少し高いような気がしますが参考にしたコード付近のprecisionになりました。

# use five people for evaluation of recommend system
for user1_id in tqdm([1, 100, 233, 666, 888]):
    cos_sim_list = []
    for user2_index in range(rating_arr.shape[1]):
        user1 = rating_arr[:, user1_id-1][:, np.newaxis]
        user2 = rating_arr[:, user2_index][:, np.newaxis]
        two_users_mat = np.concatenate((user1, user2), axis=1)
        two_users_mat = two_users_mat[~np.isnan(two_users_mat).any(axis=1), :]
        # calucalate cosine similarity between user1 and user2
        cos_sim = cosine_similarity(two_users_mat[:, 0], two_users_mat[:, 1], mean_adjustment=True)
        cos_sim_list.append(cos_sim)
        cos_sim_mat = pd.Series(cos_sim_list, index=[i for i in range(1, rating_arr.shape[1] + 1)])

    # use top 10 users of cosine similarity
    top_n = 10
    top_n_sim = cos_sim_mat.sort_values(ascending=False)[1:top_n+1]
    top_n_users = top_n_sim.index

    # test data of user1
    test_user1 = test[test['user_id'] == user1_id].sort_values(by='rating', ascending=False)

    # calculate the prediction of user1
    user1_not_rating = train_rating_mat.iloc[user1_id-1, :]
    user1_not_rating = pd.Series(np.logical_not(user1_not_rating), index=user1_not_rating.index)
    mean_r = train_rating_mat.replace(0, np.nan).mean(axis=1).drop(labels=[user1_id])
    mean_r = mean_r[mean_r.index.isin(top_n_users)]

    not_user1_rating_item = train[train.index.isin(user1_not_rating[user1_not_rating == 1].index)]
    not_user1_rating_item = not_user1_rating_item[not_user1_rating_item['user_id'].isin(top_n_users)]
    not_user1_rating_item.reset_index(inplace=True)

    ra = train_rating_mat.replace(0, np.nan).iloc[0, :].mean()
    bottom_value = np.sum(top_n_sim)
    item_id_list = []
    pred_list = []
    hits = 0

    # recommend top 10 item
    for item_id in not_user1_rating_item['item_id'].unique():
        rating_by_item = not_user1_rating_item[not_user1_rating_item['item_id'] == item_id]
        top_value = np.sum([top_n_sim[uid] * (rating_by_item[rating_by_item['user_id'] == uid]['rating'].values - mean_r[uid]) for uid in rating_by_item['user_id'].values])
        pred = ra + top_value / bottom_value
        item_id_list.append(item_id)
        pred_list.append(pred)

    # check the precision of recommend list
    pred_dict = {'item_id': item_id_list, 'pred': pred_list}
    pred_df = pd.DataFrame.from_dict(pred_dict).sort_values(by='pred', ascending=False).reset_index(drop=True)
    recommend_list = pred_df[:10]['item_id'].values
    purchase_list = test_user1['item_id'].values
    for item_id in recommend_list:
        if item_id in purchase_list:
            hits += 1
    precision_ = hits / 10.0
    precision_list.append(precision_)
print('precision: {:.2f}'.format(sum(precision_list) / len(precision_list)))

最後に

コードは以下においてあります。

とりあえずユーザーベースは理解できた気がするので次はアイテムベースに挑戦します。 間違いなどがあればコメントに書いていただけると幸いです。

参考にした以下の文献に感謝です。

情報推薦システム入門 -理論と実践-

情報推薦システム入門 -理論と実践-

pandasを用いてフラグがついている列が先頭になるように行ごとにシフトする

やりたい処理は、すべての行でフラグ1が先頭にくるようにシフトしたい。

うまいやり方かどうかわからないけど、一応うまくいっているような気がする。

サンプルコード

import numpy as np
from numpy import nan
import pandas as pd
import warnings 
warnings.filterwarnings("ignore") # 実行に関係のない警告を非表示
# 左からゼロパディング
columns_list = ["a_{0:02d}".format(i) for i in range(0, 12)]
value_arr = np.array([[nan, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
                      [nan, nan, 1, 0, 0, nan, nan, nan, nan, nan, nan, nan], 
                      [nan, nan, nan, 1, 0, 0, 0, nan, nan, nan, nan, nan] 
                     ])

時系列のデータを想定したDataFrameを作成。

フラグが1になるカラムが行ごとにバラバラなので、すべての行で先頭のカラムが1になるようにシフトしたい。

df = pd.DataFrame(value_arr, columns=columns_list)
df
a_00 a_01 a_02 a_03 a_04 a_05 a_06 a_07 a_08 a_09 a_10 a_11
0 NaN 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1 NaN NaN 1.0 0.0 0.0 NaN NaN NaN NaN NaN NaN NaN
2 NaN NaN NaN 1.0 0.0 0.0 0.0 NaN NaN NaN NaN NaN

行ごとに処理を行う。

  1. 各行ごとにフラグ1がくるインデックスを取得して、その分だけ左方向にシフトする(作成されるのはSeries)

  2. 作成されたseriesを行方向に結合する

  3. 結合してできたDataFrameをtransposeして元のDataframeと形を揃える。

se_concat = pd.concat([df.ix[i].shift(val) for i, val 
                       in enumerate([-np.where(value_arr == 1)[1][j] 
                                     for j in range(len(df))])], axis=1)
se_concat.T
a_00 a_01 a_02 a_03 a_04 a_05 a_06 a_07 a_08 a_09 a_10 a_11
0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 NaN
1 1.0 0.0 0.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 1.0 0.0 0.0 0.0 NaN NaN NaN NaN NaN NaN NaN NaN

参考にした書籍・サイト

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

<Python, pandas> 縦にずらす。

numpyのarrayで複数の要素が配列内に存在するか判定する

numpy.arrayで複数の要素をリストで渡して、真偽値行列を作成してみる。

使うのは、numpy.in1dで、 指定したarray-likeな要素がある配列内に存在するかどうかを判定して、1次元の真偽値行列を返してくれる。

配列の形を合わせたいなら、reshape(arr.shape)で形を合わせることができる。

In [1]: import numpy as np

In [2]: arr = np.array([[1, 3, 5], [2, 3, 4], [2, 3, 3]])

In [3]: np.in1d(arr, [1, 3, 5])
Out[3]: array([ True,  True,  True, False,  True, False, False,  True,  True])

In [4]: np.in1d(arr, [1, 3, 5]).reshape(arr.shape)
Out[4]: 
array([[ True,  True,  True],
       [False,  True, False],
       [False,  True,  True]])

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

LightGBMをインストールする

LightGBMをubuntu18.04にインストールします。

LightGBMとは

Microsoftが開発した勾配ブースティング(Gradient Boosting)のフレームワーク。勾配ブースティングのフレームワークには、他にXGBoostとかも割と有名なのは知っていましたが、LightGBMは知りませんでした。知ったきっかけは、Kaggleで適当に興味あるコンペのKernelを読み漁っているとよく登場してたため調べてみる気になりました。

勾配ブースティング

そもそもブースティングは、複数の弱学習器を用意して、それぞれの学習器を直列接続するイメージの手法。前の弱学習器で学習した内容を現在の弱学習器に継承しながら学習を進めていきます。

勾配ブースティングは、各ステップごとの弱学習器のでの損失関数の最小化問題に対して勾配降下法を用いるのが由来です。

勾配ブースティングの弱学習器では、決定木が採用されることが多いです。

LightGBMのインストール

LightGBMの公式の通りインストールしていくだけですが、一応手順を示します。 下準備としてsetuptoolsは準備しておきましょう。 今回は、Github経由でのインストールを採用しました。

Cmakeのインストール

Cmakeはソフトウェアのビルド自動化ツールで、Windows,mac, linuxとクロスプラットホームで使えます。 LightGBMのインストール時に必要なのでインストールします。

$ cd ~/Downloads
$  wget https://cmake.org/files/v3.12/cmake-3.12.2.tar.gz
$ tar xvf cmake-3.12.2.tar.gz
$ cd cmake-3.12.2
$ ./configure
$ make
$ sudo make install 
# インストールできているか確認
$ cmake
Usage

  cmake [options] <path-to-source>
  cmake [options] <path-to-existing-build>

Specify a source directory to (re-)generate a build system for it in the
current working directory.  Specify an existing build directory to
re-generate its build system.

Run 'cmake --help' for more information.

cmakeして上の表示がされていれば、インストールできています。 うまく行かない場合は、C++コンパイラがない可能性があるのでgccとかを入れる必要があるかも?

LightGBMのインストール

$ git clone --recursive https://github.com/Microsoft/LightGBM.git
$ cd LightGBM/python-package
$ python setup.py install
# grepしてインストールできているか確認
$ pip list --format=columns | grep -i lightgbm
lightgbm                           2.2.1 

上のように、lightgbm 2.2.1 のようにversion表示されていればインストール成功です。

参考にしたサイト

qiita.com

blog.amedama.jp

はじめてのパターン認識

はじめてのパターン認識

2つのnumpy arrayからDataFrameを作る

2つの同じ長さのnumpy arrayを用いてそれらを列に持つpandasのDataFrameを作成します。

In [1]: import numpy as np

In [2]: import pandas as pd

In [3]: arr1 = np.array(["a", "b", "c", "d"])

In [4]: arr2 = np.array([1, 2, 3, 4])

In [5]: dict_ = dict(zip(arr1, arr2))

In [6]: df = pd.DataFrame(list(dict_.items()), columns=["col1", "col2"])

In [7]: df
Out[7]: 
  col1  col2
0    a     1
1    b     2
2    c     3
3    d     4

手順は、

  1. 1つのarrayをkey、もう一つのarrayをvalueとする辞書を作成
  2. 作成した辞書をlist.items()でkeyとvalueのarrayのタプルを1つの値ごとに格納したリストに変換
  3. 2.で作成したリストをDataFrameに変換

という流れです。

参考にしたサイトは以下のサイトです。

stackoverflow.com

PandasのDataFrameから特定の値を持つ行を削除する

言われてみたら簡単なんだけど、意外に思いつかなかった。

DataFrameの特定の行に含まれている値を指定して、それ以外を抽出するイメージ。

In [1]: import pandas as pd

In [2]: df = pd.DataFrame([[1, 2], [2, 3], [3, 4]], columns=["a", "b"])

In [3]: df
Out[3]: 
   a  b
0  1  2
1  2  3
2  3  4

In [4]: df[df.a != 2]
Out[4]: 
   a  b
0  1  2
2  3  4

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

pythonで地図上の2時点間の位置関係を求める;

2時点間の緯度、経度、高度が与えられたときに簡易的に2時点間のローカルな位置関係を求めます。 緯度・経度・高度といった情報ではローカルな2時点間の関係性が分かりづらく扱いづらいので、変換します。 具体的には、2時点間の距離とある時点からもう一時点を見たときの方位角、仰角を求めることを目標とします。

以下、国土地理院のサイトと独立行政法人電子航法研究所の資料を参考に作成しています。

https://www.enri.go.jp/~fks442/K_MUSEN/1st/1st060428rev2.pdf

地球上の2時点間の距離

高校などで扱った平面世界での2時点間の距離は、簡単に求めることができたと思いますが、地球上の2時点間ではそれほど簡単には求めることができません。 なぜなら、地球が球体であり、局面上にある2時点間の距離を求めることになるからです。

ジオイドとは

地球は、時点による遠心力の影響で、赤道方向が少し膨らんだ楕円体のような形をしています。 測地学では、世界の界面の平均位置に最も近い「重力の等ポテンシャル面」を「ジオイド」と定めて、これを地球の形状ということにしています。標高はジオイドを基準として定められており、標高はジオイドから測った高さになります。

GPSを用いて標高を求めるには

現在広く活用されているGPSでは、緯度・経度・高度といった幾何学的な位置を求めることができますが、重力を考慮した設計になっていないため標高を直接求めることができません。GPSを用いて標高を求めるためには、ジオイド高がj必要になります。

標高 = 楕円体高 ー ジオイド

ジオイド高は国土地理院が提供している計算サイトで計算することができます。

ジオイド計算

WGS84

WGS84とは、GPSの基準座標系で、直交座標系です。これはECEFと呼ばれ、地球の重心を原点として地球の自転軸の北極方向をZ軸、Z軸に垂直にグリニッジ子午線の方向をX軸として、これらの軸と直行するように右手系でY軸とした直交座標系です。

ECEFから東をX, 北をY、上空をZとした地平直交座標(ENU)に変換します。 ECEFからENUへの変換は、行列による回転と原点移動で実現できます。詳しくは上記で示した電子航法研究所の資料をご確認ください。

githubにサンプルコードを提示してあります。coordinate_test.pyの9, 10行目で取り扱いたい2時点の緯度・経度・ジオイド高のサンプルデータを与えています。結果表示では、2時点間の距離およびある時点(コードではRx)からみたもう1時点(コードではSat)の方位角、仰角を求めることができます。

github.com