2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

ベイズ機械学習で解釈可能なレコメンドエンジンを作ってみた

Last updated at Posted at 2025-08-09

はじめに

こんにちは、事業会社で働いているデータサイエンティストです。

これは、レコメンドエンジンの推定結果を用いて、限られた資源をいかにユーザーに最適配分するかを説明する記事の第1回です。
今回の記事では、MovieLens のデータをもとにした解釈可能なベイズ機械学習モデルについて解説し、次回の記事で実際の資源配分手法を紹介します。

本モデルには、以下のような特徴があります。

  • 未評価の理由も推論
    ユーザーが映画を評価しなかった場合、それが「見たけど好きじゃなかった」のか、「そもそも見ていない」のかを推論
  • 次元数を自動推定
    モデルは無限次元を仮定し、その中からデータに基づいて実際に必要な次元数を決定
  • 事前分布もデータ駆動
    ユーザーと映画の埋め込みベクトルの事前分布は、固定的な正規分布ではなく、データから推定される柔軟な構造を採用

分析の結果、以下のような傾向が得られました。

  • MovieLensの評価構造は、1次元(商業志向のアメリカ映画 ↔ アート志向の映画) でほぼ説明可能。
  • ユーザーの嗜好分布は、釣鐘型(正規分布に近い形) となっている。

あらかじめお詫びしますが、私は映画にほとんど興味がなく、人生で映画館に行ったのはたった2回、しかも両方とも同じ『名探偵コナン』の作品でした笑。そのため、モデル自体は私が考案したのですが、映画に関する解釈はほぼウィキペディアと生成AIに頼っています。しかし逆に言えば、ドメイン知識がほぼゼロでも、データから傾向を自動発見できるのが、このモデル構造の強みだと考えています。

モデル構造

さて、ここからはモデルの構造について詳しく説明していきます。まずは尤度関数、つまり「ユーザー × 映画」の評価ペアがどのように生成されるかを解説し、その後で詳細な事前分布について触れていきます。

ユーザー x 映画の評価

まず直感的にイメージしてみましょう。

MovieLensだけではないですが、私たちがユーザーからの評価データを見ると、「評価がない」(欠損)ユーザー・アイテム(本記事だと映画)ペアのほとんどです。ただ、それが本当にその映画が嫌いだったのか、それとも単に見ていなかっただけなのかは分かりません。例えば、ある映画を見ていなければ評価が付かないのは当然ですが、見た上で評価しなかった場合は「好きではない」という意思表示かもしれません。

このように「評価がない」理由は複数あるため、単純に「評価がない=嫌い」と決めつけるのは危険です。

そこで本モデルでは、評価データの生成過程を2段階に分けて考えます。

  1. まずはその映画を見たか(視聴・検索)?
  2. 見た上で、どのくらい好きか(評価)はどうか?

この2段階に分けることで、たとえ評価が付いていなくても「未視聴である可能性がある」と推論でき、ユーザーの本当の好みをより正確に捉えることができます。ただし、実際に映画を見たかどうかはデータから直接は分からないため、潜在確率変数としてモデル内で推論します。

ここからは少し数式を使って、モデルの仕組みを説明します。

1. 視聴・検索段階(ロジスティック回帰)

ここでは $i$ をユーザーのID、$j$ を映画のIDとします。

ユーザーが映画を評価しない理由は大きく2つ考えられます。

  • そもそも映画を見ていない(未視聴)
  • 見たが評価しなかった(嫌い、または興味がない)

この違いを区別するために、まずユーザーがその映画を視聴・検索するかどうかをベルヌーイ分布でモデル化しています。この段階の(観測できない)視聴確率は、映画の基本的な認知度($\kappa$)とユーザーごとの探索傾向($\lambda$)、および共通の切片項($\mu$)の組み合わせで決定されます。

$$
P(視聴_{ij} = 1) = \mathrm{logit}^{-1}(\mu + \lambda_{i} + \kappa_{j})
$$

2. 評価段階(順序ロジスティック回帰)

視聴された映画に対して、評価なしを 0 として扱うと、ユーザーは 0〜5 段階などの評価を与えます。この評価は、映画の人気度($\phi$)と、ユーザーと映画の潜在ベクトルの内積に基づく好みの一致度によって決まります。

具体的には、ユーザーの潜在ベクトル ($\gamma$)と映画の潜在ベクトル($\zeta$)の次元ごとに重み付けされた($\delta$)内積が評価の中心的な指標です。

$$
評価スコア_{ij} = \phi_{j} + \sum_{d=1}^{\infty}\left( \gamma_{i,d} \cdot \zeta_{j,d} \cdot \delta_{d} \right)
$$

これは後述の事前分布の部分で詳しく説明しますが、$\sum_{d=1}^{\infty}\left( \gamma_{i,d} \cdot \zeta_{j,d} \cdot \delta_{d} \right)$ は発散せず、有限の値に収束することが数学的に証明されています。簡潔にいうと、$\delta$ は多くの次元に実質ゼロの重みをつけることで、モデルから不要な要素を自動・データドリブンに削除し、次元・変数選択を実現することができます。

では、視聴確率と評価スコアをもとに、ユーザーの評価データがどのように生成されるかを示します。

まず、潜在的な視聴フラグをサンプリングします:

$$
視聴_{ij} \sim P(視聴_{ij})
$$

もし $視聴_{ij} = 0$ の場合、視聴していないため、評価は欠損(0)となります:

$$
評価_{ij} = 0
$$

一方で $視聴_{ij} = 1$ の場合、視聴済みなので、次に順序ロジスティック分布から評価をサンプリングします:

$$
評価_{ij} \sim OrderedLogistic(評価スコア_{ij}, cutpoints_{i})
$$

簡潔さのため、$評価_{ij} \sim OrderedLogistic(評価スコア_{ij}, cutpoints_{i})$ を $P(評価_{ij})$ と書くことがあります。

評価結果は順序ロジスティックモデルによって生成され、ユーザーごとに異なるカットポイント($cutpoints$)を持つことで、評価尺度の違いや個人差を柔軟に捉えています。

生成モデルの概念としては以上ですが、残念ながらStanは $視聴_{ij}$ のような離散確率変数を扱うことができません。よって、$視聴_{ij}$ を周辺化する必要があります。

