前言

您好，我是一名在事業公司工作的資料科學家。

本文為以下文章的續篇，重點將放在「如何實際運用」已推估的貝葉斯機器學習模型，尤其是所謂的因子分析系模型。

首先大前提是，本文假設不需要報告信用區間的分析任務。因此，如果法律或內部/組織規定要求嚴格報告信用區間，請勿使用本文所介紹的手法。

接下來，作為一名在民間企業工作近5年的資料科學家，並透過學會發表等形式維持與學術界的關係，強烈感受到政治學方法論、計量經濟學等統計學領域與商業之間存在著龐大的鴻溝。其中關鍵在於，統計學基本上並不假設「運用模型」。

此處所謂「運用」，指的是持續利用一次推估的模型並對學習時不存在的新數據做出判斷。

當然，政治學領域中也有如 DW-NOMINATE 或 V-Dem 等持續更新數據與推估結果的優秀專案。然而，從整體來看，許多情況下論文發表時分析已告一段落，而隨後的運用則不被重視。

在學術界，常常是「論文發表了，好，結束」；而在商業現場，則需要在每週的經營會議上更新指標並報告，或為新用戶推薦影片，持續使用模型。因此，這種差距可說是那些從學術界轉向商業界的人最初面臨的重大挑戰之一。

本文將以因子分析模型為軸心，整理「運用到底是什麼」的議題，並解釋在考慮運用時容易陷入的陷阱，最後介紹可在實務中使用的具體運用方法。

正如標題所示，本文中介紹的運用手法，在嚴格遵循貝葉斯原則的立場下，可能會被視為對統計學的背叛，因為它重視信用區間等參數不確定性的價值。

儘管如此，這依然是我經過試行錯誤後所達成的，在現實業務中穩定運用貝葉斯因子分析模型的方法。這是優先考慮持續運行和不崩壞使用的結果，而非理論之美。

如果您非常希望將貝葉斯原則視為最優先的考量，請務必發展出能正面解決本文所指的「運用上問題」的更精緻新手法。

現在，就讓我們直接進入主題吧！

運用是什麼

無論是在學術界還是企業中，基本上在資料科學專案中「推估某個模型，並報告該模型所得到的啟發或評估其性能」的情況都是相當普遍的。使用的方法論之間也有許多共同點，在這一點上學術界與商業界的距離並不像想像中那麼遙遠。

然而，後續發展的結果卻截然不同。在學術界，獲得的啟發或性能會以論文的形式整理，一旦發表並被採納便宣告結束；而在商業界，特別是在以平均處置效果推估為目的的因果推論與效果檢驗以外的專案中，幾乎可以確定會收到以下的回應。

• 「哇，商談成約率的預測模型準確度真高。那麼，讓我們建立一個系統，以便每日更新給銷售部門」
• 「原來用戶的喜好被提取出來了。那麼我們要儘快與後端團隊討論，將其納入商品顯示邏輯中」

這就是這裡所謂的運用。能事先想像到這種後續行動並不感到驚訝，是我認為自己已經開始適應商業環境的第一個定性指標。

進入運用階段後，與工程師的溝通會大幅增加。之前有篇文章講述過如何與工程師有效溝通，如果您感興趣，請務必參考。

不過，實際上問題不僅僅是溝通。假設我們以商談成約率的預測模型為例。為了簡化問題，我們假設推估出下列模型：

$$
\hat{成約率{i}} = \hat{f}(商談回數{i}, 業種_{i})
$$

在這種情況下，當未曾接觸過的業種客戶來時，怎樣預測的問題存在，但對於這點有眾多應對方法，因此不贅述。簡而言之，從事前定義的特徵量進行預測的任務，在運用時通常不會面臨本質上的困難。

那麼，對於因子分析及其它模型來說又如何呢？每當新用戶進來時，需要在同樣的模型框架中估計該用戶的偏好所表徵的向量。若您來自學術界，您會首先想到使用包含新用戶的所有用戶數據重新推估模型。這在統計學上確實是一個非常正確的判斷。

然而，這種做法存在兩個重大問題。

第一，計算時間問題。若從零開始重新推估，將耗費大量時間。稍後將提供具體數字，但單單為了一位新用戶而重新推估模型的方法，比本文所提議的方法約需多花57倍的時間。而且，即使使用本文所提及的方法，與整體模型的整合性依然能夠保持良好。那麼，為了執行用戶推薦等可能不使用信用區間的實務任務，真的要將57倍的時間花費在僅僅是統計學的嚴謹性上嗎？你有勇氣在經營會議上向董事或高層解釋這一判斷標準嗎？請從經營的角度思考這一問題。

