警察庁オープンデータの前処理と死亡事故発生予測のための機械学習について

はじめに

先日、近所の幹線道路の信号付き横断歩道を横断しようとした際に、2メートルほど先で回避不能な速度で車が赤信号を無視して横断歩道を突っ切きるという非常に恐ろしいイベントに遭遇しました。
そこで、私は「このイベント、何回かあったら死亡事故が起きるぞ」と思いました。警察署に近隣道路の不満について電話をしてみたのですが、「警察もちゃんと考えて取り締まっていたり対策してますので」と言われ、何の成果も獲れない電話となりました。
そのようなモヤモヤの中、私は民間でできることは「客観的なデータでもって、エリア、曜日、天候、時間帯ごとの、死亡事故のリスクを見積もり、それを公開して注意喚起していくこと」なのではないかと思うに至りました。
よって、警察庁が公開している交通事故のオープンデータについて、その前処理と機械学習によるエリアの死亡事故の発生予測などについて記していきます。

前処理のためのデータやコードはこちらのGitHubにあります。
https://github.com/KamonohashiPerry/traffic_accident_analysis
スターを10個以上いただけて、データ作って良かったなと思いました。

目次

・データについて
・先行研究について
・前処理について
・EDAからわかること
・機械学習
・さいごに

データについて

2019年と2020年の交通事故に関するデータが警察庁よりオープンデータとして公開されています。
オープンデータ 2020年(令和2年)
オープンデータ 2019年(平成31年/令和元年)

どのようなデータかというと、
それぞれ2019年と2020年の間に発生した全国の事故について、発生した状況(時、場所、天候)やその結果に関するデータからなり、項目の数としては58個とオープンデータとしてはかなり豊富な情報量となっています。

データは本票、補充票、高速票の3つからなるのですが、今回は本票だけを用います。
高速道路の精緻な分析をしようと思うと、高速票を繋ぐ必要がありますが、私の目的は一般道にあるので、繋がないものとします。

先行研究について

事故の予測に関する先行研究

交通事故の分析に関する研究は以下のものがあり、参考にしました。
死亡事故を予測するという問題設定として、その特徴量としては場所、時間帯、曜日、場所かつ時間、場所かつ曜日かつ時間、天候、季節、風速・風向き、車線の数、歩道の有無、制限速度などを用いていることがわかりました。
そして、データをクラスタリングするためにDBSCANを用いたり、不均衡データであることから、アンダーサンプリングを行ったりしているようです。
最後に、アルゴリズムはSVMやロジスティック回帰、深層学習と色々と使われているようです。

・Live Prediction of Traffic Accident Risks Using Machine Learning and Google Maps

・Data-Driven Urban TrafficAccidentAnalysis and Prediction Using
Logit and Machine Learning-Based Pattern Recognition Models

・A Deep Learning Approach for Detecting Traffic
Accidents from Social Media Data

・A BASIC STUDY ON TRAFFIC ACCIDENT DATA ANALYSIS USING
SUPPORT VECTOR MACHINE

・Predicting Crash Injury Severity with Machine Learning Algorithm Synergized with Clustering Technique: A Promising Protocol

・Comparison Analysis of Tree Based and Ensembled Regression Algorithms for Traffic
Accident Severity Prediction

その他

・Taxi Demand Prediction System

・2019年04月号トピックス 警察庁のオープンデータを活用した独自AIによる犯罪発生予測

前処理について

ソースコードに関しては、GitHubに上げておりますので、
.pyなら
https://github.com/KamonohashiPerry/traffic_accident_analysis/blob/main/src/preprocess.py
https://github.com/KamonohashiPerry/traffic_accident_analysis/blob/main/src/preprocess_second.py

.ipynbなら
https://github.com/KamonohashiPerry/traffic_accident_analysis/blob/main/Traffic_Accident_Analysis.ipynb
を参照ください。

基本方針

・天候、時間帯、曜日、平日、道路の形状、スピード制限など、事故の結果(死亡の有無)を明らかに予測できそうなもの以外を特徴量にします。
・メッシュ単位での過去の事故に関する情報に関する特徴量の作成(学習を行う時点よりも過去の時点のデータのみ利用)

詳細

・事故の発生日時のdatetime化(列で年・月・日・時・分が分かれているため)
・マスタ情報(PDF)をもとにマスタのIDを名前に変換
・緯度経度のデータについて、度分秒から度への変換を行い、3次メッシュ(1km四方)のデータを作成
・メッシュ単位での特徴量として、過去の事故発生件数、死亡数などを集計
・メッシュ単位での特徴量として、過去に発生した事故に関するもろもろのデータの集計(そのメッシュで事故を起こした平均年齢や事故が発生してきた時間帯など)
・もろもろの変数のダミー変数化

EDAからわかること

EDAに関しては.ipynbのみで行っております。
https://github.com/KamonohashiPerry/traffic_accident_analysis/blob/main/Traffic_Accident_Analysis.ipynb

私がよく分析で使うのは、データの該当件数とカテゴリごとの率の可視化なのですが、以下の図の見方としては、青の棒グラフが各カテゴリの件数、赤の折れ線グラフが死亡率となっています。

都道府県


福井県、鳥取県、秋田県、島根県、岩手県などは死亡率が高そうです。

道路種別


農道(815件)や林道(88件)で起きる事故非常に死亡率が高くなるようです。

時間帯


夜間や明け方は死亡率が高くなるようです。

天候


霧の天候(294件)では死亡率が高くなるようです。

人口密集地かどうか


市街地ではないところでは死亡率が高くなるようです。

道路の形状


踏切-第一種(254件)はかなり死亡率が高くなるようです。カーブも死亡率が高くなります。そういえば実家の近くにも魔のカーブと呼ばれるところはありますし、事故も起きていました。

信号機のタイプ


点滅-3灯式(2566件)や点灯-押ボタン式(7703件)は死亡率が高いようです。

道路の幅


大きな交差点では死亡率が高い傾向にあります。単路に関しては一様な傾向ではないですが、小さい方が死亡率が高そうです。

道路の区切り


道路の区切りがペイントだけのケースは、中央分離帯がある場合よりも死亡率が高い傾向にあります。あと、中央線がポストコーン(チンアナゴみたいなやつ)だと死亡率が上がるようです。

歩道との区別


道路について歩道との区分がない場合に死亡率が高くなるようです。

速度制限


速度域が高くなるにつれ、死亡率が高くなるようです。あと、速度規制がされていないところも高くなっています。

曜日


土日の死亡率が高い。

祝日・その前日


祝日とその前日の死亡率は高い傾向にある。

以上のことから、得られた考察としては、
・土日や祝前日の道路はいつもより危ない
・夜間や明け方の道路は危ない
・制限速度の高い道路は危ない
・大きな交差点は危ない
・点滅信号は危ない
・カーブは危ない
・霧は危ない
でした。普段、出歩く際に意識しておくといいと思います。

機械学習

教師データを事故発生からの死亡有無(1or0)、特徴量を先行研究を参考にオープンデータから作れるものとして、
基本的には、LightGBMでの学習、ハイパーパラメータ探索、モデルの訓練(メッシュ番号でグループ化したクロスバリデーション)、テストデータ(学習は2020年11月まで、テストは2020年12月以降としている)での予測、日本地図への可視化を行っています。

ソースコードはGitHubに上げてあるので、ここでは記しません。

.pyなら(学習のみ)
https://github.com/KamonohashiPerry/traffic_accident_analysis/blob/main/src/train_model.py

.ipynbなら(学習だけでなく、予測と日本地図への可視化もしている)
https://github.com/KamonohashiPerry/traffic_accident_analysis/blob/main/Traffic_Accident_Analysis.ipynb

10-foldのCVをメッシュ番号でグループ化して学習させた際のAUCです。高い部類に入りますが、今回は不均衡データなので、抜群に当たっている感じではないです。


事故の死亡率は、0.9%なので、この予測の場合、2.2%が死亡すると当てれているということになります。また、予測から、全体の72%の死亡事故を当てているという結果です。
この性能に関しては今後、もっと高めていく余地はいくらでもあると思います。

こちらは、テストデータでの予測スコアの階級値(青い棒グラフ)と実際の死亡率(赤い線グラフ)を示したものです。スコアが高いものは信じてもいいかなぁと思いました。

特徴量重要度の上位10件は以下の通りです。

時間帯とか、防護策とか、信号機のタイプとか、市街地かどうかなどが予測に役立っているようです。

死亡事故発生を予測した確率を日本地図の関東エリアにプロットしたものが以下の図になります。

市街地のほうが死亡事故が起きやすいように見えますね。

東京都だけに絞ってみましたが、危険なエリアは都心に点在はしていますが、郊外のほうが死亡率が高そうです。

さいごに

交通事故にあいかけて、その怒りの感情を警察庁のデータにぶつけてみましたが、死亡事故が起きやすい条件などがある程度見えてきたと思いました。もちろん、なんとなく思っていた傾向もあったのですが、数字でハッキリと出せるのは良いですね。今後、旅行などをする際に、確率の高いエリアを警戒して行動するようにしたいです。
そして、地図への可視化は普段の仕事でほとんど行わないのですが、こうしてデータに触れて可視化の練習ができたのも有意義に思いました。
この警察庁のデータは分析するまでに工数がかかるデータとも言えるので、この記事をとっかかりとして、オープンデータの分析に挑戦して安全な生活を送る上での知見をシェアしていけると良いな、というか民間としてはこれくらいしかできないなと思っています。