まずは、全てのパターンを列挙してみましょう:

  • $評価_{ij}$ = 0(視聴したかがわからない)
    • $視聴_{ij} = 0$: $P(視聴_{ij} = 0)$
    • $視聴_{ij} = 1$: $P(視聴_{ij} = 1) \cdot P(評価_{ij}=0)$
  • $評価_{ij}$ = 1(視聴した)
    • $視聴_{ij} = 1$: $P(視聴_{ij} = 1) \cdot P(評価_{ij}=1)$
  • $評価_{ij}$ = 2(視聴した)
    • $視聴_{ij} = 1$: $P(視聴_{ij} = 1) \cdot P(評価_{ij}=2)$
  • $評価_{ij}$ = 3(視聴した)
    • $視聴_{ij} = 1$: $P(視聴_{ij} = 1) \cdot P(評価_{ij}=3)$
  • $評価_{ij}$ = 4(視聴した)
    • $視聴_{ij} = 1$: $P(視聴_{ij} = 1) \cdot P(評価_{ij}=4)$
  • $評価_{ij}$ = 5(視聴した)
    • $視聴_{ij} = 1$: $P(視聴_{ij} = 1) \cdot P(評価_{ij}=5)$

よって、視聴したかどうかがわからないのは、$評価_{ij}$ = 0 の時だけです。評価が1 ~ 5の場合は、Stanでいうと bernoulli_logit_lpmf(1 | ...) + ordered_logistic_lpmf(1 ~ 5 | ...) のように、それぞれの部分の対数尤度を足せばいいです。

問題なのは $評価_{ij}$ = 0の場合ですが、対処方法は実は簡単で、まずそれぞれのパターンの対数尤度を計算して、log_sum_exp 関数に入れれば計算できます。詳細はStanのコードをご確認ください。

この2段階モデルにより、評価が存在しない理由の切り分けと、評価の質的な差異を同時に推定可能になっています。この仕組みが、単なる「好き嫌い」では説明できないデータの複雑な構造を捉える鍵となっています。

事前分布

ここの内容は、基本的には前回の記事の内容をそのまま利用しています。

ユーザー埋め込みベクトル

ユーザー埋め込みベクトルの事前分布には、正規分布を設定することが多いですが、本当に正規分布で良いのか?という疑問に回答するため、正規分布から成すディリクレ過程で設定しています。

まずは全体のハイパーパラメーター$\alpha_{user}$をガンマ分布からサンプリングします。

$$
\alpha_{user} \sim Gamma(0.001, 0.001)
$$

次に、棒折り過程でディリクレ過程を構築します。無限にあるクラスターの中のクラスターpについてこのように必要な変数をサンプリングします:

$$
\pi_{user,p}\sim Beta(1, \alpha_{user})
$$

$$
p_{user,p} = \pi_{user,p} \prod\limits_{l=1}^{p - 1} (1 - \pi_{user,l})
$$

また、クラスター $ p_{user} $ の中心ベクトルとそのばらつきは下記のようにサンプリングされます:

$$
P_{latent, user, p_{user}} \sim Normal(0,1)
$$

$$
P_{\sigma,user, p_{user}} \sim Gamma(0.001, 0.001)
$$

単語Sのベクトルをサンプリングする際は、まず上記の棒折り過程で構築した無限次元の確率ベクトルをパラメーターとするカテゴリ分布から、インデックスをサンプリングします:

$$
\eta_{i} \sim Categorical(p_{user})
$$

次に、サンプリングされた$\eta_{i}$を基に、ベクトルをサンプリングします:

$$
\gamma_{i} \sim Normal(P_{latent, user, \eta_{i}}, P_{\sigma, user, \eta_{i}})
$$

この一連のサンプリングのプロセス分布を$G_{user}$で表します。

映画埋め込みベクトル

映画の埋め込みベクトルの構造はユーザー側と一緒です。

まずは全体のハイパーパラメーター$\alpha_{movie}$をガンマ分布からサンプリングします。

$$
\alpha_{movie} \sim Gamma(0.001, 0.001)
$$

次に、棒折り過程でディリクレ過程を構築します。無限にあるクラスターの中のクラスターpについてこのように必要な変数をサンプリングします:

$$
\pi_{movie,p}\sim Beta(1, \alpha_{movie})
$$

$$
p_{movie,p} = \pi_{movie,p} \prod\limits_{l=1}^{p - 1} (1 - \pi_{movie,l})
$$

また、クラスター $ p_{movie} $ の中心ベクトルとそのばらつきは下記のようにサンプリングされます:

$$
P_{latent, novie, p_{movie}} \sim Normal(0,1)
$$

$$
P_{\sigma, movie, p_{movie}} \sim Gamma(0.001, 0.001)
$$

単語Sのベクトルをサンプリングする際は、まず上記の棒折り過程で構築した無限次元の確率ベクトルをパラメーターとするカテゴリ分布から、インデックスをサンプリングします:

$$
\eta_{j} \sim Categorical(p_{movie})
$$

次に、サンプリングされた$\eta_{i}$を基に、ベクトルをサンプリングします:

$$
\zeta_{j} \sim Normal(P_{latent, movie, \eta_{j}}, P_{\sigma, movie, \eta_{j}})
$$

この一連のサンプリングのプロセス分布を$G_{movie}$で表します。

次元重要度

次元数は棒折り過程で定式化します。具体的には、

$$
\beta \sim Gamma(0.001, 0.001)
$$

をサンプリングし、次に棒折り過程で $\delta_{d}$ を求めます:

$$
\omega_{d}\sim Beta(1, \beta)
$$

$$
\delta_{d} = \omega_{d} \prod\limits_{l=1}^{d - 1} (1 - \omega_{l})
$$

これをd = 1からd = $\infty$まで繰り返します。

さて、詳細は筆者の論文のAppendix(Wu, 2025)を参照していただきたいですが:

次元重要度のベクトルのd番目の次元の期待値はこのように書くことができます:

$$
\mathbb{E}[\delta_{d}] = \frac{1}{1 + \beta} \left( \frac{\beta}{1 + \beta} \right)^{d-1}.
$$

$\beta$がガンマ分布から生成されているため、0より大きい実数で、$\frac{\beta}{1+\beta}$は0から1までの間の実数になります。したがって、次元重要度のベクトルは、次元を追うごとに指数関数的にゼロへ収束していく特性を持っていることがわかります。

このように、ノンパラメトリックベイズの代表的な構造である棒折り過程を用いることで、無限次元の重要度ベクトルであっても収束性とモデルの数学的健全性が保たれているだけでなく、急速に値がゼロに収束する特性は次元の変数選択の機能も果たすこともできます。

Stanコード

Stanによるモデル実装のコードはこちらです。他のあまり重要ではないパラメータの事前分布の設定は実装コードよりご確認いただければと思います:

