OpenCVで遊んだり、Kerasで「かものはしペリー」と「イコちゃん」の画像分類をするの巻

はじめに

正月は帰省もせずに暇だったので、普段の業務ではほとんど扱わない画像関連でブログを書いてみようと思います。
まず、よくありそうなOpenCVの適用や、DNNを用いた画像分類などを行います。扱う画像は当然、かものはしペリーのぬいぐるみです。ブログのタイトルやドメイン名になっている通り、私はかものはしペリーの愛好家です。
この記事の読者層は私で、いつか仕事で使うのをイメージして備忘録として残す感じです。
説明があまりありませんが、気分次第で追記するかもしれません。あしからず。

OpenCVで行うあれやこれや

インストールする

これで入ります。

画像を読み込んでRGBにする

明度値の書き換え

量子化

グレースケール

画像ヒストグラム


HSV変換について
・Hue:色の種類(例えば赤、青、黄色)
・Saturation:色の鮮やかさ。色の彩度の低下につれて、灰色さが顕著になり、くすんだ色が現れる。
・Brightness:色の明るさ。

明度変換

フィルタ処理

特徴点検出

画像の中から特徴的なポイントを抽出するアルゴリズム。

マッチング

特徴量抽出

・BFMatcher:Brute-Force matcher(総当たりで特徴点を比較している。類似度は数値が低いほど類似しているとみなす。)
・AKAZE:日本語の『風』から命名された手法。二つの画像のキーポイントを発見するために使われる。変化への耐性が強いとされる。非線形拡散方程式を近似的に解いて、非線形スケール空間の特徴量を検出している。

これがベースとなるペリーの画像。

以下では、このペリーとの類似度を計算している。

画像分類

ペリーとイコちゃんの画像からペリーかどうか判定したい。
今回はペリーの画像27枚、イコちゃんの画像20枚を用意しました。

今回はKerasにある、VGG16というモデルを使って分類を行います。ハイパーパラメータはAdamです。reluとかドロップアウト層とかは参考文献のまんまを使っています。

ファイルの読み込みです。
trainというディレクトリにperry_とかikochan_とかからなる画像ファイルがある想定です。

画像からの学習用のデータの作成。

学習の実行。


テストデータにモデルを当てはめて推論し、精度をみます。今回はテストデータとして別に、ペリー9枚、イコちゃん7枚の計16枚の画像を用意しました。
果たしてペリーとイコちゃんを識別することはできるのでしょうか。

・Accuracy:0.688
・Recall:0.857
・Precision:0.6
・F1:0.705

どうなんでしょう。まずまずなんでしょうか。

間違えたやつがどれか確認

仕事で使う機会があるようなないような画像の世界ですが、計算資源が大事だなと思いますね。

参考情報

テキスト・画像・音声データ分析 (データサイエンス入門シリーズ)
Display OpenCV Image in Jupyter Notebook.py
Python でグレースケール(grayscale)化
OpenCV で画像のヒストグラムを作成する方法
HSV色空間
色空間の変換
画像フィルタリング
OpenCV: 特徴点抽出とマッチング
キーポイントとマッチの描画関数
OpenCVのAKAZEで顔写真の類似度判定をやってみた
総当たりマッチングの基礎
OpenCV3でAKAZE特徴量を検出する
KAZE Features
小惑星画像の対応点決定を目的としたSIFTとAKAZEの性能比較
線形・非線形拡散方程式の差分解法と解の可視化
AI技術を魚種の画像分類に応用してみた!
ベクトルの内積や行列の積を求めるnumpy.dot関数の使い方
Python | Count of common elements in the lists
[解決!Python]リストの内包表記と「if」を組み合わせるには
Convert png to jpeg using Pillow
ValueError: Failed to convert a NumPy array to a Tensor (Unsupported object type int) in Python
Pythonで文字列のリスト(配列)の条件を満たす要素を抽出、置換
ImportError: cannot import name ‘adam’ from ‘keras.optimizers’
pandas.DataFrame
Pythonでenumerateとzipを組み合わせて同時に使う