参考情報

・[第3版]Python機械学習プログラミング 達人データサイエンティストによる理論と実践 (impress top gear)
・Pythonでプロットした地図を出力する
・Pythonで度と度分秒を相互に変換する

Bayesian Methods for Media Mix Modeling with Carryover and Shape Effects[A4一枚まで備忘録]

A4用紙1枚にまとめるイメージでメモを残そうという取り組みです。

今回はメディアミックスモデリングに関して仕事があった際に参照できるメモの一つとして残します。
先日、Python/STAN Implementation of Multiplicative Marketing Mix Modelに参考文献としてBayesian Methods for Media Mix Modeling with Carryover and Shape Effectsの存在を知りました。時系列の重回帰の一種じゃないかと思っていたんですが、天下のGoogle様がどのようにマーケティングの成果を測ろうとしているのか、この論文から知ってみようと思います。

ちなみに、Python/STAN Implementation of Multiplicative Marketing Mix Modelの写経したものはMarketingMixModelingに残しています。

このコードでマーケティングにおけるチャネルごとのROAS(Return On Advertising Spend)を推定することができます。

以下の図はコードを回すことで得られるものですが、チャネルごとのROASの事後分布となります。

この図もコードを回すことで得られるものですが、広告などの支出に関する売上の関数をフィットしたものになります。精度高くこれらの関数を求めれば、過剰に広告に投資している可能性のあるチャネルを見つけることが可能になります。

目次

Abstract
Introduction
Model Specification
Estimating the Bayesian Model
Attribution Metrics and Optimal Media Mix
Application to a Simulated Data Set
Impact of Sample Size on Estimation Accuracy
Choice of Priors
Application to Real Data and Model Selection
Conclusion

Abstract

ROASやmROASに関するアトリビューション分析をベイズモデリングで行い、事後分布から効果を推定している。
広告関連のデータが、売上などに対して、数ヶ月後に遅れて効果が現れるキャリーオーバー効果の想定、支出を増やしても効果がなくなる収獲低減の想定などをモデルに置いている。
モデルの識別に関してはBICを評価基準として用いている。

Introduction

  • Media mix models(MMM)はメディア支出が売上にどのように影響を与えるのかを理解するために、最適なメディア投資を行うための支出の配分を決めるために使われる。
  • 扱うデータとしては売上、物価、プロダクトの分布、様々なメディアの支出、マクロ経済や天候、季節性、競合の情報などの外部要因が含まれる。
  • RCTが難しいため、あるいは傾向スコアなどのモデリングもデータが少なく難しいため、回帰分析が行われることが多い。
  • 線形回帰だと、広告がサチっている効果(広告が飽和している状態)や収獲低減の効果も表現できない。
  • 広告にはキャリーオーバー効果と言って、売上への効果が遅れて現れるという可能性が広く信じられている。
  • 事前分布に広告に関する様々なノウハウを反映させることで、データ数に比して推定するパラメータの数が多い分析に向きあう。

Model Specification

  • モデル構築には週次データが使われるケースが多い。
  • 売上のデータ、メディア支出のデータ、コントロール変数として売上と関係してそうなデータが扱われる。コントロール変数は業界によって異なる。
  • キャリーオーバー効果
    $$ adstock(x_{t-L+1, m}, \dots, x_{t, m} ; w_{m} ,L) = \frac{\sum_{t=0}^{L-1} w_m(l) x_{t-l, m} }{\sum_{l=0}^{L-1} w_m(l)} $$

    \( w_m(l) \)は非負のウェイトでラグ期間の関数になっている。
    \( L \)はキャリーオーバー効果の最大期間。
    \( x_{t, m} \)はある期間のあるメディアへの支出。
    論文ではL=13とされている。
    ウェイトの関数について、メディアの売上に対する効果は徐々に減っていくが、効果が遅れて生じることもあるので、効果のピークをずらした表現も扱われている。
    ウェイトの関数で最初にピークがくる場合、

    $$ w_{m}^{g} (l; \alpha_{m}) = \alpha_{m}^{l}, \\ l = 0, \dots, L-1 , \\ 0 < \alpha_m < 1 $$

    ウェイトの関数で遅れてピークがくる場合、

    $$ w_{m}^{d} (l; \alpha_{m}, \theta_m) = \alpha_{m}^{(l – \theta)^2}, \\ l = 0, \dots, L-1 , \\ 0 < \alpha_m < 1 ,\\ 0 \leq \theta_m \leq L -1 $$

    となる。
    \( \alpha_m \)はメディアの効果の、ある期から次の期への維持率。
    \( \theta_m \)は遅れてピークが現れる程度に関するパラメータ。

  • 形状効果
    $$Hill(x_{t,m};\mathcal{K}_m, \mathcal{S}_m) = \frac{1}{1+ \left (
    \frac{x_{t, m}}{\mathcal{K}_m} \right )^{- \mathcal{S}_m} }, \\ x_{t, m} \geq 0 $$

    \( x_{t, m} \)はある期間のあるメディアへの支出。
    \( \mathcal{S}_m \)は形状に関するパラメータで傾きに関わるもの。
    \( \mathcal{K}_m \)は形状に関するパラメータでサチるポイントに関わるもの。
    分子に\( \mathcal{K}_m ^{\mathcal{S}_m} \)を足して引くことで、

    $$ \beta_{m} Hill_{m} (x_{t, m}) = \beta_{m} – \frac{\mathcal{K}_m^{\mathcal{S}_m}\beta_{m}}{x_{t,m}^{\mathcal{S}_m}+\mathcal{K}_m^{\mathcal{S}_m}} $$

    に変換できる。
    なお、\( \beta_m \)は各メディアの回帰係数となる。これでもって各メディア支出の売上に対する効果を見れる。パラメータの大きさによって様々な形状をとりうる。

  • キャリーオーバー効果と形状効果をモデルで一緒に表現する。
    シンプルにするためにメディア間のシナジー効果を無視している。

    $$ y_t = \tau + \sum_{m=1}^{M} \beta_m Hill(x_{t, m }^{*}; \mathcal{K}_m,\mathcal{S}_m) + \sum_{c=1}^{C} \gamma_{c}z_{t, c} + \epsilon_t $$

    \( y_t \)は売上。
    \( x_{t, m }^{*} = adstock(x_{t-L+1, m}, \dots, x_{t, m} ; w_{m}, L) \)は前述の広告支出のキャリーオーバー効果を表現したもの。
    \( \tau \)は広告支出によらないベースラインの売上。
    \( z_{t,c} \)はコントロール変数。
    \( \gamma_c \)はコントロール変数の係数。
    \( \epsilon_t \)はホワイトノイズ。

Estimating the Bayesian Model

  • 階層ベイズに拡張することもできる。
  • ハイパーパラメータはメディアの効果が非負であるとか、割合は0から1だとか色々と事前情報として設定する。

Attribution Metrics and Optimal Media Mix

  • 経路ごとのROSAやmROASを計算してメディアの最適化をしたい。
  • メディアに投資することでの予測売上と、投資なかった場合の売上の差分でもって成果を測る。
  • ROASの定義式
    $$ ROAS_{m} = \frac{\sum_{t_{0} \leq t \leq t_{1} + L – 1 } \hat Y_{t}^{m} (x_{t-L+1, m}, \dots, x_{t, m}; \Phi) – \hat Y_{t}^{m} (\tilde{x}_{t-L+1, m}, \dots, \tilde{x}_{t, m}; \Phi) }{\sum_{t_{0} \leq t \leq t_{1}} x_{t, m} } $$

    \( \hat{Y} \)は売上の予測。
    \( \tilde{x} \)はメディアの支出の変化。

  • mROASの定義式
    $$ mROAS_{m} = \frac{\sum_{t_{0} \leq t \leq t_{1} + L – 1 } \hat Y_{t}^{m} (\tilde{\tilde{x}}_{t-L+1, m}, \dots, \tilde{\tilde{x}}_{t, m}; \Phi) – \hat Y_{t}^{m} (x_{t-L+1, m}, \dots, x_{t, m}; \Phi) }{ 0.01 \times \sum_{t_{0} \leq t \leq t_{1}} x_{t, m} } $$

    \( \hat{Y} \)は売上の予測。
    \( \tilde{\tilde{x}} \)はメディアの支出の1%の変化。

  • ROASの分布は事後分布と上述の定義式から得られる。
  • メディアミックスの最適化
    予算という制約付きの極値問題を解く。ラグランジュ未定乗数法を用いる。

Application to a Simulated Data Set

  • シミュレーションしたデータを使っている。2年間の週次のデータで、売上、3つのメディア、一つのコントロール変数からなる。

Impact of Sample Size on Estimation Accuracy

  • 単一のデータセットのみで評価することが難しいので、それぞれが独立した500個のデータセットを生成し、それぞれ推定をしている。

