R Advent Calendar 2017 rvestを用いてポケモンデータをスクレイピング&分析してみた

R Advent Calendar 2017の11日目を担当するMr_Sakaueです。
今回はrvestパッケージを用いて、友人がハマっているポケモンの情報を集めてみようと思います。
もっとも、業務でWebスクレイピングする際はPythonでBeautifulSoupやSeleniumを使うことがほとんどなのですが、たまにはRでやってみようと思います。

目次
・やりたいこと
・rvestについて
・データの取得と集計と可視化と分析
・まとめ
・参考情報

やりたいこと

今回はポケモンたちのデータを集めた上で、以下の内容を行いたいと思います。

  • ポケモンのサイトから種族値を取得
  • ポケモンの種族値を標準化して再度ランキング
  • ポケモンのレア度や経験値に関する情報を取得
  • レア度や経験値と相関しそうな種族値を探る

今回扱った全てのコードはこちらに載せております。
https://github.com/KamonohashiPerry/r_advent_calendar_2017/tree/master

※種族値はゲームにおける隠しパラメータとして設定されている、ポケモンの能力値とされている。

rvestについて

rvestはRでWebスクレイピングを簡単に行えるパッケージです。ここでの説明は不要に思われますが、今回はread_html()、html_nodes()、html_text()、html_attr()の4つ関数を用いました。

基本的に以下の3ステップでWebの情報を取得することができます。

  • STEP1
    read_html()でHTMLからソースコードを取得する。(Pythonでいう、requestとBeautifulSoup)
  • STEP2
    html_nodes()でソースコードから指定した要素を抽出する。(PythonでいうところのfindAll)
  • STEP3
    html_text()やhtml_attr()で抽出した要素からテキストやリンクを抽出する。(Pythonでいうところのget(‘href’)など)

データの取得と集計と可視化

検索エンジンで検索してだいたい1位のサイトがあったので、そちらのWebサイトに載っているポケモンの種族値の一覧をスクレイピング対象とさせていただきます。

  • ポケモンのサイトから種族値を取得

以上のコードを実行すれば、こんな感じでポケモンの種族値一覧を得る事ができます。

とりあえず、種族値合計(Total Tribal Value 以下、TTV)のランキングの上位を確認してみます。知らないんですが、メガミュウツーとかいうイカつそうなポケモンが上位にいるようです。昭和の世代には縁のなさそうなポケモンばかりですねぇ。

■TTVランキング

取得した種族値を項目別に集計したり、Boxプロットを描いてみます。どうやら、攻撃の平均が高く、ヒットポイントや素早さの平均は低いようです。

■種族値のサマリー

■種族値のBoxプロット

  • ポケモンの種族値を標準化して再度ランキング

さて、攻撃の平均が高かったり、ヒットポイントと素早さの平均が低かったりしたので、各々の項目を標準化した上で、再度ランキングを作ってみたいと思います。

■標準化した種族値のサマリー

平均0、分散1にできているようです。

■標準化した種族値のBoxプロット

他よりも低かったヒットポイントと、高かった攻撃がならされていることが確認できます。

■標準化前後でのTTVランキングのギャップが大きかったものをピックアップ

ラッキーが144位ほど出世しています。攻撃が低く、ヒットポイントの高いラッキーが標準化により優遇されるようになったと考える事ができます。ポケモン大会の上位ランカーである後輩社員もラッキーは手強いですと言っていたのでまんざらでもないのでしょう。

  • ポケモンのレア度や経験値に関する情報を取得

今回のサイトには、個別にポケモン別のページが用意されており、そちらから、ゲットしやすさや経験値に関する情報を抽出します。

以上のコードを実行すれば、やや時間がかかりますが、全ポケモンのゲットしやすさや経験値のデータを抽出する事ができます。それらの情報がゲットできたら、まずは可視化します。

■ゲットしやすさのヒストグラム

ゲットのしやすさは、小さいほど捕まえる難易度が高くなっています。難易度の高いポケモンである0が多過ぎるので、このデータは欠損値が0になっているのではないかと疑われます。

■経験値のヒストグラム

経験値は、レベル100になるまでに要する経験値をさしています。ほとんどが100万程度となっているようです。

■ゲットしやすさと標準化TTVの散布図

やはり、ゲットしやすさに関してはデータに不備があるようで、コラッタ(アローラの姿)のような雑魚ポケのゲットのしやすさが0だったり、伝説のポケモンであるネクロズマが255だったりします。ただ、上限と下限のデータを間引けば右下がりの傾向が見られそうです。

■経験値と標準化TTVの散布図

経験値が多く必要にも関わらず、TTVが低い集団があります。どうやらこの集団に属するのは、「キノガッサ」・「マクノシタ」・「イルミーゼ」・「ゴクリン」・「シザリガー」などで、一回しか進化しないポケモンのようです。これらのポケモンは育てにくく、TTVの低い、コスパの悪そうなポケモンと考えることができるのではないでしょうか。(技や特性によってはバリューあるかもしれませんが。)

  • レア度や経験値と相関しそうな種族値を探る

先ほどのレア度に関しては、データがおかしそうだったので、レア度0と255に関しては除外してみます。