第二個問題則是因子的穩定性。像主題模型、因子模型、協同過濾等因子分析類模型，當學習數據變更時，推估出的因子也可能會略有變化。因此，在原有推估結果的第一次元與包含新數據後重新推估的第一次元，並不能保證它們表現出相同的趨勢。這可能是由於學習算法的隨機性質而產生，也可能是因為新數據向模型提供新信息，使得結構本身受到了重新評估。

後者實際上是因子模型的優勢，但例如「第一次元與我司的廣告收入強相關，因此我要在每週的經營會議中分享第一次元為負的用戶數量」等運用一旦開始，模型自動進行的結構重新評估，對現場來說可能成為多餘的麻煩。

更棘手的是，因子模型具有符號的不確定性，每次估計可能導致因子的正負反轉。

當然，剛從學術界來的人可能會直接拒絕這種監控因子的請求。就我自己而言，起初也曾考慮拒絕。然而，為了正面應對商業挑戰，我決定思考解決方案。不限制方法論能解決的課題，而是當必須解決的課題存在時，考慮相應的方法。本文即是這一嘗試所形成的結果。

運用的智慧：將參數當作數據處理

以下內容假設您已閱讀前提的文章，請您先行參閱！

這裡所介紹的運用的智慧是，當新用戶流入時，已推估的「維度相關」「電影相關」「用戶群集相關」的各參數將事後平均值直接固定為常數，並排除在學習對象外。基於此，僅針對新用戶的參數進行推估。

這種方式首先大幅減少了新增用戶時所需學習的參數數量。而且，維度、電影、用戶群集相關的參數都來自現有模型的常數，因此，模型的表現能力將被主動限制，這樣一來，新用戶的參數必然會在原始模型的整合性位置內進行推估。

另一方面，需要注意的是，這種方法完全忽略了原模型中的參數不確定性。若從嚴格遵循貝葉斯原則的立場來看，重視信用區間等參數不確定性的價值，這可能會被視為背叛統計學。實際上，與新用戶相關的參數的變異量過小的問題無法避免。

然而，正如稍後所示的性能驗證所示，只要著眼於參數的平均值，這種運用方式能以相當高的精度重現原始模型。因此，在基於平均值做出決策的任務中，這種手法顯得非常實用且有效。

重現實驗

首先是本文中使用的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, matrix 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;

  vector[dimension_type] dimension;

  vector[group_type] group_user;
  vector[group_type] group_user_sigma;
  matrix[group_type, dimension_type] group_user_latent;

  real search_intercept;
  vector[movie_type] search_movie;

  vector[movie_type] movie_propensity;
  matrix[movie_type, dimension_type] movie_latent;

  int N;
  array[N] int user;
  array[N] int movie;
  array[N] int result;
}
parameters {
  array[user_type] vector[dimension_type] user_latent;

  array[user_type] ordered[result_type - 1] cutpoints;

  vector[user_type] search_user;
}
model {
  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);
  }

  search_user ~ 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
  );
}

與前提文章中的Stan程式碼相比，所有「維度相關」「電影相關」「用戶群集相關」的參數均以 data 的方式傳入模型。因此，推估新用戶參數時，這些值的更新絕對不會發生。這在嚴格上違反了貝葉斯的原則，但其代價是避免新維度突現或現有維度解釋變更，在運用的角度上優雅地提升了易用性。

不過，這種看似奇特的Stan寫法，是否真的能以原模型的整合性位置推估新用戶，當然會引起疑問。為此，本文將進行隨機單元測試來驗證此點。

具體來說，我們假設已推估的「維度相關」「用戶相關」「電影相關」「用戶群集相關」的參數已經掌握，而故意隱藏用戶相關的參數，以確認本文所提方法能以多大精度重現這些參數。

首先編譯模型：

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

執行推估：