Choice of Priors

  • サンプルサイズが小さい場合、事後分布が事前分布の影響を大きく受けるため、事前分布に関して色々と考察している。
  • \( \beta \)に関する事前分布。設定次第で結構変わる。メディアの支出にかかってくるので分析に与える影響は大きい。
  • \( \mathcal{K}_m \)に関する事前分布。要はサチるポイントに関するもの。設定次第でグラフの挙動が結構変わる。

Application to Real Data and Model Selection

  • シャンプーの広告主の実際のデータを使って分析。
  • 2.5年の週次データで、反応変数として売上、メディアデータとしてTV・雑誌・ディスプレイ・YouTube・検索、コントロール変数として小売のデータ(平均オンスあたり価格)・重み付けされた全商品の流通量や広告量が扱われている。
  • BICを頼りにモデルを選択した。BICは事後分布のデータから計算できる。

Conclusion

  • メディアミックスモデルはモデルの識別によって結果が大きく変わる。
  • ベイズモデルでキャリーオーバー効果や形状効果を考慮したモデルを求め、それらの事後分布からROASやmROASを計算した。
  • データ数が少ないとバイアスの伴った推定となる。事前情報の選定にも慎重になるべき。
  • BICでモデルを識別した結果、キャリーオーバー効果や形状効果を想定した設定が選ばれた。
  • 残差に自己相関があったため、今回のモデルはまだ不完全で改善の余地がある。

[Python]Music×AnalyticsというイベントでLTをやってきましたレポート

今回はゆるめで、あるイベントのLTに参加してきた感想やら、資料やら動画やらソースコードを共有するだけの記事となっております。

ことの成り行き

Tokyo.Rでお世話になっているy__mattuさんから、「Music×Analytics Meetup Vol.4」のイベントのことを教えていただき、せっかくなのでLTをさせてもらえることになりました。

音楽の分析に関して入門レベルなので果たして参加者が盛り上がってくれる何かを提供できるか不安でした。

納期的にも、今の実力的にも過去に書いたブログ([Python]スペクトグラムで楽器の試奏データを比較(プロ奏者とも))の延長線上のことでいいなと思ったんですが、
先行研究は欲しいし、分析対象の音楽もスケール(音階)ではなく何かの曲が良いと思いました。あと、意外な情報も載せたい、分析結果を活かして便益を得たいとも思いました。

なので、LTまでに先行研究を見つけ、意外な情報としてjupyter notebookで録音する機能を見つけ、分析対象の曲にギャツビーのCMのものを採用しました。
そして、この研究から自分の楽器の演奏の腕前を向上させることができました。
これならLTとして発表しても入門レベルでも楽しいのではないかなと思いました。

LT内容

Google Presentationと事前に録音したYoutubeをこちらに載せておきます。

発表音声の文字起こし

LTを7分以内に抑えるリハーサルをした際の音声があったので、それをGoogle Speech To Text APIを使って文字起こしをしておきます。
「それではその音に近づくための研究と練習ということで江口さん自己紹介ですが名前は救えバックグランドは経済学や統計学仕事では同席部のマネージャー兼データ分析さんはこんな感じで趣味はトランペットブログ作成をやってるので良かったみてください結局どうしたんですけどその他データ分析関連のかって言いよったでどうやって勝手の音消したいなって思いますがここまでのホスピスがまずあって写真撮って進度が速い遅いである時と言うとって治ってるということですが音の長さが変わりますなのでお体変わります秘密のピストンだけなので発送の仕方ありませんのでコントロールして40種類ぐらいのレベルからは死んだことに帰って前記処理学会のそれの大会でスペクトルグラムを使ってプロのアクセスをしたのとの比較ではその方向の出来不出来を判定すると犬が行われていました鬼時間に周波数色の濃淡で信号ってなっちゃったから明日グラフとなりますがこの場合は時間が12秒まであって周波数が8000まで行ってあのね精神旺盛な強さが違ってます本当の鳴き声何もとしてるの俺も行ってここだったから体操マット耐熱耐酸プルマンちょっと寄りますよねプレイヤーの力も貸したして自分との違いを明らかにして違いがある場合とならそこによるそれが産めるのかお話していきたいと思ってます結果こんな感じえーとまずはサイトメトロ ffmpeg をインストールしてああいパイソン mrtg をインストールしたギターの音は愛想の良いしらすを呼び出してこういう形で書くことでノートブックので他のところでエロ金ボタンが現れますこれ音声は録音開始でもう一度録音するようになっていますその後に録音したファイルを保存 SMS でまあ16000hz 2のファイルに変換してまた再読み込みという形があってもはい録音環境としてはわからないぞマイクでして何楽器はちょっと古いので今回あったのがという音楽などの分析見ております玄関のために必要なセッティングってのがあるんでは基本的なセッティングで言っても実際に自分のそれかしたら私に出てきますね自分のやってみようと思いますまず聞いてみましょう第4フォルマント信号の強さがえっと Pro は赤いんですけど強いんですけど出してもらうわいえ私なってます遅く音源第4フォルマントの方がやっぱり行けないなっていうのが分かれば鼻が詰まってるから行けないっていましたねと自分ではまあ剛てる所間違いが分かりましたどうしたら強い信号になってるのかその女お腹に力を入れる行きの数量変えるとかありますもうなんかの力を入れて強いわ強い力を込めてもっと明るくならなかったにやと言われるとも出てこなくなったよ沢山入ったとこですねこれで第4子が赤くなるのはこれまで行きの量が足りてなかった自分じゃなかったじゃないかっていうので行ってみましょう早くなりましたねでスタグラムに従ってくると自分の音が帰宅したところ息を吸う量を多くすることで気分良くなることがわかりました待機するのって管楽器重要な今なってどれぐらい剃ったらいいのか分かんなかったり根性論だったりするんでこれは便利だねまあは管楽器にも分からなくなって鳴るんですけど G 計り方とかあると結婚もしないからっても行きたいのでくれてありがとうございました」

ソースコード(LT)

こちらではLTのために使ったスペクトログラムを表示するためのコードを共有します。

ソースコード(文字起こし)

Google Speech To Text APIで文字起こしをする際のコードをここに載せておきます。

GCSのバケットを作っておきます。今回はperry_voiceという名前のバケットを作っておきました。先ほど作ったaudio_only_filtered_441.wavを、このバケットにあげておきます。

APIの認証用のjsonファイルをcloud shellにアップロードしておきます。jsonファイルはGCPのコンソールの認証情報からサービスアカウントキーを作成し、「鍵を追加」の項目からjsonファイルを選択して手に入れることができます。(当然、Google Speech To Text APIのAPIを有効にしておく必要があります。)
1ヶ月につき60分までは無料なので、今回の処理も無料でできます。(GCSは数円かかると思いますが。)

文字起こし結果はGCSのperry_voiceの中にあります。

他の発表から学んだこと

深層学習を用いたJazz音楽の研究や、音楽の理論、DJプレイの自動化の研究、声の変換技術など同じ音という題材ながらも幅広く話を聞くことができました。私だけ入門な内容で申し訳ない気持ちが生じましたが、これからこの分野を深めるのもプレイヤーとして面白いなと思ったので、長く学ばせてもらおうと思います。

懇親会で学んだこと

・NVIDIAのノイズリダクション技術は文字起こしの際にやってみる価値あり
・賃貸の場合、防音の部屋は家賃相場の2~3万円アップ
・口周辺の筋肉の研究をされている方がいる(論文でその分野があるのは知っていたが、実験をやっていた方に会えた)
・みな、条件を揃えるためにマイクや集音の仕方にこだわりがあった。
・いい音の定義にやはり答えはない。

[Python]仕事で使えそうな確率まとめ


本来、正月の企画物としてアイデアを出していたものの、作成する時間がとれなかったのですが、日常的にもっと確率に従って生活できたら面白いなと思って作ってみることにしました。

二項分布

  • やりたいこと
    オフィスにおいて、チームごとにデスク周辺の掃除当番を毎日乱数で決めているとして、1ヶ月で自分がどれくらい掃除をやらなくてはいけないのかを知りたい。あるいは、1回の掃除で済む確率とか、10回掃除する確率とかを知りたい。

    これは、成功確率p(掃除当番になる確率 = 1÷チームの人数 )でのN回の独立なベルヌーイ試行となり、自分が掃除当番になった場合X=1で、そうでない場合はX=0となる。
    N回ベルヌーイ試行した際のXの合計は二項分布従うので、N回試行して、Xがkになる確率は以下のように示される。(チームの人数はHとする。)

    $$ P(Y_N = k ) = {}_N C _k \left (\frac{1}{H} \right )^k \left ( 1- \frac{1}{H} \right )^{N-k} \\
    (0 \leq k \leq N )$$
  • コード


    平均では20営業日で2回は掃除することになりますが、4回以上掃除しないといけない確率が10%以上はあることがわかります。掃除が嫌いな人は、苦しみの上限をイメージできると楽なのではないかと思うので、確率に従い掃除当番を甘んじて受け入れましょう。
  • 参考情報
    scipy.stats.binom
  • matplotlib.pyplot.grid