■ゲットしやすさと標準化TTVの散布図

やはり除外する事で、理想的な右下がりの傾向を示す散布図が得られたと思います。
さて、各種族値がレア度にどれだけ相関しているのかを分析したいのですが、その前にレア度を表す二項変数を作成します。

■ゲットしやすさが50以下であれば1、それ以外を0にする変数を作成

続いて、各種族値を説明変数として、レア度を目的変数としたロジスティック回帰モデルの推定をrstanで実行させます。

■stanコード

■rstanでロジスティック回帰を行い、推定結果を可視化するコード

■MCMCのシミュレーション結果のトレースプロット

どうやら収束してそうです。

■ロジスティック回帰の推定結果

見にくいので、推定結果を松浦さんの「StanとRでベイズ統計モデリング」にあるコードを用いて可視化します。

■推定結果の可視化

どうやら、0を含まない係数について見てみると、b3(攻撃)、b5(特殊攻撃)、b6(特殊防御)が高いほど、レア度が増す傾向があるようです。珍しいポケモンは攻撃が強いという傾向があると言えるのではないでしょうか。

まとめ

  • rvestは簡単にスクレイピングできて便利。
  • ポケモンデータは色々整備されてそうで今後も分析したら面白そう。
  • 珍しいポケモンは「攻撃」、「特殊攻撃」、「特殊防御」が高い傾向がある。
  • 経験値が必要なのにTTVの低い、コスパの悪そうなポケモンたちがいる。

それでは、どうか良い年末をお過ごし下さい!
メリークリスマス!

参考情報

データサイエンティストのための最新知識と実践 Rではじめよう! [モダン]なデータ分析
StanとRでベイズ統計モデリング (Wonderful R)
【R言語】rvestパッケージによるウェブスクレイピング その2
Receiving NAs when scraping links (href) with rvest

階層ベイズモデルの直帰率分析への適用 with rstan

松浦さんの『StanとRでベイズ統計モデリング』の8章の階層ベイズがすごくわかりやすいなぁと思いつつも、自分の持っているデータで試していなかったので、これを機に実践してみようと思います。
やや変数を追加しているくらいで大した変更点はありませんが、題材としては当ブログのアクセスログにおける直帰率に関するデータで、どのような要素が直帰率に影響を与えるのかを分析します。

目次
・モデル概要
・前処理
・推定
・結果(非階層モデルとの比較)
・参考文献

モデル概要

モデルは8章のロジスティック回帰の階層モデルに一部変数を追加していますが、ほぼそのままです。記事ごとのパラメータやリファラーごとのパラメータを想定しています。
Nは記事数でnはそのインデックス、Cはリファラーの数でcはそのインデックス、Iはログとして残っているセッションの数でiはそのインデックスとなっています。hatebuは記事のはてぶ数、stringlineは記事の行数、holidayは休日or祝日ダミー変数、daytimeは12:00~18:00なら1をとるダミー変数、revisitedは再訪問ユーザーなら1を取るダミー変数となっています。記事ごと・リファラーごとに直帰のしやすさが違う(パラメータが従う正規分布のパラメータがそれぞれ異なる)という仮定のもとに立ったモデルとなります。

$$x[i] = b_{1} + x_{記事}[記事ID[i]] \\ + x_{リファラー}[リファラーID[i]] + x_{セッション}[i] $$

$$q[i] = inverselogit(x[i]) $$

$$Y[i] \sim Bernoulli(q[i]) $$

$$x_{記事}[n] = b_{2}hatebu + b_{3}stringline[n] \\ + b_{記事間の差}[n] $$

$$b_{記事間の差}[n] \sim Normal(0, \sigma_{記事番号}) $$

$$x_{リファラー}[c] = b_{リファラー間の差}[c] $$

$$b_{リファラー間の差}[c] \sim Normal(0, \sigma_{リファラー番号}) $$

$$x_{セッション}[i] = b_{4}holiday[i]  + b_{5}divice[i] \\ + b_{6}daytime[i] + b_{7}revisited[i] $$

前処理

GAのAPIからデータを取得して1セッション1記事になるようにデータを作成しています。数ヶ月で25000件ほどデータがあったのですが、計算に時間がかかるので、データ数を2400件くらいにサンプリングしています。

推定

stanコードはこちらになります。

rstanを用いたstan実行用のRコードです。ヴァイオリンプロットで主要な係数の分布を見る処理も書かれています。

結果

係数を見る限りは、符号の向きが確かなのはb5(PCダミー)とb7(再訪問ユーザーダミー)なので、PCの方が直帰しにくく、再訪問ユーザーの方が直帰しにくいという傾向があると考えることができます。

教科書ではAUCを非階層モデルと比較していましたので、比較してみようと思います。
AUCの計算を行うためのコードもGithubに載っていましたのでそちらを使います。

80%が良いとされているAUCには程遠いですが、記事やリファラーごとの差を考慮しない非階層のものよりもAUCが高いと言えます。
ちなみに、教科書の例のAUCは80%ほどでした。
Webマーケのデータ分析においてロジットは汎用性が高いで、今回のコードを土台に色々と業務で試していこうと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

RのFactoRizationMachinesパッケージを用いたFMのページビューデータへの適用