各業界でのデータサイエンスの活用について調べてみた(随時追加)

仕事で、いろんな会社でデータサイエンスってどう使われているのですか?と聞かれることがあり、自分としてはなんとなくしか掴めていないな、知ったかぶりしたくないなと思うところがあったので、やや手厚くリサーチをしてみようと思いました。

2021/12/13の段階では9つの市場しかないですが、最終的には30市場を目指します。

【2021/11/27追記】
公開したところ、それなりにこの記事に関心を持ってくださった方が多かったようなので、少しずつ事例を埋めていこうと思います。

業界別にデータサイエンスの事例を集めてみる

・業界としては、市場規模マップの上位30の市場としてみる

事例の調べ方

各業界の売上1~2位の企業名とデータサイエンス系キーワードを掛け合わせGoogle検索ないしGoogle Scholar検索を行う。

  • Google検索を自動で行う。
    • 「かものはし製パン 機械学習」
    • 「かものはし製パン 人工知能」
    • 「かものはし製パン データサイエンス」
  • 30業界2社の3つのキーワードごとの検索結果上位30位を取得
    • 計5400の検索結果から求人などを除外
    • 不要なものを除外したリストから目視で精査、まとめ

検索自動化のソースコード

一応、Google Scholarの検索自動化のためのPythonコードを置いておきます。クローリングのお作法などは「Pythonクローリング&スクレイピング[増補改訂版] -データ収集・解析のための実践開発ガイド-」のような本で身につけておくと分析関連業務で楽になります。

総務省の白書に関してコメント

日本企業におけるデータ活用の現状」という総務省が出している2020年3月のアンケートに基づいたレポートがあるのですが、なかなか良い統計が載っていて、日本全体の会社におけるデータ活用度をある程度掴むことができると思いました。

どうやら大企業の半数近くは経営企画や製品企画、マーケティングにデータを活用していて、大企業の2割の企業で機械学習を用いたデータ活用を行なっており、大企業の5割の企業でデータ分析専門の部署があるようです。
データ分析が、かなり社会に浸透しているのだなと思われました。
そして、多くの企業の関心事がデータの質になってきているようです。データ分析で重要な観点にしっかり関心が高まっているのは良いことですね。

業界ごとの取り組み

以下では業界ごとに事例などをまとめていきます。結構長くなりそうです。

自動車・同附属品製造業

数理工学や群体数理、行動経済学など幅広い領域を研究対象としている。

  • シミュレーションを活用した制御ソフトの自動検証技術
    • サーチベーステスト(SBT):自律走行ロボットが充電ステーションへ帰還する際に、自己位置の推定に誤差が生じる課題を抱えていたが、SBTを用いることで推定誤差が生じる環境要件を特定してソフトの誤差80%減を実現していたりする。
  • 機械学習による性能テストの効率化
    • エンジンの燃焼に関するシミュレーションの実行時間の短縮を機械学習で行う。シミュレーションの効率化が積極的に行われている。
  • 歩行者保護性能の評価試験
    • ボンネットの設計効率化に向けて、CAE(computer-aided engineering)で構造解析を行う前に、ある程度の歩行者保護性能をAIで判定できるようにする取り組みが行われている。
  • Positive-confidence(Pconf)分類
    • アプローチとしては、正例の分布とその信頼度だけで、分類ミスによって生じる損失の期待値を最小化する手段。
    • ドライバーの心電情報から眠気を推定する際に用いた
  • 深層強化学習
    • 自動運転の際に、「ぶつからない」ことを学習する分散機械学習
  • 車内にカメラやマイク、脈拍センサーなどを置き、乗員の行動や会話、生体情報などを把握する開発
    • 自動運転車からデータを集めてAIで解析すると、人の行動や嗜好、健康状態、さらには感情まで推定
    • 感情エンジンを開発し、ドライバーの表情や声の調子からストレス状況を判断して安全運転をサポートしたり、ライフスタイルでの嗜好を学習して、状況に応じた選択肢を提案する。