多項分布

  • やりたいこと
    先ほどのように乱数で掃除当番を割り当て、それをランキングにし、一位を掃除大臣として称号を付与するものとする。10名いるもののうち、上位のメンバーA・メンバーBとで累積の掃除回数がそれぞれ10回、6回としたもとで次の月に掃除大臣が入れ替わる確率を計算したい。
  • コード

    計算したところ、1.2%くらいになった。このケースにおいてAさんがBさんに次の月に掃除大臣の座を奪われる可能性は低いらしい。
  • 参考情報
    scipy.special.comb

誕生日の問題

  • やりたいこと
    あまりにも有名な確率の話で、誕生日のパラドックスとも呼ばれています。誕生日が同じ人が部署内にいる確率が0.5を超える人数は何人か?
  • コード

  • 参考情報
    誕生日のパラドックス
    matplotlib.pyplot.hlines

ファーストサクセス分布

  • やりたいこと
    アルバイトの面接を受けるとして、その合格率が1/30だとして、何回くらい面接しても合格できないものだろうか。苦しみの上限を知りたい。20回、30回面接を受けた時に合格できてない確率はどれくらいだろうか。何回くらい面接を受けたら合格できていない確率を1%未満に下げれるだろうか。
    成功確率がpのベルヌーイ試行を独立に何回も行うとき、初めて成功数までに失敗した回数をXとして、初めて成功した際の試行回数をYとした場合、その確率は\( k \geq 1 \)として、

    $$ P(Y = k) = P(X=k-1)=p(1-p)^{k-1} $$

    と表すことができる。これをパラメータpのファーストサクセス分布と呼ぶ。(幾何分布とも呼ぶ。)

  • コード

    20回面接をしても合格できない確率はおよそ50%。30回面接しても合格できない確率はおよそ36%。合格できない確率が1%を切るのは面接回数が135回になった時。このレベルを追い求めると心が折れる面接回数なのかもしれない。

  • 参考情報
    弱点克服 大学生の確率・統計

負の二項分布

  • やりたいこと
    ある営業がリンゴを売る仕事をしておりタウンページで飲食店に電話をかけているとする、商談の予約が取れる確率(p)が0.05として、商談予約が5(n)回起きるまでに失敗がk回生じる確率について知りたい。いったい、この営業は飲食店にどれくらい商談をお断りをされる覚悟で臨めばいいのか。

    このような事象は負の二項分布に従うので、
    $$ P(X = k) = {}_{n – 1 + k} C _k p^{n-1}(1-p)^k $$
    という確率に従う。
    なお、この分布の期待値は
    $$ E(X) = \frac{n(1-p)}{p} $$
    に従う。

  • コード


    75回くらいの失敗で確率自体のピークがきている。期待値だと95回だが、確率のピークはもっと手前にある。

    累積確率を知りたい場合はこちらのコード


    80%のところで125回なので、5人に一人は期待値の95回を30回くらいは超えてしまってもおかしくないということになる。そして、5人に一人くらいは50〜60回で達成できるラッキーな人もいる。

  • 参考情報
    弱点克服 大学生の確率・統計

ポアソン分布

  • やりたいこと
    先ほど、営業が飲食店に売ったリンゴの中に芋虫が住み付いているケースが500回に1回の頻度で発生し、営業と顧客の間でトラブルに発展するとする。年間に1200の飲食店と取引を行うが、トラブルはどれくらい発生するものなのか。
    顧客とのトラブルが発生する回数がYだとした際に、ポアソン分布に従うとした場合、
    $$P(Y = k) = \frac{\lambda^k}{k!} e^{-\lambda} $$
    となる。ここで、
    $$ \lambda = 1200 \times \frac{1}{500} $$
    とおく。
    なお、期待値と分散はともに
    $$ E(Y) = V(Y) = \lambda $$
    となる。
  • コード


    トラブルが2回の時が確率が一番高い。0回で済む確率も0.1に近いくらいはある。


    5人に一人は4回以上のトラブルに巻き込まれる気の毒な営業がいそうですね。まぁ、リンゴの中は開けてみないとわかりませんから、顧客トラブルがあっても恨みっこなしで頑張りましょう。
  • 参考情報
    弱点克服 大学生の確率・統計

フィギュア集め問題

  • やりたいこと
    現場猫フィギュアのガチャポンに関心があるとする。N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率はどれくらいだろうか?(どの現場猫のフィギュアも出現確率は同じものとする。)
    また、N+1回ガチャポンを回したとき、N種類の全部が揃う確率はどれくらいだろうか?あるいは、m回ガチャポンを回してそろったフィギュアの数の期待値は?そして、Y回目のガチャポンではじめてN種類の現場猫の全部揃う、Yの期待値は?

    ・N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率:
    $$\frac{N!}{N^N}$$

    ・N+1回ガチャポンを回したとき、N種類の全部が揃う確率:
    $$ N \times \frac{(N+1)!}{2} \times \frac{1}{N^{N+1}} $$

    ・m回ガチャポンを回してそろったフィギュアの種類\( X_m \)の期待値:
    m回ガチャポンを回してそろったフィギュアの種類\( X_m \)を以下のように表す。
    $$X_m = Z_1 + Z_2 + \dots + Z_N$$
    m回ガチャポンを引いて、当たらないということは
    $$ P(Z_i = 0) = \left ( \frac{N – 1}{N} \right )^m $$
    ということになる。
    これをもとに\( Z_i \)の期待値を計算すると

    $$ E(Z_i) = 1 \times \left ( 1 – \left ( \frac{N – 1}{N} \right )^m \right ) + 0 \times \left ( \frac{N – 1}{N} \right )^m \\
    = 1 – \left( \frac{N-1}{N} \right )^m $$
    よって、\( X_m \)の期待値は
    $$ E(X_m) = E( Z_1 + \dots + Z_N) \\
    = N \left ( 1 – \left( \frac{N-1}{N} \right )^m \right )$$

    となる。

    ・Y回目のガチャポンではじめてN種類のフィギュアが揃う際のYの期待値:
    i-1種類からi種類になるまでのガチャポンの回す回数を\( T_i \)として、その階差がファーストサクセス分布に従うとする。
    $$ T_i – T_{i-1} \sim F_s \left ( \frac{N – i + 1}{N} \right ) $$
    その場合、Yの期待値は以下のようになる。

    $$ E(Y) = E( T_N ) \\
    = E( (T_{N} – T_{N-1}) + (T_{N-1} – T_{N-2}) + \dots + (T_{2} – T_{1}) + T_{1} ) \\
    = E \left ( F_s \left ( \frac{1}{N} \right ) + F_s \left ( \frac{2}{N} \right ) \dots + F_s \left ( \frac{N-1}{N} \right ) + 1 \right ) \\
    = N \left ( 1 + \frac{1}{2} + \dots + \frac{1}{N} \right )
    $$
  • コード

    N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率:
    0.0384
    N+1回ガチャポンを回した時、N種類の現場猫フィギュア全部が揃う確率:
    0.1152
    m回ガチャポンを回して揃った現場猫フィギュアの種類Xの期待値:
    3.3615
    Y回目のガチャポンではじめてN種類の現場猫フィギュアが揃う際のYの期待値:
    11.4166
  • 参考情報
    弱点克服 大学生の確率・統計

[Python]スペクトグラムで楽器の試奏データを比較(プロ奏者とも)


以前、トランペットのマウスピース選びで音声解析をしてみましたが、今回は楽器選びで行ってみます。
前回の調査でf1値とかでは良い音の判断がつかなかったので、スペクトログラムを出してみたいと思います。
管楽器はそこそこに高い買い物なので、色々と比較見当をしてから意思決定したいですね。


今回の楽器たちはヤマハの800GS(所持品)、Bach 180ML37SP、Bach 180ML37GL、Bach AL190GLです。ヤマハのやつは20万円くらいですが、Bachのは30万円くらいします。ちなみにBachはアメリカの楽器メーカーで、私の感覚値ではあるのですが、音大生はだいたい持っているような気がします。10万円の差は大きいのか小さいのか。そしてプロの演奏家のスペクトグラムとも比較します。

スペクトログラムとは

パワースペクトルの値の高低を画像の濃淡で表現し、縦軸に周波数、横軸に時間を取った座標上にプロットしたものをスペクトグラムと呼びます。スペクトグラムでは、周波数が低いところ(下)から、第1フォルマント、第2フォルマント、第3フォルマント・・・と呼びます。フォルマントは周りの周波数と比べて突出している部分のことを指します。

$$ S(\omega) = \int_{-\infty}^{\infty}s(t) e^{-jft}dt$$

さて、いきなりの数式ですが、フーリエ変換の式となります。\( S(\omega) \)を二乗したものがパワースペクトルとなります。(二乗するは複素平面上で絶対値を取りたいため。)
tは時刻、fは基本周波数(最も低い周波数)、jは虚数単位を表します。ここで登場する\( s(t) \)は定常波と呼ばれ、以下の式で定まります。

$$s(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty}a_n \cos (2 \pi n ft + \theta_n)$$
\( a_n \)は元の信号\( s(t) \)にどの周波数の波がどの程度の強さで含まれているかという情報を表します。
フーリエ変換のモチベーションとしては、周波数の周期関数を周波数の異なる単純な波の重み付き和で表現できることですが、トランペットの音波を色々な周波数に分解できるのは面白いです。良い奏者の周波数がどうなっているのか調べてみるのは面白そう。

