はじめに
こんにちは、事業会社で働いているデータサイエンティストです。
本記事は、下記の記事の続編です：
前回の記事では、大量の時系列データを分類するベイズモデルを紹介しました。ただし、分類モデルでは時系列同士の類似度を抽出する用途には適していないという課題が残ります。
そこで本記事では、時系列を低次元空間に埋め込み、さらにその次元数自体も推定可能な手法を提案します。
このように時系列を低次元に埋め込むことで、コサイン類似度など一般的な指標を用いて類似度を測ることが可能となり、加えて、時系列の背後にある潜在的なトレンドも同時に抽出できるようになります。
モデル説明
今回のモデルも、前回と同様にディリクレ過程の構造を用いています。 ただし、ディリクレ過程を直接分類に使うのではなく、その中間産物である棒折り過程の重みを活用し、これをソフトな次元削減に応用します。
具体的には、「A次元は使う」「B次元は使わない」といった二値的な判断ではなく、「A次元の重みは0.82」「B次元の重みは0.11」のように、各次元に連続的な重みを割り当てる形を取ります。このようにして、モデルは必要と判断された次元に高い重みを与え、不要な次元には低い重みを与えることで、次元削減に近い構造を実現しています。
では、モデルの説明に入ります。
まず、ハイパーパラメータ$\alpha$をガンマ分布からサンプリングします：
$$
\alpha \sim Gamma(0.001, 0.001)
$$
次に、棒折り過程を用いて次元pの重み$p_{p}$をサンプリングします。次元は無限にある設定なので、理論上この操作は無限回繰り返します：
$$
\pi_{p} \sim Beta(1, \alpha)
$$
$$
p_{p} = \pi_{p} \prod_{l=1}^{p - 1} (1 - \pi_{l})
$$
次に、次元pのトレンド$\mu_{d}$は一次の自己回帰モデルの形で表現します：
$$
\phi_{begin, p} \sim Gamma(0.001, 0.001)
$$
$$
\phi_{shift, p} \sim Gamma(0.001, 0.001)
$$
$$
\mu_{p, 1} \sim Normal(0, \phi_{begin, p})
$$
$$
\mu_{p, t} \sim Normal(\mu_{p, t-1}, \phi_{shift, p}), \space \text{for t > 1}
$$
続いて、時系列iの無限次元の埋め込み表現$\beta_{i}$は、正規分布からサンプリングされます：
$$
\beta_{i} \sim Normal(0,10)
$$
最後に、時系列iの時間tの値$series_{i, t}$は、このようにサンプリングされます：
$$
series_{i, t} \sim Normal(intercept_{i} + \sum_{d=1}^{\infty}\beta_{i,d}\mu_{d,t}p_{d}, \sigma)
$$
最後の$\sum_{d=1}^{\infty}\beta_{i,d}\mu_{d,t}p_{d}$を丁寧に展開し：
$$
\beta_{i,1}\mu_{1,t}p_{1} + \beta_{i,2}\mu_{2,t}p_{2} + \beta_{i,3}\mu_{3,t}p_{3} + ...+ \beta_{i,\infty}\mu_{\infty,t}p_{\infty}
$$
さらに、$\mu_{d,t}p_{d}$を$\hat{\mu_{d,t}}$と定義すると、
$$
\beta_{i,1}\hat{\mu_{1,t}} + \beta_{i,2}\hat{\mu_{2,t}} + \beta_{i,3}\hat{\mu_{3,t}} + ...+ \beta_{i,\infty}\hat{\mu_{\infty,t}}
$$
この結果、パラメータが$\beta$、独立変数が$\mu$である無限次元の線形回帰モデルとして解釈できることがわかります。通常の線形回帰モデルでは、独立変数を明示的に指定する必要がありますが、本モデルでは、データの背後に潜む傾向やトレンドを自動的に独立変数として抽出する構造になっています。
また、前回の記事と同じように、理論上は次元数は無限に存在し得ますが、現実のコンピュータは無限長の配列を扱うことができません。そのため、実装上は「無限大」をあらかじめ十分に大きな有限の数値で近似し、計算可能な範囲に収める必要があります。
モデルの実装コードはこちらです：
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 weighted_inner_product(
vector vec_a, vector vec_b, vector weight
){
vector[size(vec_a)] result;
for (i in 1:size(vec_a)){
result[i] = vec_a[i] * vec_b[i] * weight[i];
}
return sum(result);
}
real partial_sum_lpdf(
array[] real value,
int start, int end,
array[] int date, array[] int variable,
vector dimension,
matrix mu, matrix beta, vector intercept, real sigma
){
vector[end - start + 1] log_likelihood;
int count = 1;
for (i in start:end){
log_likelihood[count] = normal_lpdf(value[count] | intercept[variable[i]] + weighted_inner_product(mu[:,date[i]], beta[variable[i],:]', dimension), sigma);
count += 1;
}
return sum(log_likelihood);
}
}
data {
int date_type;
int variable_type;
int dimension_type;
int N;
array[N] int date;
array[N] int variable;
array[N] real value;
}
parameters {
real<lower=0> dimension_alpha; // ディリクレ過程の全体のパラメータ
vector<lower=0, upper=1>[dimension_type - 1] dimension_breaks; // ディリクレ過程のstick-breaking representationのためのパラメータ
vector<lower=0>[dimension_type] phi_begin;
vector<lower=0>[dimension_type] phi_shift;
real<lower=0> sigma;
matrix[dimension_type, date_type] mu_unnormalized;
vector[variable_type] intercept;
matrix[variable_type, dimension_type] beta;
}
transformed parameters {
simplex[dimension_type] dimension;
matrix[dimension_type, date_type] mu;
dimension = stick_breaking(dimension_breaks);
for (i in 1:dimension_type){
mu[i,:] = mu_unnormalized[i,:] - mean(mu_unnormalized[i,:]);
}
}
model {
dimension_alpha ~ gamma(0.001, 0.001);
dimension_breaks ~ beta(1, dimension_alpha);
phi_begin ~ gamma(0.001, 0.001);
phi_shift ~ gamma(0.001, 0.001);
sigma ~ gamma(0.001, 0.001);
for (i in 1:dimension_type){
mu[i, 1] ~ normal(0, phi_begin[i]);
for (j in 2:date_type){
mu[i, j] ~ normal(mu[i, j - 1], phi_shift[i]);
}
}
intercept ~ normal(0, 10);
to_vector(beta) ~ normal(0, 10);
int grainsize = 1;
target += reduce_sum(
partial_sum_lupdf, value, grainsize,
date, variable,
dimension,
mu, beta, intercept, sigma
);
}
generated quantities {
matrix[variable_type, variable_type] cos_sim;
{
matrix[variable_type, dimension_type] weighted_beta;
for (i in 1:variable_type){
for (j in 1:dimension_type){
weighted_beta[i, j] = beta[i, j] * dimension[j];
}
}
for (i in 1:variable_type){
for (j in 1:variable_type){
cos_sim[i, j] = (weighted_beta[i,:] * weighted_beta[j,:]')/sqrt(sum(weighted_beta[i,:].^2) * sum(weighted_beta[j,:].^2));
}
}
}
}
モデル推定
ここからは、実際のデータ分析に進みます。使用するのは、アメリカのマクロ経済に関連する時系列データです。このデータには一部欠損値が含まれていますが、本モデルの設計上、問題なく取り扱うことが可能です。
前処理は前の記事とほぼ一緒です。
data("fred_md", package = "BVAR")
ts_df <- fred_md |>
tibble::rownames_to_column(var = "date") |>
tibble::tibble() |>
tidyr::pivot_longer(!date, names_to = "variable", values_to = "value") |>
tidyr::drop_na() |>
dplyr::group_by(variable) |>
dplyr::mutate(
value = (value - mean(value))/sd(value),
date = as.Date(date)
) |>
dplyr::ungroup() |>
dplyr::arrange(variable, date)
date_master <- ts_df |>
dplyr::select(date) |>
dplyr::distinct() |>
dplyr::arrange(date) |>
dplyr::mutate(
date_id = dplyr::row_number()
)
variable_master <- ts_df |>
dplyr::select(variable) |>
dplyr::distinct() |>
dplyr::arrange(variable) |>
dplyr::mutate(
variable_id = dplyr::row_number()
)
ts_df_with_id <- ts_df |>
dplyr::left_join(
date_master, by = "date"
) |>
dplyr::left_join(
variable_master, by = "variable"
)
data_list <- list(
date_type = nrow(date_master),
variable_type = nrow(variable_master),
# 無限大を10で近似
dimension_type = 10,
N = nrow(ts_df_with_id),
date = ts_df_with_id$date_id,
variable = ts_df_with_id$variable_id,
value = ts_df_with_id$value
)
続いて、モデルのコンパイルを実施し：
m_ts_emb_init <- cmdstanr::cmdstan_model("ts_emb.stan",
cpp_options = list(
stan_threads = TRUE
)
)
変分推論を用いて、モデルを推定します：
> m_ts_emb_estimate <- m_ts_emb_init$variational(
seed = 12345,
threads = 16,
data = data_list
)
------------------------------------------------------------
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.01188 seconds
1000 transitions using 10 leapfrog steps per transition would take 118.8 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)
Success! Found best value [eta = 1] earlier than expected.
Begin stochastic gradient ascent.
iter ELBO delta_ELBO_mean delta_ELBO_med notes
100 -391383.427 1.000 1.000
200 -130630.821 1.498 1.996
300 -76558.960 1.234 1.000
400 -75998.109 0.927 1.000
500 -75540.845 0.743 0.706
600 -75662.680 0.620 0.706
700 -75530.389 0.531 0.007 MEDIAN ELBO CONVERGED
Drawing a sample of size 1000 from the approximate posterior...
COMPLETED.
Finished in 22.2 seconds.
最後に、結果を整形して保存します：
m_ts_emb_summary <- m_ts_emb_estimate$summary()
推定結果可視化
まず、推定された次元の重みを確認します：
> m_ts_emb_summary |>
+ dplyr::filter(stringr::str_detect(variable, "^dimension\\["))
# A tibble: 10 × 7
variable mean median sd mad q5 q95
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 dimension[1] 0.472 4.72e- 1 0.00220 2.28e- 3 4.69e- 1 0.476
2 dimension[2] 0.528 5.28e- 1 0.00220 2.28e- 3 5.24e- 1 0.531
3 dimension[3] 0.00000655 2.82e- 6 0.0000116 3.43e- 6 1.36e- 7 0.0000241
4 dimension[4] 0.00000382 1.33e- 6 0.00000815 1.77e- 6 3.81e- 8 0.0000145
5 dimension[5] 0.00000163 3.19e- 7 0.00000455 4.56e- 7 3.10e- 9 0.00000657
6 dimension[6] 0.000000477 4.30e- 8 0.00000147 6.28e- 8 2.67e-10 0.00000228
7 dimension[7] 0.000000197 7.31e- 9 0.00000103 1.08e- 8 1.73e-11 0.000000610
8 dimension[8] 0.0000000508 1.19e- 9 0.000000270 1.77e- 9 1.35e-12 0.000000220
9 dimension[9] 0.0000000173 2.03e-10 0.0000000896 3.01e-10 8.87e-14 0.0000000655
10 dimension[10] 0.00000000946 3.58e-11 0.000000119 5.30e-11 1.39e-14 0.0000000222
前回の記事と同じ結論で、二つの大きなトレンドがある結果になっています。また、最初の二つの次元以外の次元の重みはほぼゼロなので、しっかりソフトな次元削減がなされたことがわかります。
続いて、最初の二つの次元のトレンド$\mu$を可視化します：
m_ts_emb_summary |>
dplyr::filter(stringr::str_detect(variable, "^mu\\[")) |>
dplyr::mutate(
id = variable |>
purrr::map(
\(x){
as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])
}
)
) |>
tidyr::unnest_wider(id, names_sep = "_") |>
dplyr::left_join(date_master, by = c("id_2" = "date_id")) |>
dplyr::filter(id_1 %in% c(1,2)) |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = date, y = mean, color = as.factor(id_1))) +
ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95, fill = as.factor(id_1)), alpha = 0.3) +
ggplot2::scale_x_date(
breaks = seq(min(date_master$date), max(date_master$date), by = "5 year"),
date_labels = "%Y"
) +
ggplot2::labs(
color = "dimension",
fill = "dimension"
)
ほぼ前回の記事のモデルが抽出したトレンドですねw
最初の次元（dimension 1）は、前回の記事におけるグループ4とほぼ一致しています。また、二番目の次元（dimension 2）は、前回のグループ1に-1を掛けたものとほぼ同じ結果となっています。
「統計モデルは設定次第でどんな結果も出せる」という批判は、文系寄りの方や、因果推論・部分識別などモデルに依存しない手法を好む方々からよく聞かれます。ですが、ぜひ冷静に見てください。今回は、モデルの設定を「分類」から「埋め込み」に変更しました。それでも、抽出された隠れトレンドはほぼ同じでした。
残念ながら、「設定次第で何でも出せる」というのは、単にパッケージにデータを流し込むだけ、あるいは dplyr::filterを延々と使ってデータをいじるだけの、表層的なモデルに限った話であると言わざるを得ません。
一方で、データ生成過程に明示的な仮定を置いた緻密なモデルであれば、多少構造を調整しても、本質的に抽出される傾向が大きく変わることはないのです。
次に、時系列の埋め込み表現$\beta$を可視化します。t-SNE を使いたくなる気持ちはわかりますが、次元数がわずか2次元なので、ここではシンプルにそのままプロットすることにします。
m_ts_emb_summary |>
dplyr::filter(stringr::str_detect(variable, "^beta\\[")) |>
dplyr::mutate(
id = variable |>
purrr::map(
\(x){
as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])
}
)
) |>
tidyr::unnest_wider(id, names_sep = "_") |>
dplyr::left_join(variable_master |>
dplyr::rename(var = variable),
by = c("id_1" = "variable_id")) |>
dplyr::filter(id_2 %in% c(1, 2)) |>
dplyr::mutate(
id_2 = stringr::str_c("dimension_", id_2)
) |>
dplyr::select(var, id_2, mean) |>
tidyr::pivot_wider(names_from = id_2, values_from = mean) |>
ggplot2::ggplot() +
ggplot2::geom_text(ggplot2::aes(x = dimension_1, y = dimension_2, label = var), color = ggplot2::alpha("blue", 0.5)) +
ggplot2::xlim(c(-2, 2)) +
ggplot2::ylim(c(-2, 2))
プロットを見ると、下部に密集している時系列群と、右側にやや固まっている時系列群が確認できます。それ以外の時系列は、比較的まばらに分布していることがわかります。
最後に、よりシャープに時系列の類似度を確認するために、次元の重みでスケーリングされた埋め込み表現$\hat{\beta_{i}}$：
$$
\hat{\beta_{i, d}} = \beta_{i, d}p_{d}
$$
を用います。これにより、実際にモデルで重視された次元における、時系列同士のコサイン類似度を評価できます。
以下では、原油価格の時系列に注目して可視化します。
m_ts_emb_summary |>
dplyr::filter(stringr::str_detect(variable, "^cos_sim\\[")) |>
dplyr::mutate(
id = variable |>
purrr::map(
\(x){
as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])
}
)
) |>
tidyr::unnest_wider(id, names_sep = "_") |>
dplyr::left_join(variable_master |>
dplyr::rename(var_1 = variable),
by = c("id_1" = "variable_id")) |>
dplyr::left_join(variable_master |>
dplyr::rename(var_2 = variable),
by = c("id_2" = "variable_id")) |>
dplyr::filter(var_1 == "OILPRICEx") |>
ggplot2::ggplot() +
ggplot2::geom_point(ggplot2::aes(x = mean, y = reorder(var_2, mean)), color = "blue") +
ggplot2::geom_errorbarh(ggplot2::aes(xmin = q5, xmax = q95, y = var_2), position = "dodge", height = 0.2, color = ggplot2::alpha("blue", 0.5)) +
ggplot2::theme(axis.text.y = ggplot2::element_text(size = 4)) +
ggplot2::labs(
x = "cosine similarity",
y = ""
)
具体的に説明すると、原油価格（OILPRICEx）に非常に近い方向を持つ指標としては以下のようなものがあります：
- DTCOLNVHFNM（金融会社が保有する消費者向け自動車ローン）
- M2REAL（実質マネーサプライ）
- INVEST（投資）
- UEMP15T26（失業期間が15〜26週間の失業者数）
これらは、原油価格の上昇と同時に動く傾向がある、つまりマクロ経済における共通の需要・景気循環的な要因で駆動されていると解釈できます。
一方、原油価格と明確に反対方向に動く指標も存在します：
- ISRATIOx（企業の在庫売上比率）
- HOUST（住宅着工件数）
- FEDFUNDS（政策金利）
- GS10（10年物国債利回り）
これらは、原油価格が上がると逆に下がる傾向があるもので、例えば原油高によるインフレ圧力が在庫調整や住宅市場の冷え込み、または政策金利の上昇を招いているという構造的関係を示唆します。
「相関係数でいいじゃん」と結論
「相関係数でいいじゃん？」という意見は、一部のビジネスサイドや、リーダブルコードを重視するエンジニアから出てくるかもしれません。そこで、実際に比較してみる前に、まず一つ、残酷な事実をお伝えしなければなりません。このアメリカのマクロ経済データには、欠損値が存在します。つまり、欠損がある時系列と原油価格との相関係数は計算できません。
さらに現実のビジネス環境では、欠損値の存在に加えて、集計期間の不一致といった問題も頻繁に発生します。このような背景を考えると、相関係数ベースのアプローチは、非常に初歩的な段階で大きな壁に直面するのです。
「モデルはビジネスの現場では使えない」と主張する方へ申し上げます。残念ながら、そのような主張をされる方こそ、ビジネスの現場で何が起きているかを理解していないと言わざるを得ません。
実際の現場では、欠損値や集計期間のずれは頻発します。それらを考慮せずに「相関係数で十分」「モデルは複雑すぎる」と言うのは、現実を見ていない証拠です。もう一度、階級や年齢を忘れて、一度現場に足を運んでみてはいかがでしょうか？ 本当に必要とされているのは、そうした「現場感覚」に根ざしたモデルなのです。
本記事で提案されたモデルで推定した埋め込み表現を使えば、欠損があっても、集計期間の不一致があっても、コサイン類似度は理論的に正しく計算できます。欠損などの影響は、コサイン類似度の事後分布の信用区間の幅として現れます。
さて、話が長くなりましたが、相関係数とコサイン類似度の関係を可視化します。相関係数が定義できないデータは除外されます：
fred_md |>
ncol() |>
seq_len() |>
purrr::map_dbl(
\(i){
cor(
fred_md[, "OILPRICEx"] - mean(fred_md[, "OILPRICEx"])/sd(fred_md[, "OILPRICEx"]),
fred_md[, i] - mean(fred_md[, i])/sd(fred_md[, i]),
)
}
) |>
tibble::tibble(
correlation = _
) |>
dplyr::bind_cols(
variable_name = colnames(fred_md)
) |>
dplyr::left_join(
m_ts_emb_summary |>
dplyr::filter(stringr::str_detect(variable, "^cos_sim\\[")) |>
dplyr::mutate(
id = variable |>
purrr::map(
\(x){
as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])
}
)
) |>
tidyr::unnest_wider(id, names_sep = "_") |>
dplyr::left_join(variable_master |>
dplyr::rename(var_1 = variable),
by = c("id_1" = "variable_id")) |>
dplyr::left_join(variable_master |>
dplyr::rename(var_2 = variable),
by = c("id_2" = "variable_id")) |>
dplyr::filter(var_1 == "OILPRICEx") |>
dplyr::arrange(dplyr::desc(mean)) |>
dplyr::select(var_2, cos_sim = mean),
by = c("variable_name" = "var_2")
) |>
ggplot2::ggplot() +
ggplot2::geom_text(ggplot2::aes(x = cos_sim, y = correlation, label = variable_name), color = ggplot2::alpha("blue", 0.3))
このように、結論大きく変わりませんが、モデルの方が、欠損と集計期間の不一致だけでなく、モデルをさらに調整すれば、例えば季節性を除去して長期トレンドの類似度だけみたい、といったようなニーズにも応えることができます。
最後に、私たちと一緒にデータサイエンスの力を使って社会を改善していきたい方はぜひこちらのページをご確認ください：