Googleタグマネージャーで集めたアクセスログデータを用いて、前回と同様に記事のレコメンドにチャレンジしてみようと思います。FactoRizationMachinesパッケージという便利そうなパッケージの存在も知れたことから、今回は以前から気になっていたFactorization Machineを扱います。

【目次】
・Factorization Machine(FM)の概要
・パッケージ紹介とインストール
・サンプルデータの構造把握と前処理
・FMの実行
・結果
・参考文献

Factorization Machine(FM)の概要

  • 組み合わせ特徴量を扱う教師あり学習モデル。行列分解とSVMを合体させた手法。
  • スパースになりやすいデータの予測問題で扱う。
  • 1ユーザーのある商品に対しての評価を、1評価1行として表して、ユーザーとアイテムの交互作用の特徴ベクトルを扱う。
  • 相互作用項に関して、時間や文脈などを自由に入れられる。
  • 相互作用項を次元圧縮する際の要素数を事前に決める必要がある。
  • Matrix Factorizationよりも精度が良いとされている。特徴量エンジニアリングなどで使われているようです。(Click-Through Rate Prediction

パッケージ紹介とインストール

FactoRizationMachinesパッケージは線形SVMと2次のFMと高次のFMを実行することができ、引数で正則化項も加えることができます。現段階においては回帰のみで分類問題への適用は今後の開発となるようです。CRANから普通に
install.packages(‘FactoRizationMachines’)
でインストールします。libFMexeパッケージの場合は、libFMをインストールしてパスを指定しておく必要がありますが、このパッケージに関しては不要となります。

サンプルデータの構造把握と前処理

FactoRizationMachinesパッケージのサンプルコードにおいては、MovieLensのデータがサンプルデータとして載せられていました。ユーザーのID(整数)、映画のID(整数)、評価(整数、5段階)、日時(整数)からなるデータに対して、sparseMatrixに変換していました。

今回は、前回の投稿非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみると同じデータを使って、アクセスログデータに適用しようと思います。FactoRizationMachinesの形式に合わせるために、このブログのアクセスログも、クッキーのIDを整数に、記事のIDを整数に、閲覧回数を5段階(5以上を5に変換)に、日時を整数に変更しています。

FMの実行

デフォルトの設定c(1, 10)では線形のウェイトが有効で、2次の項の要素数が10で正則化項なしのFMを実行することになります。引数に関する詳しい情報はPackage ‘FactoRizationMachines’に書かれています。今回はサンプルを参考に正則化項ありでモデルを実行します。まず、アクセスログデータに対して、ユーザーのIDからなる整数ベクトル、記事のIDからなる整数ベクトル、セッションのあった日時のデータからなる整数ベクトルを作成し、sparseMatrix関数を用いて元データを変形し、80%のデータをトレーニングに、20%のデータをテストに割り当てます。さらに、テストデータに関して、予測値との平均二乗誤差を計算します。

結果

各モデルについての平均二乗誤差を計算しています。
線形モデルや高次元モデルよりも、2次の項を持つFMが精度が高いようです。

こちらは、この中で性能の良かった2次の項を持つFMの予測結果とテストデータの結果をプロットしたものです。4点を超える値をあまり予測できていないようです。今回はサンプルを回しただけなので、本来であれば次元の数kや正則化のセッティングをいろいろいじったり、相互作用項を新しく追加するなどして精度を高めることが必要です。

結果の比較だけでは仕事で使えないので、実際に予測した結果を取り出したいと思います。
実際に運用するとなると、ページIDを所与として、ページビュー数を0とおいて(型をそろえるため。NULLだとエラーになった)、任意のタイミング(date)を想定して、モデルにデータを適用し、評価の高いものをサジェストするスタイルになるのではないでしょうか。

この結果だと、ユーザー98に記事1を見せることに対して4.02点が与えられています。

参考文献

Factorization Machinesを今更読みました
Factorization Machines
High-order factorization machines with R #tokyor 61
Factorization Machinesのおはなし。
libFMexeを動かすまで (R Wrapper for the libFM Executable参照記事)
一歩Matrix Factorization、二歩Factorization Machines、三歩Field-aware Factorization Machines…『分解、三段突き!!』
[論文] Factorization Machines (ICDM 2010) 読んだ 22:41
Factorization machines with r
Factorization Machinesについて調べてみた

References   [ + ]

1, 2, 3, 4, 5. predict(model,data.test)-target.test)^2

非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみる

概要

2017年7月より当ブログにおいて、Google AnalyticsのカスタムディメンジョンにUser IDを取得して残せるようにしました。これにより、ユーザーの一人一人が閲覧した記事の分析を行うことが可能となっています。そこで今回は、記事の数も50近くなってきたことから、レコメンド手法の適用を行ってみたいと思います。Exampleデータでレコメンドをしても、ドメイン知識がないと、どれくらい良さそうな推薦なのかの判断が付きにくいと考えていたので、自分で運営しているサイトのデータであれば評価をつけやすいのではないかと思われます。
手法としてはNMF(non-Negative Matrix Factorization)を適用します。