Pythonでの計算&比較

音声データはFFMpegで変換とか切り取りをするなどして作成しました。作成方法は簡単で、以前のブログにも記しましたので、試してみたい方はそちらを参考にしてみてください。

まずはヤマハですが、この楽器になれている分、スムーズに吹けています。最初のほうのリップスラーは5つくらいのフォルマントがありそうです。音階が上がるにつれてフォルマントが減っているのがわかります。最後のHigh Bでは細めのフォルマントが3つあります。


次は、Bach 180ML37SPという楽器です。鳴らすのに慣れていないので、ゆっくり吹きました。最初のほうのリップスラーは5つくらいのフォルマントがありそうで、それはヤマハと同じですが、最後のHigh Bで2つくらいしかフォルマントがありません。楽器に十分息が吹き込めておらず、鳴らせていないのだろうと思われます。可視化できるって面白いですね。


次は、Bach 180ML37GLという楽器です。さっきの楽器との違いは銀メッキがなされているかどうかくらいです。最初のリップスラーは同じ感じで5つくらいのフォルマントがありそうですが、リップスラーで戻ったところで5つ目のフォルマントが消えています。ただ、最後のHigh Bで3つほどフォルマントがありそうです。慣れてきたのでしょうか?


次は、Bach AL190GLという高級モデルです。40万円くらいします。より技術力の高い職人が手掛けた楽器だそうです。ベルに一手間加わっていたり、専用パーツがあったりするみたいです。実際、30万円のよりも吹きやすい印象がありました。
最初のリップスラーは同じ感じで5つくらいのフォルマントがありそうですが、最後のHigh Bでフォルマントが2つになっています。やっぱ鳴らすの難しいのでしょうか。


最後に、プロの方の音源です。最初のBで私のグラフと明らかにフォルマントの数が違うことがわかります。楽器を鳴らせているというのはこういうレベルなのかなと思ってしまいました。High Bでフォルマントが4つはありそうです。このレベルになるにはどういう吹き方をすればいいのか、色々研究のしがいがありそう。

参考文献

イラストで学ぶ 音声認識 (KS情報科学専門書)

[Python]データ分析業務で使いそうなコードまとめ(随時更新)

はじめに

仕事で使いそうなPythonのコードを残しておくドキュメントが欲しいなと思ったので、よく使うものをこちらに貯めていこうと思います。まだ19個しかないですが、30個を目標に追記していきます。

フォーマットとしては、
1.やりたい処理
2.コード
3.参考情報のリンク
の3つを1セットにしています。
まずは、自分自身や周りで仕事をしている人が楽をできるドキュメントになればいいなと思って作っていきます。

目次

・重複削除
・階級のデータを作りたい
・再起的にリストをコピーしたい
・ピボットテーブルでアクセスログから特徴量用のデータを作りたい
・データフレームのある列に対して特定の文字列で分割し、その任意のN番目のものを抽出したい
・データフレームの行単位で割り算を行いたい
・列単位で行を足し合わせて文字列を作りたい
・グループのレコードごとにidを割り当てたい
・正規表現で任意の合致した文字列で挟まれた文字列を抽出したい
・複数列の集計値同士で除した値が欲しい
・等間隔の実数の数列を作成したい
・スペース区切りの文字列において、ある単語の前後N個の語を抽出したい
・グループごとに文字列を結合したい
・文字列の置換をdict形式のデータで行いたい
・テキストを形態素解析してBoW形式のデータにして特徴量にしたい
・訓練データとテストデータで特徴量の数が足りない時に不足した変数を追加したい
・アクセスログのユーザーごとのレコードごとの累積和を計算したい
・選択したカラムで一気にダミー変数にしたい
・グループ化して平均値を計算し、プロットし、2軸でそれぞれのデータ数も載せたい
・参考情報

データセット

全体を通じて使うデータセットを生成する関数を用意しておきます。一つ目の
データはkaggleのBoston Housingのデータセットです。

二つ目はアクセスログのデータが欲しかったので、Google Analytics APIを使ってデータを抽出したものとなります。ここではソースコードを載せますが、皆さんは当然見れません。一応、GitHubにデータを上げておきます。

三つ目はテキストデータです。

前処理

重複削除

  • やりたいこと
    重複したデータを除外したい。最初に現れたレコードだけを残したい。
  • コード
  • 参考情報
    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.duplicated.html

階級のデータを作りたい

  • やりたいこと
    任意の範囲に応じたラベルをつけた階級値を作りたい。
  • コード
  • 参考情報
    https://pandas.pydata.org/docs/reference/api/pandas.cut.html

再起的にリストをコピーしたい

  • やりたいこと
    リストをコピーした際に、コピーしたものを変更した際に、もとのものも変更されてしまうので、それを防ぎたい。
  • コード
  • 参考情報
    https://docs.python.org/3/library/copy.html

ピボットテーブルでアクセスログから特徴量用のデータを作りたい

  • やりたいこと
    集計されていないアクセスログをもとに、セッションidごとの触れたページを集計し、Bag of Words形式のデータフレームにしたい。これは機械学習の特徴量に使えるので便利。
  • コード
  • 参考情報
    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.pivot_table.html

データフレームのある列に対して特定の文字列で分割し、その任意のN番目のものを抽出したい

  • やりたいこと
    pandasのデータフレームである文字列の列に対して特定のキーワードで分割し、そのN番目のものをデータフレームの列として持たせる。今回はURLのパラメータの後ろについている文字列を抽出している。
  • コード
  • 参考情報

データフレームの行単位で割り算を行いたい

  • やりたいこと
    各行の合計で、各行を割ることでレコード単位で重み付けを行いたい時や確率を計算したいときに良く使います。データはアクセスログをBag of Word形式したものを使います。
  • コード
  • 参考情報
    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.div.html

列単位で行を足し合わせて文字列を作りたい

  • やりたいこと
    データフレームの列の要素をつなぎ合わせて文字列を作成したい。複合キーを作らないとデータを繋げない時などによく使います。
  • コード
  • 参考情報
    https://www.kato-eng.info/entry/pandas-concat-pk

グループのレコードごとにidを割り当てたい

  • やりたいこと
    アクセスログなどで、任意のユーザーの初回に触れたページや、N番目に触れたページなどを抽出しやすいように、ユーザーごとのログにidを付与したい。
  • コード
  • 参考情報
    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.factorize.html

正規表現で任意の合致した文字列で挟まれた文字列を抽出したい

  • やりたいこと
    正規表現を用いて、文字列に対して「指定した文字列と文字列の間にある表現」を抽出するための処理です。URLのデータや規則性のあるテキストデータに対して特徴量を作りたい時に使えます。
  • コード
  • 参考情報
    https://docs.python.org/ja/3/library/re.html

複数列の集計値同士で除した値が欲しい

  • やりたいこと
    複数の列で集計した値同士を用いて、さらにそこから比率のデータを作りたい。ターゲットエンコーディングを行いたい場合に使える。
  • コード
  • 参考情報
    https://stackoverflow.com/questions/45075626/summarize-using-multiple-columns-in-python-pandas-dataframe

等間隔の実数の数列を作成したい

  • やりたいこと
    パラメータチューニングをする際に、各パラメータについて等間隔の実数列を作りたいときがあるので、その時に使うコード。
  • コード
  • 参考情報
    https://www.pythoncentral.io/pythons-range-function-explained/

スペース区切りの文字列において、ある単語の前後N個の単語を抽出したい

  • やりたいこと
    テキストマイニングをする際に、関心のある単語の前後N文字を抽出したい時がある。もちろん、関心のある単語の周辺の単語というものを特徴量にすることも良いと思われる。ここではMeCabによる形態素解析により分かち書きにする関数と共に紹介する。
  • コード
  • 参考情報
    https://stackoverflow.com/questions/17645701/extract-words-surrounding-a-search-word

グループごとに文字列を結合したい

  • やりたいこと
    テキストマイニングなどで、任意のグループごとのテキストデータをスペース区切りで繋ぎたいときがあります。任意のグループごとのテキストを一気に連結できて便利です。これを形態素解析して特徴量にするのも良いと思います。
  • コード
  • 参考情報
    https://stackoverflow.com/questions/17841149/pandas-groupby-how-to-get-a-union-of-strings

文字列の置換をdict形式のデータで行いたい

  • やりたいこと
    事前に用意した置換のリストを使って、文字列の置換を一気に行いたい。
  • コード
  • 参考情報
    http://omoplatta.blogspot.com/2010/10/python_30.html

テキストを形態素解析してBoW形式のデータにして特徴量にしたい

  • やりたいこと
    任意の日本語のテキストとラベルの組み合わせがあるとして、テキストを形態素解析してTF-IDFを計算し、それを特徴量としてラベルを予測するということをやりたい。ここでは以前集めて、GitHubに載せてあるデザイナーズマンションのデータを使い、MeCabで形態素解析をしたあとに、TF-IDFを計算し、それを特徴量にしてsklearnで簡単に予測をし、リフト曲線を描いている。
  • コード
  • 参考情報
    https://scikit-plot.readthedocs.io/en/stable/metrics.html