建設

  • スマートハウス
    • 既存のIoT機器やWebサービス、建物内に設置するデータ収集サーバを連携し、個々の機器・サービスを統合するシステムを構築し、家電の遠隔操作や省エネアドバイスなどの生活サービスをユーザーの好みに応じてカスタマイズできるシステムや、音声認識による機器操作、情報配信等のサービスにつなげる。温度センサとエアコンや電動窓を連動し真夏の寝苦しさを解消するなども。
    • HEMS(Home Energy Management System)の規格のもと、IoTデバイスを起動させるケースが多い。音声認識を用いたアプリケーションや気象や様々なセンサーを基にした家電の最適化などが扱われる。
  • 施設内の映像解析
    • 施設内に設置したカメラ映像から利用者のマスク着用や施設内カフェテリアの混雑度を自動で検知する実証実験を行っている。モニターでリアルタイムの映像を確認できるが、個人情報保護の観点から人物はシルエット表示されないとのこと。なお、実証実験で利用するカメラ映像は2カ月間保存・利用し、期間終了後に消去する。
  • マンションにおける電力供給システムの最適化
    • 家庭向けの発電設備の稼働を、過去の電力使用データなどから最適制御を行い、遠隔でコントロールする。生活リズムごとに家庭をクラスタリングしている。
  • 在宅健康モニタリング
    • 在宅健康モニタリングと早期発見システム(Early Detection System = EDS) の推進に特化したプログラムを構築のために大学と共同研究

不動産

  • 不動産価格の算出のための機械学習システム
    • 所有マンションの推定成約価格を機械学習で算出して提示するシステムがWebサイトで提供されている。過去の制約事例のデータから予測している。
  • 賃料推定のための機械学習
    • 賃貸仲介データ等が少ない都心商業施設を対象に賃料を推定できるシステムが開発されている。
  • 間取りの検索システム
    • ユーザーに適した間取りの物件をレコメンドできるシステム。画像処理技術を間取りのデータに対して適用している。洋室や和室なども分類されている。

医療(医薬品の卸)

  • 機械学習を用いた、工場における医薬品の物流最適化(医薬品の卸の会社)
    • 商品(医薬品)の保管・払い出し・仕分けなどの作業を自動化する、ピッキングロボットを開発し、人間の2〜5倍の作業量をこなしている。
  • 医薬品の出荷予測
    • 過去数年分の医薬品の出荷データをもとにGCPで機械学習を行い、余剰在庫や欠品リスクを回避させようとしている。
  • 医薬品の配送の最適化
    • 医薬品の販売、物流、商品、需要トレンドのデータを用いて機械学習を行い、注文数や配送発生確率、納品時の滞在時間を予測できるようにした。