recommend.stan
functions {
  vector stick_breaking(vector breaks){
    int length = size(breaks) + 1;
    vector[length] result;
    
    result[1] = breaks[1];
    real summed = result[1];
    for (d in 2:(length - 1)) {
      result[d] = (1 - summed) * breaks[d];
      summed += result[d];
    }
    result[length] = 1 - summed;
    
    return result;
  }
  
  real partial_sum_lpmf(
    array[] int result,
    
    int start, int end,
    
    array[] int user, array[] int movie,
    
    real search_intercept,
    vector search_user, vector search_movie,
    
    array[] vector cutpoints,
    vector dimension,
    vector movie_propensity, 
    array[] vector user_latent, array[] vector movie_latent
  ){
    vector[end - start + 1] lambda;
    int count = 1;
    for (i in start:end){
      if (result[count] == 1) {
        vector[2] case_when;
        case_when[1] = bernoulli_logit_lpmf(0 | search_intercept + search_user[user[i]] + search_movie[movie[i]]);
        case_when[2] = bernoulli_logit_lpmf(1 | search_intercept + search_user[user[i]] + search_movie[movie[i]]) + 
                       ordered_logistic_lpmf(
                         result[count] | 
                         movie_propensity[movie[i]] + 
                         user_latent[user[i]] '* (movie_latent[movie[i]] .* dimension),
                         cutpoints[user[i]]
                       );
        lambda[count] = log_sum_exp(case_when); 
      } else {
        lambda[count] = bernoulli_logit_lpmf(1 | search_intercept + search_user[user[i]] + search_movie[movie[i]]) +
                        ordered_logistic_lpmf(
                          result[count] | 
                          movie_propensity[movie[i]] + 
                          user_latent[user[i]] '* (movie_latent[movie[i]] .* dimension),
                          cutpoints[user[i]]
                        );
      }
      count += 1;
    }
    return sum(lambda);
  }
}
data {
  int user_type;
  int movie_type;
  int result_type;
  int dimension_type;
  int group_type;
  
  int N;
  array[N] int user;
  array[N] int movie;
  array[N] int result;
  
  int val_N;
  array[val_N] int val_user;
  array[val_N] int val_movie;
  array[val_N] int val_result;
}
parameters {
  real<lower=0> dimension_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[dimension_type - 1] dimension_breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  real<lower=0> group_user_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[group_type - 1] group_user_breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  real<lower=0> group_movie_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[group_type - 1] group_movie_breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  vector<lower=0>[group_type] group_user_sigma;
  array[group_type] vector[dimension_type] group_user_latent;
  
  vector<lower=0>[group_type] group_movie_sigma;
  array[group_type] vector[dimension_type] group_movie_latent;
  
  vector[movie_type] movie_propensity_unnormalized;
  array[user_type] vector[dimension_type] user_latent;
  array[movie_type] vector[dimension_type] movie_latent;
  
  array[user_type] ordered[result_type - 1] cutpoints;
  
  real search_intercept;
  vector[user_type] search_user_unnormalized;
  vector[movie_type] search_movie_unnormalized;
}
transformed parameters {
  simplex[dimension_type] dimension;
  simplex[group_type] group_user;
  simplex[group_type] group_movie;
  
  dimension = stick_breaking(dimension_breaks);
  group_user = stick_breaking(group_user_breaks);
  group_movie = stick_breaking(group_movie_breaks);
  
  vector[movie_type] movie_propensity = movie_propensity_unnormalized - mean(movie_propensity_unnormalized);
  vector[user_type] search_user = search_user_unnormalized - mean(search_user_unnormalized);
  vector[movie_type] search_movie = search_movie_unnormalized - mean(search_movie_unnormalized);
}
model {
  dimension_alpha ~ gamma(0.001, 0.001);
  dimension_breaks ~ beta(1, dimension_alpha);
  
  group_user_alpha ~ gamma(0.001, 0.001);
  group_user_breaks ~ beta(1, group_user_alpha);
  
  group_movie_alpha ~ gamma(0.001, 0.001);
  group_movie_breaks ~ beta(1, group_movie_alpha);
  
  movie_propensity ~ normal(0, 10);
  
  group_user_sigma ~ gamma(0.001, 0.001);
  group_movie_sigma ~ gamma(0.001, 0.001);
  
  for (i in 1:group_type){
    group_user_latent[i] ~ normal(0, 1);
    group_movie_latent[i] ~ normal(0, 1);
  }
  
  for (i in 1:user_type){
    vector[group_type] case_vector;
    for (j in 1:group_type){
      case_vector[j] = log(group_user[j]) + normal_lpdf(user_latent[i] | group_user_latent[j], group_user_sigma[j]);
    }
    target += log_sum_exp(case_vector);
  }
  
  for (i in 1:movie_type){
    vector[group_type] case_vector;
    for (j in 1:group_type){
      case_vector[j] = log(group_movie[j]) + normal_lpdf(movie_latent[i] | group_movie_latent[j], group_movie_sigma[j]);
    }
    target += log_sum_exp(case_vector);
  }
  
  search_intercept ~ normal(0, 10);
  search_user ~ normal(0, 10);
  search_movie ~ normal(0, 10);
  
  int grainsize = 1;
  
  target += reduce_sum(
    partial_sum_lupmf, result,
    
    grainsize,
    
    user, movie,
    
    search_intercept,
    search_user, search_movie,
    
    cutpoints,
    dimension,
    movie_propensity, 
    user_latent, movie_latent
  );
}
generated quantities {
  array[user_type] vector[dimension_type] weighted_user_latent;
  array[movie_type] vector[dimension_type] weighted_movie_latent;
  array[dimension_type] real drawn_G_user;
  array[dimension_type] real drawn_G_movie;
  array[val_N] real predicted_first_score;
  array[val_N] real predicted_second_score;
  array[val_N] int predicted_result;
  vector[result_type] F1_score_by_class;
  real F1_score_macro;
  
  for (i in 1:user_type){
    weighted_user_latent[i] = user_latent[i] .* dimension;
  }
  
  for (i in 1:movie_type){
    weighted_movie_latent[i] = movie_latent[i] .* dimension;
  }
  
  {
    int sampled_group = categorical_rng(group_user);
    drawn_G_user = normal_rng(group_user_latent[sampled_group], group_user_sigma[sampled_group]);
  }
  {
    int sampled_group = categorical_rng(group_movie);
    drawn_G_movie = normal_rng(group_movie_latent[sampled_group], group_movie_sigma[sampled_group]);
  }
  
  for (i in 1:val_N){
    predicted_first_score[i] = search_intercept + search_user[val_user[i]] + search_movie[val_movie[i]];
    predicted_second_score[i] = movie_propensity[val_movie[i]] + user_latent[val_user[i]] '* (movie_latent[val_movie[i]] .* dimension);
    int coin = bernoulli_logit_rng(predicted_first_score[i]);
    if (coin == 0){
      predicted_result[i] = 1;
    } else {
      predicted_result[i] = ordered_logistic_rng(predicted_second_score[i], cutpoints[val_user[i]]);
    }
  }
  
  for (r in 1:result_type){
    int TP = 0;
    int FP = 0;
    int FN = 0;
    for (i in 1:val_N){
      if (val_result[i] == r && predicted_result[i] == r){
        TP += 1;
      }
      else if (val_result[i] != r && predicted_result[i] == r){
        FP += 1;
      }
      else if (val_result[i] == r && predicted_result[i] != r){
        FN += 1;
      }
    }
    F1_score_by_class[r] = (2 * TP) * 1.0/(2 * TP + FP + FN);
  }
  
  F1_score_macro = mean(F1_score_by_class);
}

