類似本検索システムを作りたい
ということで、本を検索すると類似している本のリストを出力するサイトを試作してみました。
https://bookrecommendst.herokuapp.com/
(herokuを他の用途で使うまではアクセス可能な予定です。)
こんな感じで、ある検索した本に対して類似度の高い本を出力してくれます。
作成にあたり、こちらのサイトを参考にさせていただいています。
実施したこと
1. データの取得
ブクログ様のサイトの情報を使用しています。 ブクログではユーザ毎にどの本を読んでいるかという情報を閲覧できるようになっているので、その情報をスクレイピングにより取得しました。 全くのランダムにユーザをピックアップするということができませんでしたので、人気タグの1ページ目にあるタグが付与されている人のIDを6万人程取得し、その各ユーザに対し登録されている本の情報を取得しました。
2. データの前処理、整形
取得した情報より、下記のようなユーザ*本の行列を作ります。
本A | 本B | 本C | ... | |
---|---|---|---|---|
ユーザ1 | 0 | 0 | 1 | ... |
ユーザ2 | 1 | 1 | 0 | ... |
... | ... | ... | ... |
スパースな行列になるので、メモリの爆発を防ぐためにscipyのlilmatrixを使用しています。 また、本の登録者が10冊未満のユーザおよび、全体で10冊未満しか登録されていない本についてはあらかじめ除外をしています。これは、レコメンドの安定性の向上および、特定のユーザが結果に大きな影響を与えるような結果を防ぐための処理になります。 この除外処理により、ユーザ数は6万→5万、ユニークな本の冊数は38万→2万となりました。 本については想定以上に多様性が大きく、今回の取得したデータ量だと結構対象外となるものが多くなってしまいました。
3. 次元削減
経験的にはスパース行列のまま類似度計算を行うよりも、次元削減をした方が意味のある特徴量になることが多いです。また、計算時間の削減にも寄与するため次元削減を行いました。 次元削減手法はSVDとNMFを選択できるようにしています。選定理由はスパース行列に対応している次元削減手法だったからです。 この処理により列数を2万→100に削減しています。
4. 類似度計算
コサイン類似度を計算します。 最初は全入力パターンを事前に計算して結果を持っておこうと思っていましたが、 下記サイトを参考に行列計算すればかなり高速にできることがわかったので、入力毎に毎回計算しています。
https://stackoverflow.com/questions/41905029/create-cosine-similarity-matrix-numpy
5. wepアプリ実装
上記の処理をstreamlitを用いてwebアプリ化し、それをherokuに実装しました。 検索は単純な部分一致検索です。herokuは無料枠です。
サンプルの確認
いくつか自分で検索してみた結果を示します。すべてSVDでの結果になります。
- know 野崎まど
順当にSFっぽいものが上位にきている気がします。 ホラーがちらほらあるのが少し気になりますね。
- ファスト&スロー(上) ダニエル・カーネマン
ちゃんと下巻がトップにでてきました。 ラインナップを見る限り、割と妥当という印象を受けます。
もう少し同名タイトルの別巻が並ぶかと思いましたがそうでもなかったです。 こういうサイトで漫画をきっちり登録している人はあまり多くない印象なので、 なんとなくですが漫画の入力に対する結果はよくなさそうです。
- 完全教祖マニュアル 架神恭介
もっと尖った結果を期待しましたが、 普通にビジネスマンが読むような本が出てきてしまいました。 こういう本を読んでいるような人はそもそも読書量が多く雑食ぎみなので、 共通項をとるとこんな感じになってしまうのでしょうか。
まとめ
ブクログで取得した情報を元に、本を検索すると類似している本のリストを出力するサイトを試作しました。 今回取得したデータ量はそれほど多くはなかったので、対象外となってしまった本も多いですし、入力に対して疑問符がつく結果となるケースもそれなりにありました。 ただ、自分が読んだことある本に対して検索をかけられるのは個人的には楽しめました。
どうしても今回のようなやりかただとみんなが共通して読んでいるものが上位に出てきてしまうので、 よりマイナーで尖っている本を探したい場合は、あえて広く皆に読まれている本を除外しておくような、 また違ったアプローチが必要になりそうかなと思いました。
今回使用したコードは下記になります。
GCPでデータ収集環境をつくる
はじめに
データ分析をするにはデータが必要ですが、常に欲しいデータが存在するとは限らず、 時には自分で取得・保存する必要が出てくると思います。 ここではGCPを用いてデータ収集環境を構築してみた備忘録を記述します。
実施概要
apiを使って日次でvtuberのデータを集めます。 全体の実行イメージ図は以下のような形です。
自宅にサーバがあれば楽なんですが、残念ながら持っていないのでGCPを使いました。 cloud schedulerを使ってvmインスタンスを起動、 日次でcronを用いてスクリプトを実行しデータを取得、 取得結果をgdriveとbigqueryに転送し、slackに結果を通知するという仕組みにしました。 cloud schedulerを使っているのは課金対策で、vmインスタンスはスクリプト実行後にシャットダウンを行います。
vmインスタンスでなくgcp functionを使う手も考えましたが、 個人的にインスタンスを使用した方が取り回し安かったのでそうしました。 動かすのは10分未満/日で使用しているインスタンスはn1-standard-1なので、 今回作った構成だとサーバ代は月100円未満です。 (vmインスタンス以外は無料枠)
データの取得
このようなvtuberのtwitterIDとyoutube channelIDのリストを元にapiを利用します。
とりあえず今回は所属がにじさんじとホロライブの人だけを対象としています。
twitterデータ取得
tweepyを使いました。詳細は公式参照。認証周りについては割愛します。 https://github.com/tweepy/tweepy
取得したいのはユーザ毎の統計情報なので以下のような感じ。
twitter_auth = tweepy.OAuthHandler(consumer_key, consumer_secret) twitter_auth.set_access_token(access_token_key, access_token_secret) api = tweepy.API(auth) results = api.get_user(id=tweetid) #任意のtweetidを入力 created_at=results.created_at #アカウント作成日 followers_count=results.followers_count #フォロワー数 friends_count=results.friends_count #フォロー数 avourites_count=results.favourites_count #いいね数 statuses_count=results.statuses_count #つぶやき数 listed_count=results.listed_count #リスト数
時間当たりのapiリクエスト数には制限があるので、以下を参考に調整します。 https://developer.twitter.com/ja/docs/basics/rate-limits
tweepyのコードで確認すると、今回のケースだとusers/showを使っているようなので900request/15minutesに収まるようリクエスト数を調整します。
youtubeデータ取得
Pythonコードサンプル:https://developers.google.com/youtube/v3/code_samples/python?hl=ja
ChannelsのAPI:https://developers.google.com/youtube/v3/docs/channels?hl=ja
こちらも認証周りについては割愛します。
実装はこんな感じ
youtube_api_key=data["YOUTUBE_API_KEY"] youtube_auth = build('youtube', 'v3', developerKey=youtube_api_key) search_response = auth.channels().list( part='statistics,snippet', id=channel_id, ).execute() viewCount=search_response["items"][0]["statistics"]["viewCount"] #再生数 subscriberCount=search_response["items"][0]["statistics"]["subscriberCount"] #登録者 videoCount=search_response["items"][0]["statistics"]["videoCount"] #視聴数 publishedAt=search_response["items"][0]["snippet"]["publishedAt"] #作成日
youtube apiを使う場合もリクエスト数の制限を考慮する必要があります。 twitter apiとは異なり、クエリ種別毎にコストの概念があります。 こちらでコストを参照すると、 今回用いるchannel,list,statistics,snippetだとコストが5であることが確認できます。 デフォルトだと1日のコスト上限は10000みたいなので、1日に取得できるユーザ数は2000が上限になります。
gcpインスタンス自動起動
cloud schedulerの公式にインスタンス起動のサンプルがあるのでそこからインスタンス名とスケジュール実行時間のみ変更。(jsなのでコードはあまり理解できていない)
gdriveへのアップロード
pydriveを使います。
参考:https://note.nkmk.me/python-pydrive-download-upload-delete/
コード例
#cronで実行するためフルパスで指定 gauth = GoogleAuth(settings_file='/home/rmizuta/get_tweet/settings.yaml') store = oauth2client.file.Storage('/home/rmizuta/get_tweet/credentials.json') credentials = store.get() gauth.credentials = credentials gauth.CommandLineAuth() drive = GoogleDrive(gauth) #フォルダ名からidを取得 folder_id = drive.ListFile({'q': 'title = "vtuber_stats"'}).GetList()[0]['id'] f = drive.CreateFile({"parents": [{"id": folder_id}]}) f.SetContentFile(filename) f['title'] = filename.split('/')[-1] f.Upload()
slackへの通知
やるだけ。異常があったときにわかればいいので処理した行列数と欠損数を通知。
bigqueryへのデータ転送
公式にpythonのapiサンプルがあるのでそれを参考にします。
https://googleapis.dev/python/bigquery/latest/index.html https://cloud.google.com/bigquery/docs/tables?hl=ja
カラム定義等のテーブルの作成は最初の取得データをGUIでインポートして自動定義し、2日目以降の取得データは下記のスクリプトで作成済みのテーブルに追加していきます。 https://cloud.google.com/bigquery/docs/loading-data-local?hl=ja#console
client = bigquery.Client() table_ref = client.dataset('vtuber_stats').table('table1') job_config = bigquery.LoadJobConfig() job_config.write_disposition = bigquery.WriteDisposition.WRITE_APPEND #WRITE_APPEND:追記 load_job = client.load_table_from_dataframe( df, table_ref, job_config=job_config ) # API request load_job.result() # Waits for table load to complete. destination_table = client.get_table(table_ref)
bigqueryにデータが入ってきているか確認します。
うまく動いているようです。 (日付けが古いのはこの文章を書くのを寝かせていたため)
まとめ
GCPを用いてデータ取得から保存を行う環境を構築しました。 最近のクラウドは機能も多く無料枠もそれなりにあるので、覚えれば結構使い勝手が良いと思います。 データが貯まったら可視化とか分析を行いたいですね。
pca+kmeansについての雑実験
はじめに
列数が多いデータセットに対してクラスタリングを行う場合にPCAで列数を次元削減してからクラスタリングをするという手法があるらしいです。 確かにPCA等で列の次元削減を行うことでノイズ成分を落とせるので、うまくいけば重要となる特徴だけを用いたクラスタリングができそうです。 近い考え方で言語処理系のデータを扱うときにBag of Words+TFIDFで作った特徴量をSVDで次元削減したデータを学習器への入力とするという手法はありますね。
下記文献の実験では上記の手法は有効となっていますが、ちょっと中身理解しきれていない部分もあるので実験ベースで確認をしていきたいと思います。 http://ranger.uta.edu/~chqding/papers/KmeansPCA1.pdf
実験
MovieLens 1M Datasetを用いて実験を行いました。 このデータセットは約6000のユーザが約4000の映画に対して5段階評価をつけたデータセットになります。
6000のユーザを映画の嗜好性が近いグループでまとめるために クラスタリングで10のグループに分割するというタスクを考えます。
異なる特徴量を用いた場合のクラスタリングの評価をどう比較しようかという問題があります。 同じ特徴量であればエルボー法で使うような残差平方和が使えますが、今回は比較したい対象同士で特徴量の数が異なります。 なので今回は映画の嗜好でクラスタリングをした場合、年齢や性別もそれなりに別れるはずだという仮定のもと、 クラスタリングには用いていないユーザに紐づく変数の性別(2カテゴリ)と年代(7カテゴリ)をダミー化した特徴の残差平方和で評価を行うことにします。
列方向の次元削減を全く行わない場合と、PCAで削減する特徴量数を変化させた場合の試行を行った結果、それぞれの残差平方和と所要時間は以下のようになりました。
手法 | 特徴量数 | 残差平方和 | 所要時間(s) |
---|---|---|---|
kmeans | 3706 | 7002 | 32.1 |
PCA+kmeans | 1000 | 7003 | 18.5 |
PCA+kmeans | 500 | 6972 | 10.1 |
PCA+kmeans | 100 | 7021 | 4.1 |
PCA+kmeans | 50 | 6986 | 2.6 |
PCA+kmeans | 10 | 6978 | 2.1 |
実験コードは下記になります。
考察
今回の結果だとPCAで次元を削減してもしなくても残差平方和はあまり変わりませんでした。 評価指標があまりよくなかったのかもしれませんが、このへんはデータ依存性が強い気がしています。 ただ、処理時間は削減した次元数に応じて小さくなるので、データ量が巨大で処理時間を気にする必要がある場合には有効に働く可能性がありそうです。
あてはまりのよい確率分布を探したい
はじめに
データを眺めていると、ある分布に対してそれが正規分布に従うのか、対数正規分布か、それともガンマ分布の方が近いのか?、というようにどの分布の当てはまりがよいかが気になることがあると思います。
これを確認する方法を探してみたところ、scipy.statsを使えばできそうだったのと、fitterというライブラリもあったので、それらを試してみた結果を記述します。
実験
scipyを使う
実装はnumpy - Fitting empirical distribution to theoretical ones with Scipy (Python)? - Stack Overflowを少しだけ修正したものです。入力に対してscipy.statsに登録されているすべての確率分布のパラメータを最尤推定した結果の平均二乗誤差を比較することで最もあてはまりのよい分布を求めます。
scipyには80以上の確率分布が存在しているのですが、全部回してしまうと100行程度のデータでも数分かかってしまいますし、知らない分布が選ばれたところで仕方がないので、ある程度選択しておくのが現実的かなと思います。今回は正規分布、一様分布、ガンマ分布、レイリー分布の4つの分布で比較を行います。
from sklearn.datasets import load_iris import pandas as pd import scipy.stats as st import matplotlib.pyplot as plt #サンプルデータの読み込み iris = load_iris() data=pd.DataFrame(iris.data, columns=iris.feature_names) #使用する分布 #正規分布、一様分布、ガンマ分布、レイリー分布 distributions=[st.norm,st.uniform,st.gamma,st.rayleigh] result=[] for distribution in distributions: y, x = np.histogram(data.iloc[:,1], bins=20, density=True) #xの値はヒストグラムの左端の値なので中心点に修正 x = (x + np.roll(x, -1))[:-1] / 2.0 params = distribution.fit(data.iloc[:,1]) # Separate parts of parameters arg = params[:-2] loc = params[-2] scale = params[-1] # Calculate fitted PDF and error with fit in distribution pdf = distribution.pdf(x, loc=loc, scale=scale, *arg) sse = np.sum(np.power(y - pdf, 2.0)) result.append((pdf,sse)) df=pd.DataFrame() df["name"]=[i.name for i in distributions] df["sumsquare_error"]=[i[1] for i in result] df
今回の分布の中では正規分布のあてはまりが最も良いようです。
for i in range(len(result)): pd.Series(result[i][0], x).plot(label=distributions[i].name) plt.hist(data.iloc[:,1],density=True,alpha=0.4,bins=20) plt.legend() plt.show()
正規分布、ガンマ分布は元の分布に対してあてはまりが比較的良いことが確認できます。
fitterを使う
上に書いたものをライブラリ化したような感じです。
from fitter import Fitter f = Fitter(data.iloc[:,1],distributions=['gamma', 'rayleigh', 'uniform','norm'],bins=20) f.fit() f.summary()
中でやっていることは同じなので結果は先ほどと同じになります。
気をつけたいこととしては、あてはまりの良い分布はbinの数によって変化することかなと思います。 fitterのbinのデフォルト値は100なのですが、irisデータはサンプルサイズが150と小さくデフォルトのbinを用いると 結果がおかしくなるため、今回の実験ではbinを20としています。
まとめ
任意のデータに対し、当てはまりの良い確率分布を導出する手法の調査を行いました。 確率分布の仮定が必要な統計モデリングを行う際や、 分布からデータ生成の仕組みなどを推測するのに役立つかもしれません。
クロス集計表とナイーブベイズの対応についてのメモ
はじめに
ある日クロス集計表とナイーブベイズを眺めていて、 そもそもこれらはどういう関係性だっけ?と思ったのでその思考メモです。
試行
以下のような所属クラス×性別×利き手の3軸のクロス集計表があるとし、 これを用いて性別と利き手がわかっている人に対して所属クラスを予測したいケースを考えます。
- 理系クラス
右利き | 左利き | |
---|---|---|
男性 | 30 | 5 |
女性 | 20 | 5 |
- 文系クラス
右利き | 左利き | |
---|---|---|
男性 | 10 | 5 |
女性 | 15 | 10 |
これを全体を足して1になるように正規化すると下記のようになり、 これは同時確率と考えることができます。
- 理系クラス
右利き | 左利き | |
---|---|---|
男性 | 0.30 | 0.05 |
女性 | 0.10 | 0.15 |
- 文系クラス
右利き | 左利き | |
---|---|---|
男性 | 0.10 | 0.05 |
女性 | 0.20 | 0.05 |
Yを所属クラス、X1を性別、X2を利き手だとすると 判別対象が女性・左利きの場合、
となるので文系クラスと予測した方が確率が高くなります。
次に、この同時確率をナイーブベイズに変換していきます。 乗法定理より
これにナイーブベイズの条件である説明変数の独立性を考えると、 X2はX1に依存しなくなるので
となり、これはナイーブベイズの式となります。 この式を元にどちらのクラスに在籍しているかを判別すると、
- Y=理系の場合
- Y=文系の場合
とこちらの場合は理系クラスと文系クラスの所属確率が等しくなります。
3軸の集計表ベースだと交互作用は考慮できますが、利き手のみや性別のみの作用を考えることができません。 逆にナイーブベイズだと特徴量の独立の仮定を置いているので、交互作用は考慮できませんが、各要素毎の影響は考慮することができます。あと軸が増えると1地点あたりのサンプルサイズが小さくなる問題も起きにくそうです。
その他、選択肢としては重回帰を使うケースも考えられそうです。 この場合は交互作用項は追加できるが多重共線性が気になるという感じでしょうか。 参考資料2によるとサンプルサイズが小さい時はナイーブベイズに優位性があるようですが、結構データ依存性が高そうです。 (結論しか読めてません)
参考資料
1.http://www3.u-toyama.ac.jp/kkarato/2016/statistics/handout/statistics-2016-23-0705.pdf
推論時の入力に未知の欠損値がある場合のlightgbmの挙動の確認
はじめに
以前にlightgbmは入力に欠損値があってもうまく学習してくれるという記事を書いたのですが、 これは学習時に欠損が存在している場合の話でした。
現実の問題を考えると、学習時とそのモデルを使った推論時では時系列の違いや環境変化の影響 により取得していたデータが欠落したりするケースもあると思います。 このように学習データには存在しなかったが、推論時の入力で欠損が入ってしまうパターンが発生した時にモデルがどのような挙動をとるかが気になりましたので、 調査・実験をしました。
参考文献
こちらでは、学習時と推論時でカテゴリカル変数列が異なる(推論時に未知のカテゴリが発生する)場合について書かれており、 カテゴリ変数であることを明示しておけばis, is notで判定されるため基本的には問題ないとしていました。
When zero_as_missing=false (default), the unshown values in sparse matrices (and LightSVM) are treated as zeros.
公式によると上記のように書いてあるので、欠損は0と扱われるっぽい?
実験
sklearnにあるボストンの住宅価格のデータを用い、 欠損なしで学習させたモデルに数値型の変数に欠損がある入力を使って推論した場合の挙動を確認してみます。
このデータセットには欠損はなく、カラムはすべて数値型になります。
boston = load_boston() X=pd.DataFrame(boston.data, columns=boston.feature_names) y=boston.target X.head()
X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.33, random_state=42) params={ 'objective': 'regression', 'random_state' : 1, "metric": "rmse", } dtrain = lgb.Dataset(X_train, label=y_train) dvalid = lgb.Dataset(X_valid, label=y_valid) bst = lgb.train(params, dtrain, num_boost_round=1000,valid_sets=[dtrain, dvalid],early_stopping_rounds=50,verbose_eval=100) #AGEをすべて0にした入力データ作成 X_valid_zero=X_valid.copy() X_valid_zero["AGE"]=0 #AGEを全て欠損にした入力データ作成 X_valid_nan=X_valid.copy() X_valid_nan["AGE"]=np.nan checkdf=pd.DataFrame() checkdf["normal"]=bst.predict(X_valid) checkdf["zero"]=bst.predict(X_valid_zero) checkdf["nan"]=bst.predict(X_valid_nan) checkdf.head()
すべて0としたケースとすべて欠損としたケースで予測値が一致しました。 一応他の変数についても同様のことを行いましたが挙動は同じでした。 学習時に欠損がなかった場合、予測時に投入するデータの欠損は0として扱われるようです。
全文コードは以下になります。
まとめ
欠損のないデータでモデルを作成した場合、欠損のあるデータを予測時に投入すると lightgbmはその欠損は0として扱うことが確認できました。 lightgbmだとこういうケースでもエラーにならずに出力してしまうので注意が必要ですね。 対応策としては、モデルの再作成や推論時の入力データの欠損値補完が良いかなと思いました。
テニスにおける疲労度の影響の定量化(spoana5の内容+α)
はじめに
「プロテニスにおいて疲れが勝敗に与える影響を定量化してみる」 という題目で下記イベントでLTしてきました。
使用した資料は以下になります。
プロテニスにおいて疲れが勝敗に与える影響を定量化してみる - Speaker Deck
内容の余談
行間的なものを何か書こうと思っていましたが正直あまり追記するようなこともなかったので使用したコードだけ載せておきます。 後述の追加検証のコードもこちらになります。
スライドに書いていない内容としては以下のような感じです。
レーティングと試合の勝敗のデータは別のデータソースから取得し選手名をキーとして紐づけているが、"Stan Wawrinka"と"Stanislas Wawrinka"のように 選手名の表記が異なっている人が何名かいたため、対応がわかる人については名前を手動修正、わからなかった人は行を捨てている。
ロジスティック回帰に入れる際、正解が1だけになってしまうのを防ぐため、勝敗を反転させたデータを追加している
疲労度を算出する際、選手・大会ごとに値を集計する必要があるため前処理が少し面倒
追加検証
twitterのコメントを読んでいたら、 汎化性能はどうなのか?とか解釈性はpdpという選択肢もあるというようなコメントを見かけたので、 この辺りの追試もしてみました。
交差検証
汎化性の確認のためロジスティック回帰で交差検証をしました。 データをkfoldで4分割し、そのオッズ比の検証をしました。
試行 | オッズ比 |
---|---|
fold1 | 0.9957 |
fold2 | 0.9969 |
fold3 | 0.9902 |
fold4 | 0.9953 |
全件 | 0.9946 |
それなりに汎化性能はあるのではないでしょうか。 全件データで出した95%信頼区間は0.9913~0.9979だったので、 信頼区間ベースでも3/4はこの枠内に収まっています。
lightgbm + pdp
次に、lightgbm+pdp(partial dependence plot)で結果の解釈をしてみます。
pbpは機械学習の解釈方法の一つです。まず機械学習モデルで学習を行い、その後入力データ全件に対し影響を見たい変数を入力データ内で一律で変化させながら予測値の平均値を算出することで、 変数が予測値に与える影響を確認する方法です。各変数は独立であるという前提をおいている点には注意が必要です。
参考資料:
ロジスティック回帰のときに使った入力データと同じものを lightgbmで学習させ、そのpdpを算出します。 4foldで行い、そのout of foldで出した結果が以下になります。
横軸が疲労度の差、縦軸が勝率に与える影響です。 一番左の点が基準点扱いとなっているため、疲労度の差が-100以下の場合はない場合と比較して勝率が0.08くらい上がり、 疲労度の差が250くらいあると勝率が0.1くらい下がるような結果です。
入力データには勝者と敗者を反転させたものを入れているので左右対称にならないのは違和感がありますが、 サンプルサイズが小さいのでこのあたりは仕方のない部分なのかなとは思います。
今回のデータの疲労度の差は実は以下のようなヒストグラムを描いています。
入力データのサンプルサイズは964あるのですが、 そのうち疲労度の差が100を超えるものは48しかありません。 疲労度の差が出にくい1回戦が全体の1/2、2回戦が全体の1/4のデータを占めるためこのような形になってしまいます。 回帰木は線型性を仮定しないので、疲労度の絶対値が大きい部分の予測値はブレが生じやすい結果になっているのではないかと思います。とはいえ全体の傾向をつかむという目的の場合はこのくらいのサンプルサイズの時でも使用できるかなと思いました。
まとめ
解釈系の話の妥当性確認は難しい。