生命保険

  • 保障設計のレコメンドシステム
    • 保障設計予測モデルを構築し、顧客の意向に基づいて保障プランをレコメンドする。過去の大量なデータを用いて、保険において重視したい保障内容や保険料の予算をもとに、保障プランを自動生成。
  • 保険の保障見直しのためのレコメンドシステム
    • 未熟練な保険のプランナーの業務を支援するアプリケーション。年齢、性別、家族構成はもちろん、保険の保障設計書などのテキストも学習用データにしている。
  • Human in the loop
    • 機械学習システムは完璧なものではないので、人間の仕事を補助するなど、飛行機における管制システムのような存在として共存するという考え方。営業職員が多い生命保険会社では重要な観点。営業職員がベテラン層の営業職員と比肩するコンサルティングをより適切なタイミングで実施しやすくするリコメンド機能などがそれにあたる。
    • 最適な提案のタイミング、最適な活動レコメンド、最適な教材・資料探しなど、営業活動の支援のための分析に取り組む事例がそれなりになる。最適な顧客訪問、最適なメッセージなどを見つける。
  • 入院日数予測モデル
    • 保険加入時の健康状態(告知事項や健診結果など)と保険給付実績を分析し、加入時の健康状態から将来の生活習慣病における入院日数を予測するモデルを開発。このモデルを用いてリスクを評価できるようになり、これまで引き受けれなかった顧客からの申し込みを引き受けることが可能になった。
  • 営業の教育支援ツール
    • AIを相手にしたロールプレイング演習などを行い、新人の育成を支援するツールなどが作られている。
  • 保険の引き受け査定の自動化
    • 健康診断結果報告書や人間ドック成績表をもとに、OCRでデータを抽出し、それらのデータをもとに機械学習を用いて自動での医務査定を行えるように訓練。これまで2〜3日かかっていた査定業務が3分以内に完了するように。

    外食

    • 配送員向けアプリ
      • 機械学習を用いて、ベテラン配達員のデータから学習し、配送員にとって効率の良いルートを提案するアプリ。
    • 機械学習を用いたクーポン配布
      • 顧客ごとの来店確率を計算し、その確率が上位の顧客にクーポンを送付するという施策を試験的に展開している。コンバージョンレートが通常の3倍になったケースもあるとのこと。天候や気温などのデータとも連動しているとのこと。
    • 注文の需要予測
      • 過去の注文履歴を分析し、宅配需要を10分単位で予測したり、配達ルートの最適化を行う取り組み。フードロスを減らすこともモチベーションとしてある。
    • レコメンドとしてのデジタルメニュー
      • 顧客の登録情報や来店時刻や立地ないし季節・天候に応じて、客ごとにメニューが変わるようなレコメンデーションを試験導入しているとのこと。

      物流

      • 強化学習を用いた船舶の航路の最適化
        • 膨大な数の航海シミュレーションを通して、徐々に最適な避航行動を学習し、様々な状況下で安全性と経済性に優れた避航操船行動を選択できるプログラムの開発などが進められている。航路作成は人間だと3時間かかるらしい。
      • 画像解析による船舶の見張り支援システム(衝突回避のため)
        • 船上で撮影された映像に対して、サーバー上で画像を解析し機械学習を実施後、その物体との距離を求めてから、遠隔で船上ソフトウェアを更新することで、利用しながら認識率等の性能が向上する仕組み。
      • 船体の異常検知の自動化
        • IoTで航海データ・機関データを取得し、船舶の運航管理者が船舶の異常を検知できるように機械学習の導入を進めている。これまでは人間が24時間体制で管理していた。
      • 業務量の予測
        • 宅急便の拠点における、数ヶ月先の業務量を予測するため、機械学習モデルを毎月アップデートし、需要に応じた効率的な経営資源の最適配置とコスト適正化を推進することを目指している。
      • 作業員の業務日報のOCR
        • 紙の作業日報の情報のシステム入力を自動化するために、機械学習を用いたOCR技術を適用している。

      BtoC-EC

      • 商品の非構造化データを構造化データにする、あるいは商品情報の自動生成、自動判別
        • 機械学習(CNN、RNN)や自然言語処理をもちいて、商品から情報を抽出したり、レビューを分析したり、商品のカテゴリを分類したりしている。商品の画像のデータからカテゴリを予測するなども行われている。
        • 出品時に登録した商品の写真画像から、商品の内容などの入力候補を提示する機能などが開発されている。ディープラーニングとk近傍法などを組み合わせた事例がある。
      • 潜在的に購買する可能性があるユーザーの推定
        • 機械学習を用いて金融サービスを使って購買をしそうなユーザーを分類するモデルを作成し、予測スコアが一定以上のユーザーに介入するなどを行っている。
        • 機械学習を用いた独自のアルゴリズムで消費行動を解析することで、購買の見込みがあるユーザーを抽出するなどし、商品の広告配信に活かしている。
      • 動画コンテンツの翻訳
        • アテンションの伴ったRNNでの機械学習を行い、中国語を英語に変換している。
      • 商品推薦システム・パーソナライズドクーポンアルゴリズム
        • 商品データ、購買履歴データ、閲覧履歴データを対象とし、古典的な相関ルール抽出、あるいはMatrix Factorizationを適用させたり、商品やユーザなどを固定長のVector に変換して距離を計算するEmbeddingなども行い、RNNなどの深層学習アルゴリズムを適用している。
      • チャットボットによるサービス利便性の向上
        • ユーザーが知りたい問い合わせの多くにチャットボットが回答している。FAQを情報源としてNLPなどを用いてその中から一番近い答えを返している。
      • 多椀バンディットアルゴリズムを用いたABテストの最適化
        • T.B.A
      • 在庫数の最適化のための商品需要予測
        • T.B.A
      • 異常検知
        • T.B.A
      • 商品価格の推定
        • 出品された商品に対して、Factorization Machinesというレコメンドシステムでよく利用される手法の一種であるDeepFMを用いた価格の予測が行われている。機能によっては端末のアプリ側で計算を行うことから、エッジコンピューティングのためのモデル軽量化なども進められている。
      • 不正検知
        • 違反出品検知システムとして、Human-in-the-Loopを実践。人間が介入するため、判定するべき対象とその他との境界に位置するデータに対してのアノテーションを継続的に入手している。アルゴリズムとしてはDNNを用いている。

      電力

      • メッシュ単位での電力需要予測
        • メッシュ単位での気象予測データを活用して電力需要の予測精度を高めるように、スパースモデリングやアンサンブル学習の機械学習技術が使われている。解釈可能性も重視されるので、一般化加法モデルなども利用されている。なお、電力需要量は、カレンダー条件(時間帯,曜日など)と気温や日射強度などの気象条件からおおむね推定できるとされている。機械学習系の手法もあれば、周期性に着目し、状態空間モデルで予測する事例もある。
      • 機械学習による石炭火力発電所における燃焼調整の体系化
        • 運転データやオペレーションとメンテナンス(O&M)の知識をもとに、設備の燃焼効率最適化を機械学習で行う。
      • 機械学習を用いたダム操作最適化システム
        • ニューラルネットワークを利用したダム操作最適化システムの研究がされている。雨量予測、河川流量予測、ダム操作最適化の三つの領域での最適化がなされている。最適化操作のために強化学習も採用されている。学習期間の報酬を最大化するようにネットワークの重みを修正している。
      • 異常検知
        • ヘリコプターやドローンなどで撮影した送電線の動画データを基に、軽微な損傷や異物の付着といった異常を見つけ出すために深層学習が使われている。誤判定のチェックに多少の人手をかけても、十分な業務時間の短縮効果を得られるとのこと。
        • 工場設備の温度変化のデータを画像化して、CNNで特徴量を作り、不具合や異常の分析を行うなどしている。
      • 太陽光の発電量予測
        • T.B.A
        • T.B.A
      • 大量のデータを駆使しての設備の劣化予測診断
        • 高クロム鋼配管溶接部のクリープ損傷評価を深層学習(AlexNet)やサポートベクターマシンを用いて行っている。学習用のデータが少ないことから、転移学習も行っている。
        • T.B.A
      • 燃料プロセスの最適化
        • ボイラ燃焼調整の最適化等
        • 熟練技術者が半日程度かけてつくっている燃料の運用計画を数分で行えるように数理最適化などの取り組みがされている。
      • 生活リズムの変化を察知
        • スマートメーターにある電気使用量の時系列変化を捉えて、独居している人の異常を親族に通知する仕組みを構築している。
        • 高齢者を介護する負担を減らすための実証研究などが進められている。
      • 電力の需要予測
        • 特徴量としては、カレンダー情報・気象情報それらを組み合わせた相互作用項で、重回帰モデルと一般化加法モデルに対してアンサンブル手法を適用している。
        • k-nearest neighborでの距離的な離れ具合でもって、異常な需要を当てるなどの取り組みがされている。それと重回帰ベースのモデルとを合わせて予測している。

      電気通信

      • 記事リンク
      • 論文リンク

      今さら”Recommendations as Treatments: Debiasing Learning and Evaluation”を読みソースコードを眺めるなど

      はじめに

      最近、『施策デザインのための機械学習入門〜データ分析技術のビジネス活用における正しい考え方』が流行っていると思いますが、私も少しずつ読んでいます。

      2019年ごろから勉強会などで知る機会のあったCountefactual Machine Learningですが、仕事で機械学習を使えば使うほどにその必要性を強く感じます。
      書籍は体系的にまとまっていて非常に良いのですが、せっかく関心を持ったのでその源流と思われる論文やソースコードを漁ってみようと思いました。

      今回、深ぼるのはこちらの論文です。
      Recommendations as Treatments: Debiasing Learning and Evaluation
      https://arxiv.org/abs/1602.05352

      情報を漁る

      arXivの方はAppendixついてなかったのでこちらの方がいいと思います。
      https://www.cs.cornell.edu/people/tj/publications/schnabel_etal_16b.pdf

      著者のスライドもありました。
      https://www.cs.cornell.edu/~schnabts/downloads/slides/schnabel2016mnar.pdf

      まずは読む

      以下に読んだ際のメモを箇条書きします。

      • レコメンドのデータは大体偏っている。偏ってアイテムの評価がされる。
      • レコメンドに因果推論のアプローチを適用する。
      • 傾向スコア(ランダムっぽさ)を用いた学習の重み付け
      • バイアスの最小化?のアプローチ
      • 行列分解による欠損値補完に対する重み付けをした。
      • 色々なデータで試した。
      • IPS(Inverse Propensity Score)で重み付けしたら誤差の不偏推定量を手に入れれる。
      • 評価方法はMAE(Mean Absolute Error)とかDCG(Discounted Cumulative Gain)など
      • 比較は一様に評価のデータが観察されると仮定したナイーブなモデルとで行う。
      • 行列分解の次元数dとか、正則化項のラムダはハイパーパラメータとして扱う。
      • 行列分解かつIPSで重みづけしたモデルはロバストに良い性能になった。

      既出の参考資料を読む

      ソースコード読む

      コーネル大学のサイトに論文のソースコードがあったので、そちらの解読を試みます。
      https://www.cs.cornell.edu/~schnabts/mnar/

      そんなにソースコードは多くなかったです。

      train.sh

      学習のバッチ処理

      train.py

      python ../bin/train.py –lambdas 0.008,25 –numdims 40 –ratings data/train.ascii –propensities data/propensities.ascii

      • コマンドライン引数たち
        • lambda:正則化パラメータ。デフォルトは’0.008,0.04,0.2,1,5,25,125′
        • numdims:行列分解の次元数。デフォルトは’5,10,20,40′
        • ratings:訓練用データの指定を行う(ASCII formatの行列)
        • propensities:傾向スコアのデータの指定を行う(ASCII formatの行列)
        • metric:デフォルトは’MSE’、MAEかMSEを選べる

      MetricsのMSE/MAE関数を呼び出し、
      Expt3のlearn関数を使って学習を実行する。

      Metrics.py

      • MSE関数
        • ITEMWISE_METRICSを呼び出している
      • MAE関数
        • ITEMWISE_METRICSを呼び出している
      • ITEMWISE_METRICS関数
        • MSEかMAEかでnumpyのデータの形式を指定している。
        • SET_PROPENSITIES関数を呼び出している
        • 観測された評価と予測した評価の差分にIPWを掛け、観測誤差とする
        • 累積の観測誤差を出す
        • IPWの合計値で累積の観測誤差を割ったglobalNormalizerを計算して返す
        • ユーザー単位での観測誤差の推定値をユーザー数で割ったuserNormalizedMetricを計算して返す
        • アイテム単位での観測誤差の推定値をアイテム数で割ったitemNormalizedMetricを計算して返す
      • SET_PROPENSITIES関数
        • IPWがあればそれを使い、ない場合は出現割合を返している
        • 観察されたデータを考慮して補間を行う(numpy.ma.arrayのmask)

      Expt3.py

      Expt2のINIT_PARAMS、MF_TRAIN、FINAL_TRAINを実行する
      MFのPREDICTED_SCORESを実行する

      • 処理の流れ
        • もろもろの引数で学習の際の設定を受け取り
        • 傾向スコアの逆数からIPWを計算し
        • そのスコアを学習用データで補間し
        • ランダムに並び替える
        • 4-foldsでのクロスバリデーションを行う
        • 学習を並列処理を行う
        • 予測スコアから一番いいスコアのモデルを選ぶ
        • 行列分解での機械学習を行う
        • 行列分解での予測値を計算する
        • 予測結果を返す

      Expt2.py

      MFのGENERATE_MATRIX、PREDICTED_SCORESを実行

      • MF_TRAIN関数
        • IPWを受け取って、上限下限などを調整
        • それをMFのGENERATE_MATRIX関数に適用
      • FINAL_TRAIN関数
        • MF_TRAIN関数と似たような処理
        • 交差検証してベストだったパラメータに適用する想定
      • INIT_PARAMS関数
        • 評価行列の欠損値を補間している
        • 行列特異値分解を行っている
        • その結果をstartTuple(初期値)として与えている

      MF.py

      • PREDICTED_SCORES関数
        • バイアスがあるかどうかで分岐
        • バイアスがあれば、元のスコアにユーザーのバイアス、アイテムのバイアス、グローバルなバイアスを足し合わせる
        • バイアスがなければそのまま返す
      • GENERATE_MATRIX関数
        • 指標がMSEかMAEかで分岐
        • MetricsのSET_PROPENSITIES関数を実行
        • ユーザー単位での標準化やアイテム単位での標準化、全体での標準化のための計算をする
        • 標準化の方法で分岐、傾向スコアを標準化する
        • バイアスのモードで分岐
        • 学習のための初期値としてのデータを設定
        • 目的関数の最適化(最小化)を行う
        • その結果を返す
      • Mat2Vec関数
        • 行列をパラメータのベクトルにする
      • Vec2Mat関数
        • パラメータベクトルを入力する
        • ユーザーベクトル、アイテムベクトル、ユーザーバイアス、アイテムバイアス、グルーバルバイアスを返す

      test.sh

      テストのバッチ処理

      test.py

      python ../bin/test.py –test data/test.ascii –completed completed_ratings.ascii

      • コマンドライン引数たち
        • test:テスト用の評価データを指定する(ASCII formatの行列)
        • completed:評価を埋めたい行列のファイル名を指定

      MSEやらMAEの結果を返す。

      所感

      まず論文の内容ですが、大体のビジネスデータはRCTな状況で生み出されたものではないので、偏って観測されるのは当然ですが、その結果としてレコメンドのアルゴリズムの性能が落ちてしまうということが論文に記されていました。性能を落とさないためにも、傾向スコアを用いた重み付けによる機械学習が大事に思いました。ただ、傾向スコアの推定自体の精度も大事に思われるので、銀の弾丸ではないですが、バイアスから抜け出るヒントをいただけたという気持ちです。
      そして、思ったよりもソースコードが長いと思ったのと、普段あまり使わないNumpyの関数を知ることができたのもよかったです。
      理論の理解だけでなく、それを実現するためのソースコードを眺めるのも大事ですね。今後も論文をただ読むだけでなくソースコードの海に飛び込んでみようと思います。

      参考リンク

      numpy.reciprocal
      maskedarray.generic
      numpy.ma.getmask
      numpy.clip
      scipy.sparse.linalg.svds
      numpy.ma.filled
      scipy.optimize.minimize

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

      はじめに

      先日、近所の幹線道路の信号付き横断歩道を横断しようとした際に、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は結構ブレてそうです。実際の音源でもわかりますね。

      ピッチの比較