前処理

ここからは前処理に入ります。まずデータを取得して、元の行列形式の横持ちから縦持ちに変換します:

data("MovieLense", package = "recommenderlab")

movie_df <- as(MovieLense, "matrix") |>
  as.data.frame() |> 
  tibble::tibble() |>
  dplyr::mutate(
    user_id = dplyr::row_number()
  ) |>
  tidyr::pivot_longer(!user_id, names_to = "movie", values_to = "point") |>
  tidyr::replace_na(
    list(
      point = 0
    )
  )

次に、映画マスターを作成して元のテーブルに LEFT JOIN して、映画 ID を付与します:

movie_master <- movie_df |>
  dplyr::select(movie) |> 
  dplyr::distinct() |>
  dplyr::arrange(movie) |>
  dplyr::mutate(
    movie_id = dplyr::row_number()
  )

movie_df_with_id <- movie_df |>
  dplyr::left_join(
    movie_master, by = "movie"
  )

次のステップに入る前に、まず評価の点数の分布をご確認ください:

> movie_df_with_id$point |> table()

      0       1       2       3       4       5 
1469760    6059   11307   27002   33947   21077 

大半というか、ほぼ全てが 0 (評価なし、欠損)に埋め尽くされています。欠損は情報としての重要性は低いと思われるため、ダウンサンプリングを実施します:

set.seed(12345)
movie_df_with_id_downsampled <- movie_df_with_id |>
  dplyr::filter(point > 0) |>
  dplyr::bind_rows(
    movie_df_with_id |>
      dplyr::filter(point == 0) |>
      dplyr::slice_sample(n = sum(movie_df_with_id$point > 0))
  )

次に、検証用の行 ID を残します:

set.seed(12345)
test_id <- c(
  sample(which(movie_df_with_id_downsampled$point == 0), 250),
  sample(which(movie_df_with_id_downsampled$point == 1), 250),
  sample(which(movie_df_with_id_downsampled$point == 2), 250),
  sample(which(movie_df_with_id_downsampled$point == 3), 250),
  sample(which(movie_df_with_id_downsampled$point == 4), 250),
  sample(which(movie_df_with_id_downsampled$point == 5), 250)
)

学習データの量を確認しましょう:

> nrow(movie_df_with_id_downsampled[-test_id,])
[1] 197284

ダウンサンプリングと検証データの除外を行った結果、約20万件のデータがあります。

最後に、モデルをコンパイルして、

m_init <- cmdstanr::cmdstan_model("recommend.stan",
                                  cpp_options = list(
                                    stan_threads = TRUE
                                  )
                                  )

8 つのコアを利用した変分推論でモデル推定を実施します:

> m_estimate <- m_init$variational(
     seed = 12345,
     threads = 8,
     data = list(
         user_type = max(movie_df_with_id$user_id),
         movie_type = nrow(movie_master),
         result_type = max(movie_df_with_id$point) + 1,
         dimension_type = 20,
         group_type = 10,
         
         N = nrow(movie_df_with_id[-test_id,]),
         user = movie_df_with_id$user_id[-test_id],
         movie = movie_df_with_id$movie_id[-test_id],
         result = movie_df_with_id$point[-test_id] + 1,
         
         val_N = nrow(movie_df_with_id[test_id,]),
         val_user = movie_df_with_id$user_id[test_id],
         val_movie = movie_df_with_id$movie_id[test_id],
         val_result = movie_df_with_id$point[test_id] + 1
     )
 )
------------------------------------------------------------ 
EXPERIMENTAL ALGORITHM: 
  This procedure has not been thoroughly tested and may be unstable 
  or buggy. The interface is subject to change. 
------------------------------------------------------------ 
Gradient evaluation took 0.118352 seconds 
1000 transitions using 10 leapfrog steps per transition would take 1183.52 seconds. 
Adjust your expectations accordingly! 
Begin eta adaptation. 
Iteration:   1 / 250 [  0%]  (Adaptation) 
Iteration:  50 / 250 [ 20%]  (Adaptation) 
Iteration: 100 / 250 [ 40%]  (Adaptation) 
Iteration: 150 / 250 [ 60%]  (Adaptation) 
Iteration: 200 / 250 [ 80%]  (Adaptation) 
Iteration: 250 / 250 [100%]  (Adaptation) 
Success! Found best value [eta = 0.1]. 
Begin stochastic gradient ascent. 
  iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
   100      -314416.202             1.000            1.000 
   200      -256408.094             0.613            1.000 
   300      -237280.195             0.436            0.226 
   400      -227671.368             0.337            0.226 
   500      -221809.193             0.275            0.081 
   600      -217867.794             0.232            0.081 
   700      -215126.504             0.201            0.042 
   800      -213123.680             0.177            0.042 
   900      -211455.421             0.158            0.026 
  1000      -210166.275             0.143            0.026 
  1100      -209136.137             0.043            0.018 
  1200      -208228.375             0.021            0.013 
  1300      -207479.388             0.014            0.009   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  285.8 seconds.

285 秒で終わりました!

推論結果をデータフレイム形式で保存します:

m_summary <- m_estimate$summary()

推定結果

では、早速モデル推定の結果を確認しましょう!

予測精度

まず、F1スコアの事後分布から確認しましょう:

m_estimate$draws("F1_score_by_class") |>
  tibble::as_tibble() |>
  `colnames<-`(as.character(c(0:5))) |>
  dplyr::mutate(
    rid = dplyr::row_number(),
    dplyr::across(
      dplyr::everything(),
      ~ as.numeric(.)
    )
  ) |>
  tidyr::pivot_longer(!rid, names_to = "type") |>
  dplyr::bind_rows(
    m_estimate$draws("F1_score_macro") |>
      tibble::as_tibble() |>
      `colnames<-`("macro") |>
      dplyr::mutate(
        rid = dplyr::row_number(),
        dplyr::across(
          dplyr::everything(),
          ~ as.numeric(.)
        )
      ) |>
      tidyr::pivot_longer(!rid, names_to = "type")
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_density(ggplot2::aes(x = value, fill = type), alpha = 0.5) + 
  ggplot2::labs(
    title = "posterior distributions of F1 scores"
  )

f1_scores.png

正直、F1スコアから見ると、非常に悪い結果になっています。ただ、本当にそうなのか?他の観点からの予測結果の可視化も実施しましょう。まずは、評価と評価スコアの関係です:

m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^predicted_second_score\\[")) |>
  dplyr::bind_cols(
    ans = movie_df_with_id_downsampled$point[test_id]
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_violin(ggplot2::aes(x = as.factor(ans), y = mean), trim = FALSE, fill = "#1E90FF", alpha = 0.7) +
  ggplot2::geom_jitter(ggplot2::aes(x = as.factor(ans), y = mean), width = 0.2, alpha = 0.3, color = "black") + 
  ggplot2::labs(
    x = "評価",
    y = "評価スコア"
  ) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

second_score.png

評価が高ければ高いほど、評価スコアも高くなる傾向があります。

では、純粋に評価と評価スコアの回帰モデルを組んだらどうなるのかも確認しましょう:

> m_summary |>
     dplyr::filter(stringr::str_detect(variable, "^predicted_second_score\\[")) |>
     dplyr::bind_cols(
         ans = movie_df_with_id_downsampled$point[test_id]
     ) |>
     lm(ans ~ mean, data = _) |>
     summary()

Call:
lm(formula = ans ~ mean, data = dplyr::bind_cols(dplyr::filter(m_summary, 
    stringr::str_detect(variable, "^predicted_second_score\\[")), 
    ans = movie_df_with_id_downsampled$point[test_id]))

Residuals:
    Min      1Q  Median      3Q     Max 
-5.1035 -0.9351  0.0383  1.0249  4.5176 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.81831    0.04293   42.35   <2e-16 ***
mean         0.55611    0.01960   28.38   <2e-16 ***
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

Residual standard error: 1.378 on 1498 degrees of freedom
Multiple R-squared:  0.3496,	Adjusted R-squared:  0.3492 
F-statistic: 805.3 on 1 and 1498 DF,  p-value: < 2.2e-16

係数が 0.55611 で、決定係数も 0.3492 になっています。よって、評価スコアと評価には一定程度の正の関係があるといえるでしょう。

最後に、評価と評価スコアのペアを評価スコアの昇順で並べて、その際の評価の分布状況を確認しましょう:

m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^predicted_second_score\\[")) |>
  dplyr::bind_cols(
    ans = movie_df_with_id_downsampled$point[test_id]
  ) |>
  dplyr::arrange(mean) |>
  dplyr::mutate(
    rid = dplyr::row_number()
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = rid, y = mean, color = ans)) +
  ggplot2::geom_errorbar(ggplot2::aes(x = rid, ymin = q5, ymax = q95, color = ans), alpha = 0.3) +
  ggplot2::scale_color_gradient(
    low = "yellow",
    high = "blue",
    limits = c(0, 5)
  ) + 
  ggplot2::labs(
    x = "評価スコアの順位",
    y = "評価スコア",
    color = "評価"
  ) +
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

sort.png

ここではより明確に見えますが、完璧とは言えないものの、評価スコアが低いほどデータが黄色く(評価が低い)、評価スコアが高いほどデータが青くなる(評価が高い)傾向が視覚的に確認できます。

よって、完璧な精度ではないものの、モデルはある程度ユーザーの視聴選好を適切に推定できたと考えられます。

また、冒頭で紹介したF1スコアについてですが、これは元々順序関係を考慮しない多クラス分類向けの評価指標であり、本記事で扱うような順序データの分類精度評価には必ずしも最適ではありません。機械学習の分野では順序データを意識した評価指標の開発はまだ十分とは言えませんが、政治学の方法論研究に携わる身としては、こうした評価指標の改良・開発に期待を寄せています。

ユーザーと次元の全体像

ここではまず、次元の重要度を可視化しましょう:

m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^dimension\\[")) |>
  dplyr::mutate(
    id = dplyr::row_number()
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_bar(ggplot2::aes(x = as.factor(id), y = mean), 
                    stat = "identity", fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::geom_errorbar(
    ggplot2::aes(x = as.factor(id), ymin = q5, ymax = q95),
    width = 0.2
  ) + 
  ggplot2::labs(
    title = "次元重要度の事後分布",
    x = "次元ID",
    y = "重要度") +
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

dimension.png

おーい!!!!!!ほぼ一次元でデータを全部説明できるやんwwwwwww

少しつまらない(?)結果ですが、これはビジネス的にも非常に意義のある発見です。まず、この次元の正体を明らかにすることで、マーケティング戦略に有益な示唆を提供できます。

また、エンジニアリングの観点から見ると、数百あるいは数千次元のモデルではなく、実際には一次元で十分であれば、モデルのファイルサイズが大幅に軽くなり、実装が容易になるだけでなく、推論(予測値を出すこと)の速度も自然と向上します。

では、ユーザーと映画の第一次元の分布を確認しましょう:

g_hist_movie <- m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^movie_latent\\[")) |>
  dplyr::mutate(id = variable |> purrr::map(\(x){as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])})) |> 
  tidyr::unnest_wider(id, names_sep = "_") |> 
  dplyr::filter(id_2 == 1) |>
  ggplot2::ggplot() + 
  ggplot2::geom_histogram(ggplot2::aes(x = mean), fill = ggplot2::alpha("blue", 0.5)) + 
  ggplot2::labs(
    title = "movie first dimension"
  )

g_hist_user <- m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^user_latent\\[")) |>
  dplyr::mutate(id = variable |> purrr::map(\(x){as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])})) |> 
  tidyr::unnest_wider(id, names_sep = "_") |> 
  dplyr::filter(id_2 == 1) |>
  ggplot2::ggplot() + 
  ggplot2::geom_histogram(ggplot2::aes(x = mean), fill = ggplot2::alpha("blue", 0.5)) + 
  ggplot2::labs(
    title = "user first dimension"
  )

library(patchwork)

g_hist_movie | g_hist_user

user_movie_hist.png

ここで改めて強調したいのは、本モデルではユーザーと映画の潜在ベクトルの事前分布に、単純な正規分布のような釣鐘型の分布を仮定しているわけではなく、データから最適な事前分布を推定している点です。

したがって、今回の可視化でユーザーの第一次元の分布が単峰型(釣鐘型)になっているのは、事前分布の制約によるものではなく、純粋にデータから推論された結果と捉えて問題ありません。

ビジネス的に注目すべきは、この単峰型のユーザー嗜好の分布が示す意味です。政治分野でよく見られる二峰型(二極化)とは異なり、ユーザーの映画の好みが一つの中心に集約されているため、マーケティング戦略やリソース配分を検討する際には、中央寄せのアプローチが最適であると考えられます。つまり、極端に特定のジャンルや種類のコンテンツを嫌うユーザーは非常に少なく、多くのユーザーに受け入れられるコンテンツを重視する施策が効果的だと言えるでしょう。

映画のパラメータの解釈

映画認知度

まず、映画の基本的な認知度($\kappa$)が一番高い映画と低い映画をそれぞれ確認しましょう。

まず上位の映画は:

> m_summary |>
     dplyr::filter(stringr::str_detect(variable, "^search_movie\\[")) |>
     dplyr::bind_cols(
         m = movie_master$movie
     ) |>
     dplyr::arrange(-mean) |>
     dplyr::select(m, mean, q5, q95)
# A tibble: 1,664 × 4
   m                            mean    q5   q95
   <chr>                       <dbl> <dbl> <dbl>
 1 Scream (1996)                7.11  6.10  8.02
 2 Contact (1997)               6.86  5.94  7.80
 3 Dante's Peak (1997)          6.71  5.84  7.60
 4 Evita (1996)                 6.51  5.43  7.54
 5 Saint, The (1997)            6.35  5.44  7.19
 6 Liar Liar (1997)             6.27  5.64  6.90
 7 English Patient, The (1996)  6.07  5.49  6.74
 8 Game, The (1997)             5.90  5.18  6.63
 9 Murder at 1600 (1997)        5.72  4.96  6.48
10 Air Force One (1997)         5.67  5.09  6.22
# ℹ 1,654 more rows
# ℹ Use `print(n = ...)` to see more rows

下位の映画は:

> m_summary |>
     dplyr::filter(stringr::str_detect(variable, "^search_movie\\[")) |>
     dplyr::bind_cols(
         m = movie_master$movie
     ) |>
     dplyr::arrange(mean) |>
     dplyr::select(m, mean, q5, q95)
# A tibble: 1,664 × 4
   m                                             mean    q5   q95
   <chr>                                        <dbl> <dbl> <dbl>
 1 Mostro, Il (1994)                            -4.39 -6.06 -2.70
 2 Wend Kuuni (God's Gift) (1982)               -4.38 -6.13 -2.72
 3 Some Mother's Son (1996)                     -4.33 -5.48 -3.15
 4 Tainted (1998)                               -4.16 -5.64 -2.72
 5 Homage (1995)                                -4.11 -5.73 -2.58
 6 Man from Down Under, The (1943)              -4.03 -5.71 -2.34
 7 Terror in a Texas Town (1958)                -3.98 -5.86 -2.09
 8 Baton Rouge (1988)                           -3.94 -5.60 -2.15
 9 Eye of Vichy, The (Oeil de Vichy, L') (1993) -3.89 -5.77 -1.94
10 Someone Else's America (1995)                -3.88 -5.36 -2.28
# ℹ 1,654 more rows
# ℹ Use `print(n = ...)` to see more rows

値が高い映画のリストには、『スクリーム (1996)』Scream (1996))、『コンタクト (1997)』Contact (1997))といった、当時の流行りや話題作が並んでいます。これらは、多くの人が劇場に足を運んだり、話題にしたりした映画でしょう。

興行収入はウィキペディアによると、『スクリーム (1996)』1.7億ドルに達しています。

一方、値が低い映画には、『Il Mostro (1994)』『Wend Kuuni (God's Gift) (1982)』のような、よりニッチな作品が含まれています。どれほどニッチなのかというと、両方ともまず日本語のウィキペディアのページがございません

映画人気度

次に、映画の人気度($\phi$)を確認しましょう。同じく、上位の映画は:

> m_summary |>
     dplyr::filter(stringr::str_detect(variable, "^movie_propensity\\[")) |>
     dplyr::bind_cols(
         m = movie_master$movie
     ) |>
     dplyr::arrange(-mean) |>
     dplyr::select(m, mean, q5, q95)
# A tibble: 1,664 × 4
   m                                                       mean    q5   q95
   <chr>                                                  <dbl> <dbl> <dbl>
 1 Schindler's List (1993)                                 4.58  4.31  4.84
 2 Wrong Trousers, The (1993)                              4.46  4.09  4.84
 3 Close Shave, A (1995)                                   4.43  4.01  4.83
 4 Shawshank Redemption, The (1994)                        4.31  4.06  4.56
 5 Casablanca (1942)                                       4.28  4.01  4.56
 6 Wallace & Gromit: The Best of Aardman Animation (1996)  4.19  3.75  4.62
 7 Usual Suspects, The (1995)                              4.09  3.86  4.32
 8 Star Wars (1977)                                        4.06  3.89  4.22
 9 Rear Window (1954)                                      3.93  3.65  4.20
10 Godfather, The (1972)                                   3.88  3.68  4.09
# ℹ 1,654 more rows
# ℹ Use `print(n = ...)` to see more rows

下位の映画は:

> m_summary |>
     dplyr::filter(stringr::str_detect(variable, "^movie_propensity\\[")) |>
     dplyr::bind_cols(
         m = movie_master$movie
     ) |>
     dplyr::arrange(mean) |>
     dplyr::select(m, mean, q5, q95)
# A tibble: 1,664 × 4
   m                                                  mean    q5   q95
   <chr>                                             <dbl> <dbl> <dbl>
 1 Witness (1985)                                    -4.82 -6.28 -3.28
 2 Girl in the Cadillac (1995)                       -4.75 -6.24 -3.21
 3 They Made Me a Criminal (1939)                    -4.62 -6.24 -3.13
 4 Touki Bouki (Journey of the Hyena) (1973)         -4.47 -5.83 -3.19
 5 Police Story 4: Project S (Chao ji ji hua) (1993) -4.41 -5.97 -2.90
 6 Good Morning (1971)                               -4.41 -5.94 -2.90
 7 King of New York (1990)                           -4.30 -6.29 -2.48
 8 Spanish Prisoner, The (1997)                      -4.25 -5.83 -2.80
 9 Further Gesture, A (1996)                         -4.23 -5.77 -2.64
10 Visitors, The (Visiteurs, Les) (1993)             -4.21 -5.51 -2.94
# ℹ 1,654 more rows
# ℹ Use `print(n = ...)` to see more rows

まず、上位にランクインした映画は、純粋にやばいです。『シンドラーのリスト(1993)』(Schindler's List (1993))、『スター・ウォーズ(1977)』(Star Wars (1977))、『ゴッドファーザー(1972)』(Godfather, The (1972))など、映画にまったく興味がない私でも、歴史に名を残す名作として認識している作品が並んでいます。

ウィキペディアのリンクを貼るのが面倒なくらい、説得力のある推定結果だと言えるでしょう。

映画ベクトルの第一次元

最後に、モデルが推定した最も重要度の高い第一次元の正体を確認するため、第一次元の値が最も高い映画と最も低い映画を抽出してみましょう。

まず、上位の映画は:

> m_summary |>
     dplyr::filter(stringr::str_detect(variable, "^movie_latent\\[")) |>
     dplyr::mutate(id = variable |> purrr::map(\(x){as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])})) |> 
     tidyr::unnest_wider(id, names_sep = "_") |> 
     dplyr::filter(id_2 == 1) |> 
     dplyr::bind_cols(m = movie_master$movie) |> 
     dplyr::arrange(-mean) |> 
     dplyr::select(m, mean, q5, q95)
# A tibble: 1,664 × 4
   m                                                           mean    q5   q95
   <chr>                                                      <dbl> <dbl> <dbl>
 1 Before the Rain (Pred dozhdot) (1994)                       3.66  2.29  4.99
 2 Burnt By the Sun (1994)                                     3.49  2.58  4.41
 3 Pather Panchali (1955)                                      3.20  1.54  4.90
 4 Thirty-Two Short Films About Glenn Gould (1993)             3.20  2.30  4.10
 5 Aparajito (1956)                                            3.17  1.57  4.78
 6 Living in Oblivion (1995)                                   2.89  1.98  3.84
 7 Top Hat (1935)                                              2.89  1.50  4.35
 8 Horseman on the Roof, The (Hussard sur le toit, Le) (1995)  2.89  1.73  4.11
 9 Down by Law (1986)                                          2.78  1.95  3.60
10 Beat the Devil (1954)                                       2.78  1.15  4.41
# ℹ 1,654 more rows
# ℹ Use `print(n = ...)` to see more rows

下位の映画は:

> m_summary |>
     dplyr::filter(stringr::str_detect(variable, "^movie_latent\\[")) |>
     dplyr::mutate(id = variable |> purrr::map(\(x){as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])})) |> 
     tidyr::unnest_wider(id, names_sep = "_") |> 
     dplyr::filter(id_2 == 1) |> 
     dplyr::bind_cols(m = movie_master$movie) |> 
     dplyr::arrange(mean) |> 
     dplyr::select(m, mean, q5, q95)
# A tibble: 1,664 × 4
   m                         mean    q5    q95
   <chr>                    <dbl> <dbl>  <dbl>
 1 Star Kid (1997)          -2.81 -4.88 -0.897
 2 Radioland Murders (1994) -2.80 -4.15 -1.44 
 3 That Old Feeling (1997)  -2.32 -3.37 -1.31 
 4 Three Wishes (1995)      -2.26 -3.47 -1.09 
 5 Prefontaine (1997)       -2.22 -4.25 -0.185
 6 War, The (1994)          -2.09 -3.18 -0.917
 7 Prophecy, The (1995)     -2.07 -3.02 -1.17 
 8 Sixth Man, The (1997)    -2.05 -3.13 -1.03 
 9 Just Cause (1995)        -2.05 -2.78 -1.35 
10 Dangerous Minds (1995)   -2.05 -2.71 -1.41 
# ℹ 1,654 more rows
# ℹ Use `print(n = ...)` to see more rows

この次元で高い値を持つ映画には、『ビフォア・ザ・レイン (1994)』Before the Rain (Pred dozhdot) (1994))、『太陽に灼かれて
(1994)』
Burnt By the Sun (1994))、『大地のうた (1955)』Pather Panchali (1955))といった、ヨーロッパやアジアのアート映画、あるいは歴史的な名作が多く見られます。