訓練データとテストデータで特徴量の数が足りない時に不足した変数を追加したい

  • やりたいこと
    訓練データとテストデータの変数の数が違う時に、そのままではモデルを使えないので、補うために使う。
  • コード
  • 参考情報
    なし

アクセスログのユーザーごとのレコードごとの累積和を計算したい

  • やりたいこと
    アクセスログや購買データなどで、レコードごとの累積の値を計算したいときに使う。時系列に従いどんどん足されていく。
  • コード
  • 参考情報
    https://ohke.hateblo.jp/entry/2018/05/05/230000

選択したカラムで一気にダミー変数にしたい

  • やりたいこと
    ダミー変数を一気に生成したい。名前も一気に付けたい。
  • コード
  • 参考情報
    なし

グループ化して平均値を計算し、プロットし、2軸でそれぞれのデータ数も載せたい

  • やりたいこと
    EDAの時に、グループ化して平均値などを計算するだけでなく、該当するデータ数も2軸でプロットしたい。
  • コード
  • 参考情報
    なし

参考情報

前処理大全[データ分析のためのSQL/R/Python実践テクニック]
はじめてのアナリティクス API: サービス アカウント向け Python クイックスタート
カモノハシペリー
ペリー|フィニアスとファーブ|ディズニーキッズ公式

[Python]音声解析でトランペットの録音データを可視化してみる

はじめに

私はデータ分析以外の趣味としては、トランペットの演奏があります。個人でプイプイ吹いているだけのしょうもない趣味ですけどね。
先日、楽器屋に新しいマウスピース(楽器を口に当てる金属部分、振動させる部分)を買いに行ったのですが、吹いているうちによくわからなくなったので、録音した音声波形などを可視化することで、いいなと思う音がどんな波形をしているのかを明らかにしたいと思います。

進め方

  • 録音データを用意(iPhoneのムービー)
  • FFmpegで動画から音声のみを切り出し、サンプルレートを変え、モノラルにする
  • Pythonで基本周波数(F0)などを可視化

データの準備

今回はチューニングB♭(ベー)の音3種類(Yamaha 11C4、Bach 5C、Bach 3C)と、練習曲(ルパン三世のフレーズ)の音2種類(Yamaha 11C4、Bach 5C)の計5つの音源を扱います。

マウスピースのイメージ

  • Yamaha 11C4:オールマイティー(これまで使っていたやつ)
  • Bach 5C:少し深くて、少し内径が小さい
  • Bach 3C:少し深くて、少し内径が小さい(5Cと3Cの違いは当て心地?)

指標の説明

  • 音声波形:波動の振動の大きさを表す非負のスカラー量を縦軸に、横軸を時間とした場合のグラフを音声波形と呼ぶ。
  • ピッチマーク:単位波形を抽出する基準点でエポックマーカーと呼ばれる。音が鳴った際に1をとり音がない際は0を取る。
  • 基本周波数:信号を正弦波の合成で表したときの最も低い周波数成分の周波数で、音程に該当する。
  • 相関係数:ピッチマークの検出に使う指標でNCCF(Normalized Cross-Correlation Function)というものを使って計算している。自己相関のようなもので-1から1の間を取る。

音声ファイルのあるディレクトリで以下のコマンドを実行します。(私の環境はMacです。)

音声波形などの可視化〜チューニングB♭(ベー)〜

ここでは単純にチューニングベーの音を比較します。まず音源ですが、以下で実際に聴けるようになっています。(耳を痛めないように、最も小さい音で聴いてから調整してください。)
しかし、慣れないマウスピースで音程が不安定になっているのが悲しいですね。

・Yamaha 11C4のチューニングベー

・Bach 5Cのチューニングベー

・Bach 3Cのチューニングベー

以前書いたブログと同様に、pyreaperというREAPERのラッパーとなっているライブラリを用います。REAPERはGoogleが開発していて、C++で書かれています。

音声波形の比較


やはり慣れている分、11c4が安定してそうな波形に見えますね。3Cは結構ブレてそうです。実際の音源でもわかりますね。

ピッチの比較


発音のあったところで1となる指標です。今回は伸ばしているだけなので、何の面白味もないです。

基本周波数(F0)の比較


音程をつかさどるF0に関しては3Cが安定してそうです。聴いたらブレブレなのに。

相関係数の比較


いつも使っている11c4がブレにくいようで、5Cが比較的に上下動が多いように見えます。

音声波形などの可視化〜練習曲〜

結局、マウスピースは5Cを購入して、その後スタジオで練習曲の音を比較しました。先ほどと同様に音源が以下で実際に聴けるようになっています。(耳を痛めないように、最も小さい音で聴いてから調整してください。)私は、実際に録音して、マウスピースでこんなに変わるんだと驚きを隠せないです。チューニングベーを比較するだけではわかりにくいんですかね。私としては、5Cの音の方が好きです。

・Yamaha 11C4での練習曲

・Bach 5Cでの練習曲

音声波形の比較


5Cのほうが上下の大きさが対称に近いような印象を受けます。いい音を捉える参考になるんですかね?

ピッチの比較


おそらく、11c4の音源の息継ぎのタイミングがやや長くなった結果、ピッチに差が出たのではないかと思います。

基本周波数(F0)の比較


比較が難しいですが、最後の伸ばしは5Cのほうが安定しているようです。

相関係数の比較


5Cのほうが上下動をしているように見えます。いいなと思う音の方が上下動しやすいのでしょうか。気になるところ。

追記:スペクトログラムの比較

どうやら音色を考えるうえで、縦軸に周波数、横軸に時間をとったスペクトログラムというものを可視化するのも良いらしいです。参考文献に従って行ってみます。

5Cよりも11C4が赤い線が1つ多いように見えます。幅広く音が鳴っていてモサっとしてしまっているのでしょうか。私としては5Cの方が音が鋭くて好きです。

おわりに

音声波形で音色の比較がうまくできるかどうか試してみたのですが、良いと思われる音とそうでない音で相関係数や音声波形などのデータの可視化において違いが生じやすいのかなという印象でした。そして音色をつかさどるとされるスペクトログラムにおいても、違いがありそうでした。
もっと音源を集めていっていい音の研究をしていこうと思います。あと、音楽の理論の方も細々と学んでいきます。

参考文献

IPythonデータサイエンスクックブック ―対話型コンピューティングと可視化のためのレシピ集
pyreaperで音声データのピッチを掴むためのメモ書き(F0の抽出)
Pythonで音声解析 – 音声データの周波数特性を調べる方法
基本周波数についてのまとめ
イラストで学ぶ 音声認識 (KS情報科学専門書)

BLSTMを用いた文書分類でデザイナーズマンション予測に再挑戦

はじめに

仕事で深層学習を使うことはないのですが、扱える技術の幅を広げていきたいなと思い、BLSTMを用いた文書分類についてkerasでの簡単なサンプルコードをやってみようと思います。データは以前集めた、某不動産紹介サイトの賃貸マンションの設備に関する文書とそれがデザイナーズマンションかどうかのラベルを使います。そして文書内の単語からその文書がデザイナーズマンションかどうかを予測します。前回はAUCで83%だったので、それを超えれると良いですね。

目次

単純なRNNとは

  • モチベーション
    • フィードフォワード型のニューラルネットワークではうまく扱うことができない時系列データをうまく扱えるようにすること。
  •  特徴
    • 入力が互いに関係している(多層パーセプトロンの際は置かれていない仮定)
      • 直訳すると循環するニューラルネットワークとなる。
    • 最初の文の単語が2つ目、3つ目の単語に影響を与える可能性を考慮。
    具体的な関数の形としては、 $$ h_t = \tanh (W h_{t-1} + U x_t) \\\ y_t = softmax(V h_t) $$ で与えられる。
    \( h_t \)は隠れ層の状態を、\( x_t \)は入力変数を、\( y_t \)は出力ベクトルを表している。(他のアルファベットは重み行列)

LSTMとは

  • モチベーション
    • 単純なRNNの勾配消失問題を解決するために提案された手法。
  • 特徴
    • 単純なRNNに置き換えることで性能が大幅に向上することも珍しくない。
      • 時系列データ、長い文章、録音データからなる長期的なパターンを取り出すことを得意としている手法。
    • 勾配消失問題に強い
    • Long Short-Term Memory:長短期記憶
      • 長期依存性を学習できるRNNの亜種。
      • 入力ゲート、忘却ゲート、出力ゲート、内部隠れ層という4つの層が相互に関わり合う。重要な入力を認識し、それを長期状態に格納すること、必要な限りで記憶を保持すること、必要なときに記憶を取り出すことを学習する。
        • 入力ゲート:内部隠れ層のどの部分を長期状態に加えるかを決める。
        • 忘却ゲート:長期状態のどの部分を消去するか決める。
        • 出力ゲート:各タイムステップで、長期状態のどの部分を読み出し、出力するかを決める。
$$
i = \sigma (W_i h_{t-1} + U_i x_t) \\
f = \sigma (W_f h_{t-1} + U_f x_t) \\
o = \sigma (W_o h_{t-1} + U_ox_t ) \\
g = \tanh (W_g h_{t-1} + U_g x_t) \\
c_t = (c_{t-1} \otimes f ) \otimes ( g \otimes i) \\
h_t = \tanh (c_t) \otimes o
$$ i:入力ゲート(input)
f:忘却ゲート(forget)
o:出力ゲート(output)
\( \sigma \):シグモイド関数
g:内部隠れ層状態
\(c_t \):時刻tにおけるセル状態
\(h_t \):時刻tにおける隠れ状態