m_nu_estimate <- m_nu_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,

         dimension = m_summary |>
             dplyr::filter(stringr::str_detect(variable, "^dimension\\[")) |>
             dplyr::pull(mean),
         group_user = m_summary |>
             dplyr::filter(stringr::str_detect(variable, "^group_user\\[")) |>
             dplyr::pull(mean),
         group_user_sigma = m_summary |>
             dplyr::filter(stringr::str_detect(variable, "^group_user_sigma\\[")) |>
             dplyr::pull(mean),
         group_user_latent = m_summary |>
             dplyr::filter(stringr::str_detect(variable, "^group_user_latent\\[")) |>
             dplyr::pull(mean) |>
             matrix(ncol = 20),

         search_intercept = m_summary |>
             dplyr::filter(stringr::str_detect(variable, "^search_intercept")) |>
             dplyr::pull(mean),
         search_movie = m_summary |>
             dplyr::filter(stringr::str_detect(variable, "^search_movie\\[")) |>
             dplyr::pull(mean),

         movie_propensity = m_summary |>
             dplyr::filter(stringr::str_detect(variable, "^movie_propensity\\[")) |>
             dplyr::pull(mean),
         movie_latent = m_summary |>
             dplyr::filter(stringr::str_detect(variable, "^movie_latent\\[")) |>
             dplyr::pull(mean) |>
             matrix(ncol = 20),

         N = nrow(movie_df_with_id_downsampled[-test_id,]),
         user = movie_df_with_id_downsampled$user_id[-test_id],
         movie = movie_df_with_id_downsampled$movie_id[-test_id],
         result = movie_df_with_id_downsampled$point[-test_id] + 1
     )
)

所需時間分別僅為 1.7 秒和 5.3 秒。相比於在前篇文章中從零開始推估所花的 285.8 秒，兩者之間的差距顯而易見。

即使在本方法上採取較為苛刻的條件（較慢的 5.3 秒）來作比較，將現有的各參數（維度、電影、集群）固定常數，僅有限地推估新用戶的參數，這種做法實現了對全數據重新推估方法的53.92倍的壓倒性計算時間縮減。

那麼，推估的參數與原模型是否相符呢？我們首先檢查藝術電影的用戶：

m_aml_summary

用戶的喜好可被彙整為「藝術取向的電影」（正向）與「商業取向的電影」（負向），原模型的結論是，僅對藝術取向的電影前十作品滿分的用戶，其第一次元得分如預期偏向正向，與原模型結果相符。

那麼商業電影的用戶又是如何呢？

m_uml_summary

商業取向的電影得分偏向負向，但對商業取向的電影前十作品滿分的用戶的第一次元得分則如預期偏向負向，同樣與原模型結果相符。

由此可見，至少在本模型中，此方法的有效性得以證明。然而，這不能保證在所有因子分析模型中都能直接應用。

在其他模型中採取本方法時，務必事先驗證推估參數與原模型估計值的充分一致性。此外，也強烈建議事先使用類似本次實施的「商業電影派 vs 藝術電影派」的推估行為可預測的「極端虛擬數據」進行模擬。

總結

本文介紹了一種為了「作為系統運用」而部分推估新用戶的手法，固定現有參數為常數。

這種手法在某種意義上背叛了貝葉斯統計學本應重視的「不確定性的傳播」和「參數的同時推估」等大原則。從重視統計學嚴謹性的人員的角度來看，這可能是一種不完善的手法。

然而，實務所需的並不是「花費50倍以上計算資源換來的100分的精準性」，而是「在數秒內回應結果，與昨天今天保持一致的80分的穩定性」也是事實。

此外，我們常聽到「貝葉斯的真正價值在於能夠量化不確定性」的說法，但換句話說，這種認知已經過時。當代貝葉斯統計學的本質在於對眼前的問題靈活建構模型，鮮明地可視化數據中的結構，這種建模的自由度才是真正的關鍵。為了理解其一端，請務必觀看這段影片的48分鐘後的內容：

當然，這種「運用的智慧」並不是萬能的。若數據的分布發生戲劇性變化，則可能需重新檢討（再學習）固定參數本身。

重要的是，不應限制方法論能解決的問題，應考慮到若有必須解決的問題，則需針對該問題思考相应的方式。

希望本文能成為那些勇敢從學術界踏入商業浪潮，並在理想與現實之間徘徊的資料科學家的武器。

讓我們在堅守理論的同時，繼續在現場戰鬥。

最後，對於有意與我們共事的人，請務必查閱以下連結：

原文出處：https://qiita.com/Gotoubun_taiwan/items/8010f0d62025db8f3d81