【目次】
・前処理
・データ確認
・行列因子分解(Matrix Factorization)
・レコメンド結果の確認
・感想
・参考文献

前処理

ユーザー単位の閲覧した記事データを作るためにGAのAPIで取得したデータを整形し、集計します。このようなデータを手に入れるには自身のブログにおいて、GTM(Google Tag Maneger)を設置して、GTMのコンテナにカスタムJS変数などを格納する必要があります。詳しくはこちらのブログを参照ください。( 2016年の新定番!ユーザーエクスプローラーをもっと活用するための簡単な方法

このような形で行にユーザーが、列に記事が、各要素にアクセスの有無が入るようになっています。本来は評価データを用いますが、評価データがないので、アクセスしたかどうかの1-0データとして扱います。

データ確認

以下より、行列の要素のうち、値が入っている要素は全体の2.9%であることがわかります。参考文献は2.6%をスパースと言っていたことから、今回のデータもスパースであると考えることができます。

行列因子分解(Matrix Factorization)

行列因子分解はユーザーごとのアイテム購買履歴などのデータからなる行列の因子分解を行うもので、レコメンドにおいては因子分解の結果として、ゼロ以上の要素を持つ場合に推薦対象とするなどの活用がされています。非負値行列因子分解は要素が非負となるように制約を付けて行列を分解するもので、正の値を取ることから推薦に使いやすいとされています。

今回はNMFパッケージを用いて、ほぼデフォルトな設定で結果を眺めてみますが。本来であれば以下のような考慮が必要となります。最適化などをがっつり進めようとするとクロスバリデーションなどを行ったり、計算リソースが結構かかったりすると思われます。

  • アルゴリズムの選択
    • brunet,KL,lee,Frobenius,offset,nsNMF,ls-nmf,pe-nmf,siNMF,snmf/r,snmf/lがある。
    • デフォルトではbrunetとなっている。
  • シーディングメソッドの選択(初期値の決め方)
    • 固定値(none),ランダム(random),ica,nndsvdがある。
    • デフォルトではランダムとなっている。
  • 実行回数の設定
    • 100〜200回を選ぶことが多い。
  • 並列処理
    • 処理量が多いので最適化などを考えると並列処理も重要。
  • 因子分解ランクの推定
    • ランクの数が多いと残差は減るが、過学習になりやすい。クロスバリデーションを行う必要がある。

レコメンド結果の確認

各ユーザーに対して推薦されたアイテムのうち、比較的値の大きなものに関してコメントをしたいと思います。

ユーザー1
「XGBoostやパラメータチューニングの仕方に関する調査」を閲覧したユーザー1に対して、「XGBoostのパラメータチューニング実践 with Python」を推薦できている一方で、「OpenCV&Pythonで画像の類似度を計算させる〜イケメンの顔比較」など関係性の薄そうな記事も推薦してしまっています。

ユーザー2
「人工知能学会全国大会2017のWebマーケティングで参考になりそうな研究9選」を閲覧したユーザー2に対して、「Kaggleで使われた特徴量エンジニアリングとアルゴリズムまとめ」を推薦できている。

ユーザー3
「洋楽の歌詞データでDoc2vecを実行してみる」を閲覧したユーザー3に対して、「LDA(潜在的ディリクレ配分法)まとめ 手法の概要と試行まで」を推薦できている一方で、「ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる」と関係性の薄そうな記事も推薦している。本来、近いと思われる、「Word2Vecでクラシックの楽曲情報をコーパスとして類似度を出してみる」を推薦できていない。

感想

  • 結果に関して、納得のいくものといかないもので半々な印象でデフォルトの設定では満足できない。
  • NMFはクロスバリデーションなどにより過学習を避けた形での最適なランク数やアルゴリズムの選定、初期値の設定などが必要と思われるので、そこらへんの取り組みも今後行いたい。
  • FM(Factorization Machine)は回帰項を追加できるので、今後挑戦してみたい。
  • 1-0のデータではなく評価のデータが欲しい。

参考文献

岩波データサイエンス Vol.5
非負値行列因子分解(NMF)によるレコメンドのちょっとした例
An introduction to NMF package
非負ッ値行列・ファクトリゼーション ~ 半群じゃないから難しくないもん!

ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる

ベイジアンネットワークの知見が無かったので、調べた情報をまとめています。一応、載せているスクリプトでRを用いて予測するということができます。

【目次】
・ベイジアンネットワークとは
・ベイジアンネットワークの用途
・ベイジアンネットワークの推定のステップ
・Rでの実行例
・おまけ:Webサービスへの応用例
・参考文献

ベイジアンネットワークとは

・複数の確率変数の間の定性的な依存関係をグラフ構造によって表し、個々の変数の間の定量的な関係を条件付確率で表した確率モデル。前提として、有向非循環(グルグルと回らないグラフ)となっているグラフの構造を持つものに限定している。
・入力となる変数と出力となる変数はモデルの中では区別されない。
・時間という明確な因果関係などをモデルに組み込みやすいので、系列データなどを扱うケースが多い。モデル設計者がデータが生成されるプロセスを考慮しながらモデルを組んでいける。
・非循環性とd分離の仮定のみによって導かれる現在考えられる最も自然な離散モデルであり、現在の様々なモデルの中でも、最も表現力と予測力を持つモデルとされる。

ベイジアンネットワークの用途

・物体追跡
・ジェスチャ認識
・Webサイトなどのレコメンデーションサービス
・広告配信
・メルマガなどの配信最適化
など画像処理問題から消費者行動問題までいろいろな分野で活用されているようです。

ベイジアンネットワークの推定のステップ

モデルの選択や推定は学習 (learning)と呼ばれ、以下の2段階のステップを踏みます。

Step1:構造学習(structure learning)…データからネットワーク構造を学習する。
Step2:パラメータ学習(parameter learning)…Step1で学習した構造によって意味付けした分布のパラメータを学習する。

ベイジアンネットワークの推定は条件付き確率やMAP(Maximum a posteriori)推定量などを用います。構造学習の際に、事前知識を導入することができます。

Rでの実行例

・データ
何千もの主要な免疫系細胞から取り出した11種類のリン酸タンパク質とリン脂質の測定値からなるデータです。以下のサイトよりダウンロードできます。
Supporting Online Material Causal Protein-Signaling Networks Derived from Multiparameter Single-Cell Data

・パッケージ
ベイジアンネットワークを実行できるパッケージはいろいろありますが、参考文献に従いbnlearnを用います。加えて、ネットワーク図の可視化のためにRgraphvizを用います。Rgraphvizがcranでは対応したものがなかったので、こちらを参考にして、インストールしました。
Provides plotting capabilities for R graph objects

・ネットワークの図示

・データ確認

・構造学習
データセットから、ネットワーク構造を推定します。推定するための手法は様々あるようですが、今回はヒルクライムアルゴリズムを用います。

構造学習より得られたネットワーク構造を可視化します。因果関係がダメな感じのネットワークが出来ていないかチェックする必要がありますね。機械的に作るだけでなく、背景的知識も考慮するという進め方が推奨されています。

・パラメータ学習
先ほど、データセットから求めた構造について、以下のコードでパラメータ推定します。パラメータの推定や構造推定は試行錯誤するところなので、こんな簡単に済む訳では無いようです。

・予測
テストデータを用意して、データセットから求めた構造・パラメータを用いて、任意の変数の予測を行います。予測したい変数の値は欠損していたらエラーで予測してくれないですが、全部0にしておけば回ります。運用上は予測したい変数の値はわからないはずなので。

全然45度線上に乗っていないので、あまり精度は高くないようです。

cpquery関数と言って、条件付き確率を予測できる関数もあるようです。
Perform conditional probability queries

おまけ:Webサービスへの応用例

参考文献には、Webサイトの閲覧データをもとに、ユーザーの行動を予測するというモデルの紹介がなされていました。ただ、Webサイトの場合、ページ数や商品の数・ユーザーの数も膨大なことから、ネットワーク構造の推定が難しいようです。そこで、ネットワークを作るに際して、変数を確率的潜在意味解析(pLSA)で扱いやすい数に絞るなどの工夫をされていました。自社のサイトに関しても適用する際に、変数の多さは待った無しだと思うので、いざ適用する際は情報圧縮技術を駆使したいですね。

参考文献

ベイジアンネットワーク技術 ユーザ・顧客のモデル化と不確実性推論

確率的グラフィカルモデル

Learning Bayesian Networks in R an Example in Systems Biology

グラフィカルモデル入門

ベイジアンネットワークを用いたWeb レコメンデーションシステムの開発

PRML8章

Billboard100位以内の楽曲の歌詞情報にLDAを適用してみた

目次

・はじめに
・データ収集
・Rによる分析
・LDAの結果
・参考文献

はじめに

前回の投稿でBillboardの週次洋楽ランキングデータをWebスクレイピングで取得し、楽曲の消費サイクルのような順位の挙動を確かめることができました。(某洋楽ヒットチャートの週次ランキングデータをBeautiful Soupで集めてみた)今回は、歌詞の情報を用いて順位データとつなぐことにより、どのような単語の入っている洋楽がBillboardにおいてTop10に入る傾向があるのかをLDAを行うことで確かめたいと思います。

データ収集

残念なことに、Billboardのサイトに歌詞の情報は載っていません。そこで、洋楽の歌詞が取り上げられている某サイトをPython(3系)でWebスクレイピングし、名寄せを頑張って順位データと歌詞データを繋ぎます。

幸いなことに某サイトのURLに規則性があったので、アーティスト名からなるURLを生成し、そのURLをWebスクレイピングして楽曲のリストを集め、今回のBillboardのランキングに入った楽曲のみに絞ります。

楽曲をランキングに含まれるもののみに絞ったら、歌詞詳細ページを取得します。

うまいこと歌詞情報を手に入れることができました。ざっと947曲です。

Rによる分析

ここから、Rにてテキストマイニングを行いたいと思います。まず、tmパッケージを用いて、不要語(stop word)を除去します。具体的にはtheとかyouとかを除外しています。

続いて、LDAを実行できるtopicmodelsパッケージで扱えるようにするために、テキストデータに以下の処理を施します。

あとは以下のコードでLDAを実行するだけです。トピック数はアドホックに20としています。研究者の方、いい加減ですみません。

LDAの結果

まずは推定されたトピックごとの上位10単語をみてみます。トピック1はラブソングとかでしょうか。トピック17にパリピっぽい単語が、トピック18にスラングが含まれていますね。

見ずらいので、行を一つにまとめて、トピックにidを割り振ります。

最後に、BillboardでTop10に入ったかどうかのデータを作っておき、そのデータと各歌詞を繋ぎ、各歌詞ごとに割りふられた確率が最大のトピックで集計をします。

BillboardのTop10ランクイン割合の高いトピックTop3
「one,ooh,call,cause,gettin,born,day,makes,came,stand」
「better,world,whoa,run,light,things,find,show,see,waiting」・・・明るい感じ?
「stop,just,hands,put,party,crazy,live,lights,play,see」・・・パリピぽい

BillboardのTop10ランクイン割合の低いトピックTop3
「wanna,want,take,rock,see,kiss,come,make,body,tonight」・・・欲求系?
「feel,heart,life,away,just,break,real,enough,every,find」・・・癒し系?
「hey,said,old,every,woo,left,told,nothing,daddy,sweet」

あまり洋楽を聴かないので、得られたトピックの解釈が中々できないのがもどかしいです。ただ、スラングの歌詞を含む歌詞はそんなにランクイン割合が悪いわけではなさそうですね。洋楽をもっと聴いて、前処理などもう少し工夫してリベンジしたいですね。

参考文献

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

Pythonクローリング&スクレイピング -データ収集・解析のための実践開発ガイド-

モダンなRによるテキスト解析topicmodels: An R Package for Fitting Topic Models

某洋楽ヒットチャートの週次ランキングデータをBeautiful Soupで集めてみた

はじめに

知人より、洋楽の流行りに疎いのでキャッチアップしたいという要望があり、某洋楽ヒットチャートの週次ランキングとTop100のデータを大量に集めてみようと思うに至りました。今回は深い考察を行うには至っていませんが、簡単にRにて集計・可視化を行います。

データ収集

Webスクレイピング対象の某洋楽ヒットチャートの週次ランキングは今週の順位・先週の順位・アーティスト名・曲名・詳細ページへのリンクなどが載せられおり、毎週土曜日更新されています。サイト内から導線はありませんが、URLのパラメータに法則があるため、うまく収集できます。今回は2010年8月〜2017年6月の約7年分のデータを集めます。

URLのリストをCSVで読み込み、BeautifulSoupでタグを指定して抽出します。

データ取得後は簡単にpandasのstr.replaceで整形すると、以下のような結果になります。今週の順位と先週の順位が引っ付いてしまっています。

ここから横着してRで整形し、各週の順位データを作成しました。

データ確認

データ構造はこのような形です。

まず、どんな楽曲やアーティストがランキングに入っているのかを簡単に確認してみます。 ほとんど聞いたことない人の名前、曲名ですが。。

続いて、2010年8月〜2017年6月の間に100位以内に入った数を楽曲ごとにヒストグラムにしてみます。べき乗分布な形かと思いきや、20回前後で盛り上がっているのが気になりますね。

中央値が9週間なので、意外と長い期間Top100には入っています。

続いて、100位以内に入った数をアーティストごとにヒストグラムにしてみます。こちらはべき乗分布のような形になっています。

Top10入りの楽曲の実態

Top10に入っている楽曲のみに絞って、ヒストグラムを描いてみます。

Top10に入ったら、10週近くは10位以内に含まれるようです。上位はすぐに取って代わられるのかと思いきや、人気が人気を呼ぶとかなのでしょうか。確か、某ECサイトの方が、生キャラメルは売れるから売れたんだとか言っていた気がします。

順位の推移

100位以内にランクインした回数が最も多かった楽曲のTop10に関して、時系列プロットをしてみます。


初回に上位にランクインして、後は下がるだけの楽曲や、じわじわとランキングを上げていく楽曲などが観察されています。楽曲の消費のサイクルみたいなものがあるのでしょう。

今後について

せっかく面白そうなデータが手に入ったので、リンク先も辿って、どのような楽曲やアーティストの特徴が人気に繋がりうるのか見てみるのも良いですね。あと、もう少し洋楽聴いてみようと思います。私はクラシック音楽とジャズしかウォークマンに入っていないので。

参考文献

Pythonクローリング&スクレイピング -データ収集・解析のための実践開発ガイド-

ラーメン二郎の某飲食店レビューサイトデータに対して共分散構造分析をしてみる

データ

ラーメン二郎に関して、某飲食店レビューサイトのデータをWebスクレイピングしたもので、「料理・味」・「サービス」・「雰囲気」・「CP」・「酒・ドリンク」の評価項目に関して、1~5の実数値が割り振られています。ラーメン二郎の店舗数40店のうち、欠損のなかった37店舗の評価データとなります。

記述統計量・ヒストグラム

まずは記述統計量とヒストグラムを見てみます。

まぁ、お酒を飲むところではないので、drinkは低いですよね。しかしながら、総じて3点台なのは驚きです。昔行ったことがあって、雰囲気は決して良くはないはずなので。

共分散構造分析

今回、検証したい仮説は「ラーメン店としての質が二郎愛につながるかどうか」です。

ラーメン店としての質に繋がりそうな評価項目
「料理・味」
「CP」
「酒・ドリンク」

二郎愛に繋がりそうな評価項目
「料理・味」
「サービス」
「雰囲気」

加えて、料理・味とCPは関係していそうな項目なので、その点もパスにおいて考慮しておきます。

以下の図は仮説のイメージです。

パスを描写する際の細かい指定は参考文献を参照されると良いと思います。
さっそく、モデルを作って実行します。

こちらが推定結果です。

lavaan (0.5-23.1097) converged normally after 100 iterationsとあるので、適切に推定されたようです。自由度が2とデータ数が少ないためかなりギリギリな推定となっています。ただ、推定すべき母数の数よりもデータ数が一応多い状態ではあります。二郎インスパイア系の店のデータも集めた方がいいかもしれないですね。

モデルの評価

モデルの評価として適合度と母数の推定に関して見ていきます。

適合度
・適合度指標であるCFI(Comparative Fit Index)が1なので、適合度に関しては良さそうです。
・同じく適合度指標であるTLI(Tucker-Lewis Index)が0.998なので、適合度に関しては良さそうです。
・0.05以下であれば当てはまりが良いとされるRMSEAは0.028なので、当てはまりは良さそうです。
・0に近いほど良いとされるSRMRは0.024となっています。

母数の推定

f2 ~ f1はラーメン店としての質と二郎愛の関係を想定したものですが、p値が0.47と全然だめでした。ラーメン店としての質が二郎愛に繋がるという仮説は正しいとは言えないです。

参考文献にもあるよう、係数の解釈を行いやすくするために標準化推定値を求めます。

これによると、f1(ラーメン店としての質)に対してはコストパフォーマンスが最も影響を与えるようです。味よりもコストパフォーマンスが勝っているという考察になりますが、それはそれで面白いですね。他方、f2(二郎愛)に対しては味よりも雰囲気・サービスが影響を与えるようです。

semPlotパッケージを用いてパス図を出力してみます。

変数名が3文字に省略されているのですが、修正方法がパッと見つからなかったので、そのまま載せています。

初歩の初歩ですが、一通りの進め方がわかったので、今後も共分散構造分析にチャレンジしてみたいと思います。

参考文献

共分散構造分析 R編―構造方程式モデリング

RStanで学部時代の研究を振り返ってみる

研究概要

大学時代に実験経済学で行った実験結果のデータがUSBに入っていたので、振り返って分析などをしてみたいと思います。

研究目的
 ピア効果に関して、競争相手が自分よりも秀でた人がいいのか劣った人がいいのかを確かめる。

実験方法
・1分間で100マス計算を2セット解いてもらう。(めちゃ速い人には3枚目も渡した)
・実験開始後、実験対象のクラスによって、途中で「平均的なクラスは○○マスまで進んでいます!」とアナウンスします。アナウンスすることで、競争相手のレベルを知り、焦るなり余裕を感じるなりしてもらおうという計画です。
なお、対照群はアナウンスをしていません。アナウンス内容は「平均告知(18秒)、上告知(15秒)、超上告知(12秒)、下告知(20秒)」と4パターンとなります。
・計算が間違っているものは加点しません。

実験対象
 某国立大学の経済学部生の1~2年の必修科目履修者217名(先生に交渉して授業の開始5分を頂いて実験を行いました。)
 内部進学やスポ専などがない分、計算能力的にある程度近い集団ではないかと思われます。

検証方法
 アナウンスごとに100マス計算の点数の水準が変わりうるのかを回帰分析などで判断。

データ可視化

以下、実験カテゴリごとの略記です。
下告知(20秒)・・・slow_20
上告知(15秒)・・・fast_15
超上告知(12秒)・・・fastest_12
平均告知(18秒)・・・average_18
対照群・・・baseline

データ構造の確認です。

平均値、中央値、標準偏差、サンプルサイズを出してみます。

中央値で見てみると、baselineに対してわずかですが点数に違いがありそうに見えます。

実験種別で点数に関するヒストグラムと確率密度関数を確認してみます。


baselineが多峰性がありそうなのが気になります。average_18は低そうに見えますね。

RStanで重回帰

『StanとRでベイズ統計モデリング』にあるコードを参考にしています。正規分布を事前分布にした線形回帰モデルです。
被説明変数が点数、説明変数が実験種別のダミー変数だけからなります。

結果

traceplot(fit)でMCMCのサンプリング結果を確認する。

収束しているように見えます。

以下推定結果ですが、残念ながらベイズ予測区間において符号の逆転が起きていないものはなかったので、アナウンスによる効果があるとは言えないようです。ただ、slow_20の係数がおしいですね。少なくとも他の実験種別よりも、アナウンス効果があるかもしれないという考察に止まりそうです。

分布でも見てみます。アナウンス効果が0を確実に超えているとは言えないですね。

結局、学部時代のレポートと結論は変わらないのですが、係数が0よりも大きい確率という観点で結果に向き合えたのは良かったと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)

