異常検知手法を用いたバズり検知

この記事は[レコチョク Advent Calendar 2021](https://qiita.com/advent-calendar/2021/recochoku) の 13 日目の記事となります。

今回は、異常検知手法を用いて、「アーティストのバズり検知ロジック」を作る過程を記事にしたいと思います。

## はじめに

現在、 TikTok やサブスクサービス等で楽曲が**バズり**、一気にスターダムまで駆け上がる例が出てきています。
バズりを素早く検知することができれば、次の施策に向けて予算を確保したり、企画を立てたり、素早く動き出すことができます。

貴重なチャンスを見逃さないためには、バズりの検知を自動化することが望ましいです。
アーティストにとって**バズる**とは、楽曲が通常時よりも多く聞かれるようになることであり、音楽に関心の高いユーザが多い音楽サブスクサービスでその動きが顕著に出ると考えられます。
今回は弊社サブスクサービスにおける再生数の推移を時系列データとし、異常検知手法を適用することでバズりの自動検知を試みようと思います。

## 時系列データの異常検知

一口に異常検知と言っても、検知したい異常の形は様々です。
今回は[こちら](https://www.albert2005.co.jp/knowledge/machine_learning/anomaly_detection_basics/anomaly_detection_time)を参考に 3 つのパターンを紹介したいと思います。

### 外れ値検知

- 普段は起こり得ないようなデータ点を検知する手法
- データの値自体が異常な点を見つける

![外れ値の例](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2269973/fc7fc91a-0c5a-591f-8bc9-5008016ae645.png "外れ値の例")

### 異常部位検出

- 異常が起きている部分時系列を検出する手法
- 周期性が崩れている部分等を見つける

![異常部位の例](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2269973/e9a6d85c-2f99-24b3-de0d-287ef74ead9c.png "異常部位の例")

### 変化点検知

- 時系列データのパターンが急激に変化する箇所を検知する手法
- 検索数やツイート数の急激な増加等を見つける

![変化点の例](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2269973/4f8d169b-c32f-ad77-a960-f07cb9936c3e.png "変化点の例")


今回のケースでは **変化点検知** が当てはまるかと思います。
サブスクでの再生数が急激に増加する箇所を「バズり点」とみなし、これを検知することを目標とします。

## 時系列予測モデル

変化点検知を行うにあたり、まず予測モデルを用意します。
信頼性のある予測モデルを用意し、予測値と実測値のズレが大きい場合に「変化があった」とみなすことで変化点を検知します。

今回は Facebook の提供している将来予測ライブラリ Prophet を用いて予測モデルを作成します。

[公式ドキュメント](https://facebook.github.io/prophet/)

こちらのモデルの優れた点は以下の 5 つです。([こちら](https://qiita.com/skillup_ai/items/7673d96dd95a33f91894)より引用)

> - 統計の知識がなくてもモデルを作れる
> - ドメイン知識を取り入れやすい
> - 特徴量エンジニアリングが必要ない
> - 欠損値があっても問題ない
> - 予測結果を解釈しやすい

イベント情報を与えることができるため、モデル作成時には、再生数のデータだけでなく、楽曲の配信開始日のデータも入力として与えます。
これにより、配信開始日から数日は再生数が上昇する といったようなトレンドを考慮したモデルになります。

## バズり検知ロジック

今回のバズり検知ロジックの流れは以下のとおりです。

1. 予測モデルを作る
2. 予測値と実測値から異常度を計算する
3. 異常度がしきい値以上であれば、その時点をバズり点とする

2 の異常度には、MPE の考え方を用います。（MPE については[こちら](https://analysis-navi.com/?p=2875)を参照）

理由としては、上昇のみ検知したいという点と、人気度のフェーズによってバズりと呼べる上昇幅が異なるのを均等化したいという点です。
（普段 10 再生のアーティストが 110 再生されるのと、普段 1000 再生のアーティストが 1100 再生されるのは、同じ +100 でも意味合いが異なる）
MPE は予測値の正確性を測る指標なので実測値が分母ですが、今回は予測値を正とみなした上で、実測値のズレを測りたいので、予測値を分母とします。

式は以下のとおりです。

```math
異常度 = \frac{実測値 - 予測値}{予測値}\times100
```

3 のしきい値に関しては今後検証を行っていく中でチューニングするものなので、今回は異常度をプロットし、どのくらいが妥当なのかを把握するにとどめます。

## 実装&amp;検証

弊社サブスクサービスにおける、某アーティストの再生数データを使用します。
（※再生数は加工を行っているため、実際の値とは異なります）

検証期間は 2018 年 1 月〜12 月で粒度は日次です。

各時点で 過去 100 日間の実測値を用いてモデルを作成し、最新の値を予測します。
（例: 2017/9/23〜2017/12/31 の実測値から 2018/1/1 の予測値を算出）

検証期間の異常度の推移をプロットし、実際の再生数や事実（その時本当にバズっていたか）と照らし合わせて考察を行います。

### 環境情報

- Mac OS Big Sur 11.6
- Python 3.8.6

使用するデータの中身を確認します。

```python
import pandas as pd

fname = "tsad.csv" # 再生数データ
fname_release = "tsad_release.csv" ## 配信開始日データ

df = pd.read_csv(fname, index_col="day_date", parse_dates=True)
df.head()
```

日次の再生数が入っています。

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }

</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>play_cnt</th>
    </tr>
    <tr>
      <th>day_date</th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2015-10-06</th>
      <td>4.0</td>
    </tr>
    <tr>
      <th>2015-10-07</th>
      <td>0.0</td>
    </tr>
    <tr>
      <th>2015-10-08</th>
      <td>0.0</td>
    </tr>
    <tr>
      <th>2015-10-09</th>
      <td>0.0</td>
    </tr>
    <tr>
      <th>2015-10-10</th>
      <td>0.0</td>
    </tr>
  </tbody>
</table>
</div>

```python
df_release = pd.read_csv(fname_release, index_col="open_date", parse_dates=True)
df_release.head()
```

楽曲の配信開始日情報が入っています。
イベントに重みはつけられないので、実際に必要なのは日付の要素のみです。（この楽曲の影響度を大きく とかはできない）

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }

</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>music_title_name</th>
    </tr>
    <tr>
      <th>open_date</th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2015-08-05</th>
      <td>楽曲名</td>
    </tr>
    <tr>
      <th>2015-08-05</th>
      <td>楽曲名</td>
    </tr>
    <tr>
      <th>2015-08-05</th>
      <td>楽曲名</td>
    </tr>
    <tr>
      <th>2015-08-05</th>
      <td>楽曲名</td>
    </tr>
    <tr>
      <th>2015-08-05</th>
      <td>楽曲名</td>
    </tr>
  </tbody>
</table>
</div>

```python
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = [16, 9]
plt.rcParams["font.size"] = 20
plt.tight_layout()

df["2018"].plot(y="play_cnt")
```

検証期間の再生数推移を可視化します。
2 月頃からじわじわ伸びて、7 月に大きく上昇、10 月末に更に上昇という感じに見えます。

![play_cnt](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2269973/675086ba-9f89-e529-4db5-e2785a9c2c2a.png "検証期間の再生数推移")


検証期間の予測値および異常度を取得します。

```python
import pandas as pd
from fbprophet import Prophet
from statistics import mean
from datetime import datetime as dt

fname = "tsad.csv"
fname_release = "tsad_release.csv"

_df = pd.read_csv(fname, index_col="day_date", parse_dates=True)

df_release = pd.read_csv(fname_release)

# Prophetの形式に変換
holidays = pd.DataFrame(
    {
        "holiday": "release",
        "ds": pd.to_datetime(df_release["open_date"]),
        "lower_window": 0,
        "upper_window": 7,  # 配信開始後7日間再生数に影響があるとする
    }
)

strdt = dt.strptime("2018-01-01", "%Y-%m-%d")  # 検証開始日
enddt = dt.strptime("2018-12-31", "%Y-%m-%d")  # 検証終了日

datelist = pd.date_range(start=strdt, end=enddt, freq="D")

test_length = 1  # 各時点で最新の値のみ取得

predict = pd.DataFrame([])

for date in datelist:
    # Prophetの形式に変換
    df = _df[:date].reset_index()
    df = df[["day_date", "play_cnt"]].rename(
        columns={"day_date": "ds", "play_cnt": "y"}
    )

    df["floor"] = 0  # 予測値の下限を0に設定
    df_train = df.iloc[-(test_length + 100) : -test_length]  # 過去100日のデータから予測
    df_test = df.iloc[-test_length:]

    prophet_model = Prophet(seasonality_mode="multiplicative", holidays=holidays)
    prophet_model.fit(df_train)

    future = prophet_model.make_future_dataframe(periods=test_length, freq="D")
    pred = prophet_model.predict(future)

    df_test["Prophet Predict"] = pred.iloc[-test_length:]["yhat"].values
    df_test["Prophet Predict Lower"] = pred.iloc[-test_length:]["yhat_lower"].values
    df_test["Prophet Predict Upper"] = pred.iloc[-test_length:]["yhat_upper"].values

    df_test["Anomaly"] = (
        mean((df_test["y"] - df_test["Prophet Predict"]) / df_test["Prophet Predict"])
        * 100
    )

    predict = predict.append(df_test)

predict.to_csv(f"predict_{strdt}_to_{enddt}_{fname}", index=False)

```

実測値と予測値をプロットします。

```python
df = pd.read_csv("predict_2018-01-01 00:00:00_to_2018-12-31 00:00:00_tsad.csv", index_col='ds', parse_dates=True)
df[["y", "Prophet Predict"]].plot()
```

大局的な動きは一致しているかと思います。
前日までのデータ使って予測しているので、大きな動きがあった場合もすぐに軌道修正できているようです。
ただ、後半はかなりズレが目立ちます。

![predict](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2269973/650f412f-b7af-8502-9324-65c8de1f5102.png "実測値と予測値")


異常度をプロットします。

```python
df[["Anomaly"]].plot()
```

予測値を上回った場合は正、下回った場合は負の値となります。
大きく上振れている点をバズり点とみなします。
かなり大きく上振れるので、しきい値は 200 くらいが妥当なのかなと思います。

1~2 月にバズり点が集中しており、その後、7 月に 1 点、11 月に 1 点といったとこでしょうか。
次項で下図の3箇所に焦点を当てて、考察を行います。

![anomaly](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2269973/25763d9c-67df-6539-000c-4852044c42c6.png "異常度")


## 考察

1 月後半に異常度が 400 を超えた点がありますが、この日は某アーティストが TV 番組で大きく取り上げられて話題になった日でした。（[1]の最大値）
再生数の推移だけでは前後と比較して突出しているようには見えないので、**異常検知のロジックだからこそ検知できた点**と言えるかと思います。

7 月後半の突出した点は新曲の配信開始日と重なっていました。[2]
対象の曲はそれまでと比較して弊社でも格段に売れ行きがよかったので、配信開始日の影響を考慮した上でも、予測値を大きく上回ったのだと思われます。
この曲のヒットによって再度バズり、人気がもう 1 段階上昇したと解釈できるかと思います。

11 月の点に関しては、イベントは確認できませんでしたが、10 月後半に新曲が配信されており、その影響か再生数が急激に伸びています。[3]
年末に向けて何らかの特集がされていたのかもしれません。

一方で、1~2 月は全体的に異常度が高い点が頻出しており、この原因が気になります。[1]

1~2 月部分の予測値と実測値を拡大すると以下のとおりです。
1 月は実測値と予測値で上昇と下降が逆転している箇所が多く、予測値が低いほど異常度が高くなるため、バズり点が多くなったようです。
2 月は実測値が明確にスパイクしており、頻繁にバズりを繰り返していることが見て取れます。

このように予測値が低い時点では、バズり点が生まれやすいこと、質の異なるバズり点が混在することがわかります。

![predict_12](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2269973/ee40ac9a-a0df-767a-c429-e4607da71c7f.png "1~2月の実測値と予測値")


## 今後の課題

バズり検知を自動化するためには、多くのアーティストのデータで試し、**汎用的なしきい値を決める**というのが 1 番の課題です。
バズり点が定義されたデータで検証を行い、バズり点に対応する時点でどの程度の異常度となるかを確認しながらチューニングしていこうと思います。

また、上記はモデルの予測精度が保証されていることが前提なので、予測モデルのチューニングもしていく必要があります。
具体的には下記が考えられます。

- 学習データ期間のチューニング（今回は過去 100 日分で学習）
- イベント影響期間のチューニング (今回は 7 日間で設定)
- イベント情報追加
  - TV やライブの出演情報を取得してモデルに組み込む

その他、しきい値を 5 日連続して超えたらバズりとする等、バズりの定義を見直すことも考えられます。

## むすび

今回はサブスクの再生数推移データに異常検知ロジックを適用することでバズり検知の自動化を試みました。

バズりという言葉自体が曖昧ということもあり、はっきりと「バズり検知できた」と言える結果かは怪しいですが、バズり点がしっかり定義された検証データを用意してチューニングを行っていけば、それなりに使えるものができるかなと思いました。

ここまで読んでいただきありがとうございました！

明日の[レコチョク Advent Calendar 2021](https://qiita.com/advent-calendar/2021/recochoku) は 14 日目 「Web3.0 関連の情報収集に活用しているサイト 5 選」です。お楽しみに！

## 参考文献

- https://www.albert2005.co.jp/knowledge/machine_learning/anomaly_detection_basics/anomaly_detection_time
- https://qiita.com/skillup_ai/items/7673d96dd95a33f91894
- https://analysis-navi.com/?p=2875
- https://www.salesanalytics.co.jp/datascience/datascience006/

---
この記事はレコチョクのエンジニアブログの記事を転載したものとなります。

https://techblog.recochoku.jp/7925