逆に低い値を持つ映画には、『スター・キッズ (1997)』Star Kid (1997))や『笑撃生放送! ラジオ殺人事件 (1994)』Radioland Murders (1994))といった、より商業的で主流なアメリカのジャンル映画が並んでいます。

この結果から、モデルは「アート志向の映画」と「商業志向の映画」という、人々の映画の好みを分ける重要な軸をデータから自動で発見したと解釈できます。

では、最後に、映画認知度と映画人気度と映画ベクトルの第一次元の散布図を出して、相関があるかを確認しましょう:

g_search_propensity <- m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^search_movie\\[")) |>
  dplyr::select(認知度 = mean) |>
  dplyr::bind_cols(
    人気度 = m_summary |>
      dplyr::filter(stringr::str_detect(variable, "^movie_propensity\\[")) |>
      dplyr::pull(mean)
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_point(ggplot2::aes(x = 認知度, y = 人気度), color = ggplot2::alpha("blue", 0.5)) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

g_search_latent <- m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^search_movie\\[")) |>
  dplyr::select(認知度 = mean) |>
  dplyr::bind_cols(
    第一次元 = m_summary |>
      dplyr::filter(stringr::str_detect(variable, "^movie_latent\\[")) |>
      dplyr::mutate(id = variable |> purrr::map(\(x){as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])})) |> 
      tidyr::unnest_wider(id, names_sep = "_") |> 
      dplyr::filter(id_2 == 1) |>
      dplyr::pull(mean)
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_point(ggplot2::aes(x = 認知度, y = 第一次元), color = ggplot2::alpha("blue", 0.5)) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

g_propensity_latent <- m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^movie_propensity\\[")) |>
  dplyr::select(人気度 = mean) |>
  dplyr::bind_cols(
    第一次元 = m_summary |>
      dplyr::filter(stringr::str_detect(variable, "^movie_latent\\[")) |>
      dplyr::mutate(id = variable |> purrr::map(\(x){as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])})) |> 
      tidyr::unnest_wider(id, names_sep = "_") |> 
      dplyr::filter(id_2 == 1) |>
      dplyr::pull(mean)
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_point(ggplot2::aes(x = 人気度, y = 第一次元), color = ggplot2::alpha("blue", 0.5)) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")


library(patchwork)

g_search_propensity | g_search_latent | g_propensity_latent

parameters.png

あまり明確な関係があるとはいえなさそうですね。

このモデルはどんなドメインで利用できるか?

では、このモデルはどんなドメインで活用できるのでしょうか?まずは、モデルがどのようにパラメータの情報を取得しているかを理解しましょう。

本モデルは協調フィルタリングの変形であり、他の協調フィルタリングモデルと同様に、同じユーザーが複数のアイテムを評価し、かつ同じアイテムが複数のユーザーに評価されるという情報の相互作用を利用して、パラメータを推定しています。

具体的にデータを集計してみると以下のようになります:

> movie_df |>
     dplyr::mutate(
         point_flag = dplyr::case_when(
             point > 0 ~ 1,
             TRUE ~ 0
         )
     ) |>
     dplyr::summarise(
         movie_seen = sum(point_flag), .by = user_id
     ) |>
     dplyr::summarise(
         user_count = length(unique(user_id)), .by = movie_seen
     ) |>
     dplyr::arrange(dplyr::desc(user_count)) |>
     dplyr::mutate(
         user_ratio = user_count/sum(user_count),
         user_cumulative_ratio = cumsum(user_ratio)
     ) 
# A tibble: 279 × 4
   movie_seen user_count user_ratio user_cumulative_ratio
        <dbl>      <int>      <dbl>                 <dbl>
 1         20         31     0.0329                0.0329
 2         23         22     0.0233                0.0562
 3         22         22     0.0233                0.0795
 4         21         21     0.0223                0.102 
 5         27         18     0.0191                0.121 
 6         26         18     0.0191                0.140 
 7         25         17     0.0180                0.158 
 8         29         17     0.0180                0.176 
 9         24         17     0.0180                0.194 
10         33         16     0.0170                0.211 
# ℹ 269 more rows
# ℹ Use `print(n = ...)` to see more rows

一人当たりの評価数の分布を見ると、最も多いユーザーは 20 作品に評価をつけており、そうしたユーザーが 31 人存在します。さらに累積比率で見ると、上位 20% のユーザーは最低でも 20 本以上の映画を評価していることが分かります。

映画側に何人のユーザーが評価したかの分布も確認しましょう:

> movie_df |>
     dplyr::mutate(
         point_flag = dplyr::case_when(
             point > 0 ~ 1,
             TRUE ~ 0
         )
     ) |>
     dplyr::summarise(
         user_count = sum(point_flag), .by = movie
     ) |>
     dplyr::summarise(
         movie_count = length(unique(movie)), .by = user_count
     ) |>
     dplyr::arrange(dplyr::desc(movie_count)) |>
     dplyr::mutate(
         movie_ratio = movie_count/sum(movie_count),
         movie_cumulative_ratio = cumsum(movie_ratio)
     )
# A tibble: 272 × 4
   user_count movie_count movie_ratio movie_cumulative_ratio
        <dbl>       <int>       <dbl>                  <dbl>
 1          1         136      0.0817                 0.0817
 2          2          67      0.0403                 0.122 
 3          4          63      0.0379                 0.160 
 4          3          58      0.0349                 0.195 
 5          5          51      0.0306                 0.225 
 6          7          44      0.0264                 0.252 
 7          6          39      0.0234                 0.275 
 8         10          33      0.0198                 0.295 
 9          9          33      0.0198                 0.315 
10          8          29      0.0174                 0.332 
# ℹ 262 more rows
# ℹ Use `print(n = ...)` to see more rows

ユーザー一人あたりの評価数に比べるとやや少ないものの、一人のユーザーだけが評価している映画は全体の約8%に過ぎません。つまり、大半の映画には複数のユーザーから評価がついている状態です。

では、このように、よって、この MovieLens のデータは、同じユーザーが複数のアイテムを評価し、かつ同じアイテムが複数のユーザーに評価されるという条件を満たしていると言えよう。

ではこの条件を満たしていないドメインを見てみましょう。

  • 不動産
    単刀直入に聞きますが、あなたは生涯で何件の不動産を購入できますか?おそらく多くても一件程度でしょう。そのため、同じユーザーが複数のアイテムを評価するという条件は満たされません。

  • 求人市場
    こちらも率直に質問します。あなたはアルバイトや正社員など雇用形態を問わず、生涯で何回仕事を変えますか?おそらく1回以上はあるものの、MovieLensのユーザーのように20回も転職することは稀でしょう(ジョブホッパーすぎますね💦)。不動産よりは多いものの、やはり同じユーザーが複数のアイテムを評価する条件を満たしづらいです。

  • 営業の担当顧客リスト
    営業は非常に属人的な職種です。そのため、会社としては営業の顧客担当リストを頻繁に変えたくないのが本音です。つまり、同じユーザー(営業担当者)が複数のアイテム(顧客企業)を評価(商談結果)するものの、同じアイテムが複数のユーザーに評価される条件は満たされません。顧客は担当営業の変更を滅多に経験しないためです。

以上の説明からわかるように、本モデルの前提となる

  • 同じユーザーが複数のアイテムを評価する
  • 同じアイテムが複数のユーザーに評価される

という条件は、MovieLens の映画評価データのようなドメインでは十分に満たされています。この相互の評価情報がモデル推定の鍵となっており、協調フィルタリングの効果的な適用を可能にしています。

一方で、不動産や求人市場、営業の担当顧客リストのように、ユーザーやアイテムの評価の偏りが大きく、どちらか一方の条件を満たしにくいドメインでは、このモデルの適用は難しいと言えます。

したがって、本モデルを活用するには、ユーザーとアイテムの間で十分な評価の重複が存在し、双方の情報が豊富に得られるデータ環境が不可欠です。 本モデルを実際のビジネスに応用する際は、まずこの基本的なデータ構造を確認することを強く推奨します。

結論

いかがでしたか?

本記事では、レコメンドエンジンの推定結果を活用し、事業における資源の最適配分を考えるための基盤として、MovieLens データを用いた解釈可能なベイズ機械学習モデルを構築しました。このモデルは、ユーザーの未評価の理由を「そもそも見ていない」か「見たけど評価しなかった」かに切り分け、データから自動的に次元数を決定するノンパラメトリックなアプローチを採用しています。

推定結果は、私たちのモデルが以下の重要な示唆を導き出したことを示しています。

  • 単一の支配的な潜在次元の発見
    モデルは無限の次元を仮定しましたが、結果として評価構造の大部分がたった一つの次元で説明できると判明しました。この次元は、「商業志向の映画」と「アート志向の映画」という、人々の映画の好みを分ける最も根本的な軸を捉えていると解釈できます。この知見は、モデルの軽量化や推論の高速化といった技術的なメリットに加え、コンテンツ戦略をシンプルかつ効果的に策定するための重要なビジネス的示唆をもたらします。

  • ユーザー嗜好の一様性
    ユーザーの潜在嗜好は、政治分野でよく見られるような二極化ではなく、単一の釣鐘型の分布を示しました。このことは、極端な嗜好を持つユーザーが少なく、多くのユーザーに受け入れられる普遍的なコンテンツが重要であることを示唆しています。リソースを特定のニッチに偏らせるのではなく、幅広い層にアピールする「中央寄せ」のコンテンツ戦略が有効である可能性が高いでしょう。

  • 「人気」と「質」の分離
    モデルは、映画の「認知度」と「人気度」を明確に区別しました。例えば、高評価のリストには『シンドラーのリスト』や『スター・ウォーズ』といった時代を超えた名作が並び、認知度の高いリストには『スクリーム』や『コンタクト』といった当時のヒット作がランクインしました。この分離により、単なる人気に惑わされず、本当にユーザーの心を掴むコンテンツを特定する洞察が得られます。

本モデルは、ドメイン知識が乏しい私でも、データから本質的な構造を自動で抽出し、ビジネスに直結する示唆を提供できることを示しました。次回の記事では、今回得られた知見と推定結果を土台に、実際に限られたリソースをユーザーにどう配分するか、具体的な手法と戦略についてさらに深く掘り下げていきます。ご期待ください。

最後に、私たちと一緒に働きたい方はぜひ下記のリンクもご確認ください:

参考文献

Wu, Tung-Wen. (2025). Estimating Dynamic Dimensionality in Roll Call Data: An Application to UNGA Voting. Poster presented at the Japanese Society for Quantitative Political Science (JSQPS) Summer Meeting, Kyoto, July 2025.

2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?