PythonやRを用いたアルゴリズム取引・株価分析まとめ

はじめに

まわりでシステムトレードや株価の機械学習による予測などに関心が高まってきたので、私も少し調べてみようと思いPythonやRで行われた分析・実装の事例を集めてみました。自分の資産を突っ込む気にはなれないですが、事例を知っておくだけ知っておきたいですね。

調査法

Google検索において以下のクエリで上位に表示されたサイトを中心にまとめました。
「python 機械学習 株価」
「python 機械学習 為替」
「python アルゴリズム取引」
「python machine learning stock price」
(同様にRも調べました。)

アルゴリズムトレードの理論から学ぶ

システムトレードで億万長者になるぞ! coursera で Computational Investing Part I を受けた
こちらはCourseraで開かれていたシステムトレード向けの講座を受けた方のブログです。Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Riskはもっとも重要なアルゴリズムトレードの本とされているようです。出版年が1999年と結構古いですが。
Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Risk (McGraw-Hill Library of Investment and Finance)

Rを用いた株価データのモニタリング

Stock Analysis using R
こちらはアルゴリズム取引とかではないですが、Rで株価を簡単に取得できるquantmodライブラリについて紹介されています。

以上のコードを実行しただけで時系列データがボックスプロットで描写されました。非常に便利そうです。複数のグラフを並べてディスプレーで見てみたいですね。図ではテクニカル分析で使われるボリンジャーバンドを一発で出してくれています。傾向を掴むための分析に役立つのではないでしょうか。