Bidirectional LSTMとは

  • モチベーション
    •  従来のRNNの制約を緩和するために導入。
  • 特徴
    • ある特定の時点で過去と将来の利用可能な入力情報を用いて訓練するネットワーク
      • 従来のRNNのニューロンをフォワード(未来)なものとバックワード(過去)なものとの2つに分ける。
        • 2つのLSTMを訓練している。
      • 全体のコンテキストを利用して推定することができるのが良いらしい。
        • 文章で言うと、文章の前(過去)と後(未来)を考慮して分類などを行うイメージ。
      • 人間もコンテキストを先読みしながら解釈することもある。
    • BLSTMは全てのシークエンスの予測問題において有効ではないが、適切に扱えるドメインで良い結果を残す。
    • 1997年とかなり歴史があるものらしい
出典:Deep Dive into Bidirectional LSTM

Bidirectional LSTMで文書分類

今回は、BLSTMを使って、マンションの設備に関するテキスト情報から、そのマンションがデザイナーズマンションかどうかを予測します。全体のソースコードはGoogle Colabを御覧ください。

データ

データを確認してみると、
 
今回のデータがテキストとラベルからなっていることがわかる。なお、データ数は1864件あり、そのうち16%がデザイナーズマンションです。

前処理

Google Colab上で形態素解析を行うために、MeCabをインストールする。
 
これでMaCab NeologdをGoogle ColabのPythonで実行できます。
テキストデータの前処理を行うために名詞だけを抽出してリスト形式で返す関数を定義します。
 
ここで、語彙に関するデータを作成します。
 
続いて、単語とそのインデックスからなるdict形式のデータを作成します。
 
最後に、これまでに作成した語彙とそのインデックスのデータを用いて、実際に学習に使う入力データと出力データを作成します。
 

実行

先程作成したデータを訓練データとテストデータに分けます。
ここで、AUCを計算するための関数を定義します。(まるまる拝借しました)
続いて、kerasを用いて、ネットワークを構築し、コンパイルを行います。
  テストデータでの予測精度を確認します。
実際のラベルとテキストデータを見比べてみます。
LSTMの場合は、BLSTMに関する記述(1行)をネットワークから除外すれば良いので簡単に試すことができます。 せっかくなので比較してみたところ、AUCは0.841でした。BLSTMは0.844だったので、若干上回る程度のようです。 以前扱った記事ではデザイナーズマンション分類はAUCが0.83程度だったので、LSTMを使ったほうが精度が少しだけ高いようです。

追記

TensorBoardの検証のloglossを見てみます。 エポックに対してloglossが非常に荒れています。学習曲線としてはセオリー的にアウトなようです。一方で、検証のAUCはエポックに従って高まるようです。 AUCは良くなっているけどloglossが増え続けている。loglossが下がらないと過学習している可能性が高いので、これは過学習しているだけなのだろう。

参考文献

 

2019年に読んだデータ分析系の本の振り返り(21+1冊)

はじめに

2020年、あけましておめでとうございます。年末に自分自身を振り返ろうと思ったのですが、結局データ分析と勉強しかしていないわけで、書籍を振り返ろうと思うに至りました。私の知り合いのデータサイエンティストはだいたい全冊持っているであろうと思われますが、良い本だと思うので思い出していただければ幸いです。

1.『ベイズモデリングの世界』(岩波書店)

基本的に階層ベイズモデルを使って、個体ごとの異質性を考慮した分析手法が提案されています。前半はオムニバス形式で様々な先生がモデルの適用について執筆されており、後半では伊庭先生による階層ベイズモデルの講義になっています。途中でスタイン統計量による縮小推定の話があげられ、柔軟なモデリングのためには「階層化した方が少なくとも望ましい推定量が得られる」という数学的証明を捨てることもやむを得ないと書かれています。

2.『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』(コロナ社)

この本はトピックモデルの教科書というよりも、ベイズ推定の教科書という側面が強い印象があります。途中で出てくる数式は流し読みするのは難しく、最低2冊以上のノートが別途必要になると思います。一度でもLDAのパラメータを導出してみたいという方には良い教科書だと思います。疑似コードが提供されているので、それをもとにRやPythonでコーディングしていけば、一番シンプルなLDAが非常に短い行で実行できてしまうことに驚かれるかもしれません。人間が手を動かして推定アルゴリズムを導出しているからこそ、短いコードで済むということを実感できるはずです。

3.『構造的因果モデルの基礎』(共立出版)

グラフィカルなアプローチで因果推論を扱っている書籍です。Judea Pearl流の因果推論アプローチについて記すことを目的に書かれています。基礎と書かれていますが決して簡単ではありません。ただ、扱われる数学のレベルとしては確率と線形代数がわかれば大丈夫だと思われます。余談ではありますが、1章の相関関係と因果関係の事例紹介で「おむつとビールの話」が都市伝説ではなくきちんと記事としてWall Street Journalという雑誌に掲載されていたことが明らかにされています。

4.『現場で使える!PyTorch開発入門 深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』(翔泳社)

PyTorchを触ったことがないが、深層学習の手法について知っている層を対象とした本です。6章まではGoogleのColabで動かせるのでGoogleに課金することなく深層学習による回帰、CNN、GAN、RNN、Encoder-Decoderモデル、ニューラル行列因子分解をPyTorchで試すことができます。写経したものはこちら。転移学習や高解像度化や画像生成、文章のクラス分類、文書生成、機械翻訳などもできるので、PyTorchでこれくらいの量をコーディングしたらこれくらいのことができるのかという学びや、他の人の書いたPyTorchコードを読みやすくなるなどの便益は十分にあると思いました。

5.『作ってわかる! アンサンブル学習アルゴリズム入門』(シーアンドアール研究所)

会社で行っているPythonもくもく会用に買った本で、scikit-learnを使わずに機械学習のアルゴリズム(アンサンブル系)をコーディングするための本です。pythonのコードについて逐次、細かい解説が行われているわけではないので、1行1行自分でコメントを加えながら写経をしていけば力が付くという本かなと思われます。sklearnはそれはそれで素晴らしいですが、こういう本でフルスクラッチで修行できるのはいいですね。

6.『数理統計学―基礎から学ぶデータ解析』(内田老鶴圃)

統計検定1級を合格された方のブログで紹介されていた教科書です。理系の大学生レベルの数学知識があれば、数理統計学の基礎を学べると思います。中心極限定理の証明や、様々な分布の期待値や分散、様々な分布の性質について数式を用いてしっかり理解することができます。数式もほどよく端折られているので、無論ですがノートが数冊必要になります。各章毎にある練習問題も解くことで力が付くと思います。日本の大学の授業の教科書がこれだったらジェノサイド(再履修者の大量発生)が起きるんだろうなと思ってしまった。

7.『44の例題で学ぶ統計的検定と推定の解き方』(オーム社)

統計の検定に関してだけ扱った珍しい本です。第3部までは統計学の普通の教科書ですが、それ以降であらゆる検定の例題が44件も載せられています。パラメトリックな検定から、ノンパラメトリックな検定まで幅広く扱われています。一番気にいっているのは仮説検定法の分類の表です。これさえあれば、どのデータに対してどの検定を行えばいいかが一目瞭然です。

8.『わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する』(技術評論社)

機械学習の原理を手早く数式を交えて学べる本です。かゆいところに手が届いていると言うか、既出の教科書では捨象されがちな、条件付き確率における2変数以上の条件づけでの表現に紙面を割いていたりしてくれるのが嬉しいです。ある程度数学の話はわかるが、だいぶ忘れているビジネスパーソンには大変にありがたいコンテンツと言えると思います。ベイズ線形回帰に関しても行列を用いた、わかりやすい導出方法が紹介されています。またコラムで紹介されている、測度論にどう向き合えばいいかの著者の見解は参考になります。

9.『Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

R言語とstanを用いてベイズ統計学を入門レベルから学べる本です。各トピックごとにそれなりの紙面が割かれています。例題も豊富にあるので、線形回帰・MCMC・情報量基準・階層ベイズモデルまで、ベイズ統計学を基礎から応用までしっかりと学べると思います。youtubeで著者の講義も配信されているので、留学気分を味わえます。

10.『scikit-learnとTensorFlowによる実践機械学習』(オライリージャパン)

2019年に日本で開かれたML SummitでTFの開発者がおすすめしていた教科書です。前半部分で機械学習の入門から応用までをわかりやすい説明で学ぶことができます。数式は少ないですが、図とソースコード(Python)がちりばめられており、手を動かして理解を進めることができます。後半部分はTensorFlowを用いた深層学習の基礎を同様に手を動かして学ぶことができます。ただ、TFのバージョンも変わってきているので前半の説明をアテにして読むのも良いと思います。

11.『AIアルゴリズムマーケティング 自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)