MSFT_stockPrice

RSIで株価の連動性を見る

過去のデータからビッグデータ分析で株価を予測する
ここでは、主にpandasを用いて株価の分析を行っています。RSI(Relative Strength Index)という「一定の期間変動幅の中でどれ位株価が上昇しているのか、下落しているのかをはかるもの」を計算して日経平均との相関を見ています。データさえあれば、個別株の市場連動性を見るぶんにはpandasで十分に分析できそうです。

アルゴリズムトレードのシステム構成

自動トレードボット
こちらでは、仮想通貨取引のための自動トレードボットの作成のための手順などが書かれています。
工程としては
・仮想通貨価格データ取得
・バックテストの実施
・明日の価格予想
・学習パラメータの最適化
・結果のメール送信
・APIを使った成行トレード
・ジョブのスケジューリング
などとなっているようです。
完全自動化しようと思うと、作るのは大変そうですね。

WEB屋の自分が機械学習株価予想プログラムを開発した結果
こちらはよりシンプルで
「チャートを見る限り、誰がどうみても今日値上がりする銘柄」を検索して、毎朝Slackに通知してくれるシステム
をPythonで作成されています。最終的に資産を溶かす形になっているようです。

アルゴリズムトレード向けのPythonライブラリ

pythonのアルゴリズムトレードライブラリ
こちらでは
zipline( http://www.zipline.io/
PyAlgoTrade( http://gbeced.github.io/pyalgotrade/ )
pybacktest( https://github.com/ematvey/pybacktest )
backtrader( https://github.com/mementum/backtrader )
というライブラリが紹介されています。

遺伝的アルゴリズムで為替の自動売買

pythonと遺伝的アルゴリズムで作るFX自動売買システム その1
遺伝的アルゴリズムでFX自動売買 その2 進化する売買AIの実装
遺伝的アルゴリズムでFX自動売買 その3 OandaAPIで実際に取引
こちらは、外国為替の取引をPythonで自動化させた試みです。最終的に負けてしまってはいますがシステム構成についても詳しく書かれているので、勉強になりました。データはOANDAというサービスが提供しているAPIを用いて取得し、GA(遺伝的アルゴリズム)を使って為替の売買タイミングを決めているようです。実際の売買にもOandaAPIというものを利用して完全に自動化させています。

決定木による株価リターンの予測

機械学習で未来を予測する – scikit-learn の決定木で未来の株価を予測
こちらはPythonでscikit-learnを用いて、決定木による株価の予測をされています。目的変数としてはリターンインデックスをおいています。前処理においてpandasが使われています。

ランダムフォレストによる株価の予測

Pythonで機械学習を使った株価予測のコードを書こう
こちらでは、ランダムフォレストを用いた機械学習で、ETFなどのデータを特徴量にして個別の株価予測を行っています。

アカデミックサイドでの研究事例

Stock Price Forecasting Using Information from Yahoo Finance and Google Trend
UC Berkeleyの経済学部生の研究です。こちらはR言語で、ヤフーファイナンスとGoogleトレンドの情報を用いて株価を予測している研究です。これまでの時系列手法よりもパフォーマンスが良いそうです。

最後に

事例を集めてみて、どこまでのレベルのものに手を出すべきか悩ましいと思いました。完全自動化は維持するコストがかかりそうですし、ロスが増大しないか心配です。自分としては情報を自動で取得し、リターンが発生する確率の高そうな銘柄をサジェストしてくれるレベルで十分な気がします。

追記

人工知能学会の全国大会2017でAIを用いた株式運用に関する研究がなされていました。深層学習を用いたものが多そうです。
株価変動パターンの類似性を用いた株価予測
深層学習と高頻度データを用いた株式注文状況の推定
LSTM を用いた株価変動予測
深層学習を用いた株価予測の分析