マーケティングへのデータサイエンスの適用に関する珍しい書籍です。ソースコードはついていないですが、業務で使う際のアイデアが手に入ることもあります。一般的な回帰、生存時間分析、オークション、アトリビューション分析、アップリフトモデリング以外にも、情報検索やレコメンデーションやトピックモデルなどマーケティングながら学際的なトピックも扱われています。レコメンドなどで使われる、ランク学習に関して詳しく書かれた書籍をあまり知らないので、この本はその点においてもありがたい本でもあります。

12.『入門 統計的因果推論』(朝倉書店)

ほぼ全ての章でグラフィカルなアプローチで因果推論を扱っています。例題も豊富なので、一つ一つ丁寧にやれば理解が捗ります。おそらく、例題の多さを含め一番丁寧にd分離性、do演算子、バックドア基準、フロントドア基準に関する説明をしてくれている本なのかなと思いました。グラフでの因果推論に関して初めての人でも、確率さえ知っていれば読み進めることができるはずです。また、途中で操作変数法の紹介もされ、経済学出身者としては読みやすい。ただ、傾向スコアのくだりや、DIDなどのくだりはあまり出てきません。あと、やってないですが章末の練習問題に対するSolution Manualが提供されているようです。

13.『実践 ベイズモデリング -解析技法と認知モデル-』(朝倉書店)

ベイズモデリングを様々な事例に適用する方法がオムニバス形式で記された本です。ワイブル分布や異質性を考慮した二項分布、無制限複数選択形式のアンケートデータに対する手法、トピックモデル、項目反応理論などが扱われています。マーケティングの実務で使える事例が多いように感じました。こちらはサポートサイトでRコードとstanコードが提供されています。あと、appendixにあるプレート表現の見方も参考になります。

14.『機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習などで用いるベイズ推論を扱った教科書です。入門とありますが、入門者は書かれた数式をそのまま見ていても頭に入らないのではないでしょうか。手を動かしてなんぼの本だと思います。ノート2冊は絶対に必要です。たぶん、数式の展開を丁寧に記すと倍以上の厚みの本になると思います。各々のモデルに関してグラフィカルモデルが記されているのや、サンプルコードとしてGitHubにJuliaで書かれたソースコードが提供されているのも良いです。

15.『その問題、数理モデルが解決します』(ベレ出版)

物語形式で、様々な問題に対して数理モデリングのアプローチが紹介されています。途中でマッチング理論やゲーム理論やオークションなども登場することから、経済学出身者も喜ぶ内容かもしれません。社会人になってからナッシュ均衡という言葉が書かれた本は中々出会って来なかった。

16.『ヤバい予測学 ― 「何を買うか」から「いつ死ぬか」まであなたの行動はすべて読まれている』(CCCメディアハウス)

2013年と結構古い本ですが、データ分析を様々な事象に対して適用した事例紹介本です。アップリフトモデリングへの言及もあり、こういったものに関して日本は何年も遅れてブームが来るんだなという実感を与えてくれた本でもありました。appendixに分析事例が147個ほどあげられているのも参考になります。

17.『たのしいベイズモデリング2: 事例で拓く研究のフロンティア』(北大路書房)

主にstanを用いたベイズモデリングによる分析事例が1と2で38本もオムニバス形式で載っています。ほとんどの事例で階層ベイズモデルが扱われています。2では若干マーケティングに近い内容の題材も扱われ、データサイエンティストの人にも嬉しい内容かもしれません。もちろんデータとstanとRのコードがサポートサイトで提供されています。

18.『カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)』(共立出版)

状態空間モデルで時系列予測を行うための手法が記されている本です。RのKFASパッケージが全面に渡って扱われています。トレンドを考慮したり、カレンダー効果を追加したり、共変量を追加したりなど様々なアプローチが紹介されコードも伴っているわけですから、業務でも抜群に役に立ちました。

19.『機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)』(SBクリエイティブ)

自分のいる会社で最低限の数学がわかると思われる若いメンバーに買ってもらうように言っている本です。微積分・線形代数だけでなく、カルシュ・キューン・タッカー条件(最適化数学)に関しても扱ってくれているので、ここで出てくる数学がわかれば大体の論文に立ち向かえると思います。さらに、Pythonの基礎もこれで学ぶことができるので一石二鳥な素敵な本ですね。また、最後の方でスクラッチでアルゴリズムを書くパートがあり、こちらも勉強になります。

20.『機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)』(オライリー・ジャパン)

機械学習における前処理の指針を与えてくれる本です。Pythonのコードが提供されています。例えばですが、「テストデータにだけある、新しい単語は取り除いてしまえばいい」などの細かいアドバイスが何気に嬉しいです。「Effectコーディング」「特徴量ハッシング」「ビンカウンティング」「バックオフ」「leakage-proof統計量」などは読むまで知らないところだったので勉強になりました。

21.『データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング』(オライリージャパン)

データ分析の仕事をする上で最低限必要な知識を幅広く抑えることができる本です。数式は少なく、ところどころ出てくるコードはR言語です。参考文献などがブログだったりするため厳密さがめちゃあるわけではないですが、業務で使う分には問題ないと思います。分類問題において、AUCなどの評価指標だけでなく、予測値自体の探索的分析のすすめなどが書かれており、参考になりました。また、特徴量エンジンとしてのk-NN法の話も面白いと思いました。

[+α]『プログラマのためのGoogle Cloud Platform入門 サービスの全体像からクラウドネイティブアプリケーション構築まで』(翔泳社)

Google Cloud Platformを初めて触るデータ分析者にはちょうど良い本です。説明もわかりやすいので、いきなりアカウントを作ってドキュメントを解読するよりかは戸惑いは減るはずです。この本を土台に、GCS・GCEを駆使してML系のAPIを呼び出して使うなどの最低限の操作は私でもできるようになりました。GCPの画面や機能もどんどん変わっていくので書籍を買ってもアレなんですが、歴史的な背景も若干記述されているので、それはそれで勉強になります。ただ、エンジニアにこの本を買うべきか聞いた際にネガティブな意見があったのですが、たぶん現役プログラマからすると簡単過ぎるからなんだろうなと思います。

終わりに

2019年もぼちぼち勉強できましたが、2020年もこれまで同様にノートとペンを大事にする勉強を続けていき、コーディングも分析ももっともっと数をこなして会社や社会に求められるようなデータ分析官を目指していこうと思います。あぁ、英会話などの勉強をする時間を作るのが難しい。

[Python]機械学習などでテキストデータを特徴量にする際のソースコード集

テキストデータの特徴量化について

仕事ではテキストデータを多用するので、機械学習などで扱うためにテキストデータを特徴量にするためのアプローチを色々と整理してソースコードを残しておきたいと思います。今回はあくまでも私の知っているものだけなので、網羅性はないかもしれませんが悪しからず。
(2019/08/18 追記)Stackingをカジュアルに行えるvecstackというモジュールを用いた予測も試してみました。下の方の追記をご覧ください。

アプローチ

テキストデータを特徴量にする際のアプローチとしては、以下の3つが良く使っているものとなります。
・単語ベース
・クラスタ、トピック、分散表現ベース
・文書間の類似度ベース

今回扱うデータ

ひょんなことから、昨年10月くらいに取りためたマンションの施設情報のテキストです。

緑色が印象的な某不動産紹介サイトをクローリングしました。全部で1864件ほどの文書数となります。

加えて、デザイナーズマンションかどうかのフラグを作成しました(17%くらいがデザイナーズマンションの割合)。これでもって、マンションの施設情報からデザイナーズマンションかどうかを分類できるかチャレンジしたいと思います。
ここにデータを置いていますので、興味のある方はご利用ください。

今回扱うモデル

ランダムフォレストです。10foldsクロスバリデーションによるAUCの結果を各手法のスコアとして扱います。

こちらは、任意の手法に関して10foldsクロスバリデーションで実行し、AUCのグラフを生成してくれるソースコードです。主にscikit-learnのサイトに載っているものです。引数のclassifierをsklearnの任意のモデルのインスタンスで渡せば動きます。

単語ベース

シンプルに単語をそのまま特徴量にするというものですが、文書によっては単語数が多すぎて収集がつかないと思います。そこで単語を簡単に選択できるDocumentFeatureSelectionというパッケージを利用します。

このパッケージでは
・TF-IDFベースの特徴量選択
・PMI(Pointwise Mutual Information)ベースの特徴量選択
・SOA(Strength of association)ベースの特徴量選択
・BNS(Bi-Normal Separation)ベースの特徴量選択
を行うことができます。

まずは今回のベースラインとして、単語のカウントベースでの特徴量を扱いたいと思います。
その前に、GitHubに上がっているデータに対して以下のように簡単な前処理をしておきます。

ようやくベースラインの予測となります。以下のコードを実行すると、ROCが描かれた図がJupyter上で表示されます。

AUC82%というのはベースラインとしてはなかなか強敵なのではないでしょうか。

さて、本題の特徴量選択パッケージの適用をするためのソースコードを以下に記します。

以上のソースコードを実行すれば、tf_idf_scored_df、pmi_scored_df、soa_scored_df、bns_scored_dfにスコアを付与された単語のリストが手に入ります。

ここでは各スコアに関してアドホックに閾値を設けて、特徴量として利用することにします。

TF-IDFベースの特徴量選択

PMIベースの特徴量選択