はじめに

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

本記事では、前回紹介したガウス過程近似法におけるハイパーパラメータ調整の効果を、実際のガウス過程との比較を通じて検証し、実証分析に携わる方々の参考となる知見を提供することを目的とします。

シミュレーションデータ

本記事では、こちらのRコードで生成したランダムウォークを対象に、ガウス過程およびその近似法による推定を行い、その性能を比較します：

set.seed(123)
gp_df <- rnorm(600) |>
  cumsum() |>
  tibble::tibble(
    y = _
  ) |>
  dplyr::mutate(
    x = dplyr::row_number(),
    x_std = (x - mean(x))/sd(x),
    y = y - mean(y)
  )

まずは形を確認しましょう：

gp_df |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = x, y = y)) + 
  ggplot2::labs(
    title = "データ"
  ) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

まさにランダムウォークらしい形ですね！もちろん、シードを変更すれば生成されるランダムウォークの形も変わり、それに伴って性能比較の結果も異なります。しかし、本記事では紙幅の都合上、シード間の比較は行いません。興味をお持ちの読者は、ぜひご自身でシードを変更し、結果の変化を確認してみてください。

Stanコード

ガウス過程モデルは以下のように実装します。本記事では、説明変数が時間に相当すると想定しているため、time_type というデータを渡しています。時系列データの場合、これは観測値の数 N と同一になります。ただし、同一の時点で複数の観測対象を持つパネルデータへの拡張を容易にするため、時間の数と観測値の数をあえて別の変数として定義しています。

gp_org.stan

data {
  int time_type;
  array[time_type] real time_seq;
  
  int N;
  
  array[N] real y;
}
parameters {
  real<lower=0> rho_quad;
  real<lower=0> alpha_quad;
  real<lower=0> sigma_quad;
  vector[time_type] eta_quad;
  
  real intercept;
  real<lower=0> sigma;
}
transformed parameters {
  vector[time_type] f_trend;
  {
    matrix[time_type, time_type] L_K;
    matrix[time_type, time_type] K = gp_exp_quad_cov(time_seq, alpha_quad, rho_quad);

    for (n in 1:time_type) {
      K[n, n] = K[n, n] + sigma_quad;
    }

    L_K = cholesky_decompose(K);
    f_trend = L_K * eta_quad - mean(L_K * eta_quad);
  }
}
model {
  rho_quad ~ inv_gamma(5, 5);
  alpha_quad ~ normal(0, 1);
  sigma_quad ~ normal(0, 1);
  eta_quad ~ normal(0, 1);
  
  intercept ~ normal(0, 10);
  sigma ~ inv_gamma(0.1, 0.1);
  
  y ~ normal(intercept + f_trend, sigma);
}

ガウス過程の近似法は、以下のように実装します。詳細については前回の記事（リンク）をご参照ください。概念的には、ガウス過程を複数の基底関数によって近似する手法です。本記事では、この近似に用いる関数の数 K をハイパーパラメータとして扱い、K の増加に伴って近似性能がどのように向上するかを検証します。

gp_approx.stan

functions {
  vector lambda(
    matrix S, real L, int j
  ){
    int dimension = dims(S)[2];
    vector[dimension] output;
    for (i in 1:dimension){
      output[i] = ((pi() * S[j, i])/(2 * L))^2;
    }
    return output;
  }
  
  real phi(
    vector x, matrix S, real L, int j
  ){
    int dimension = dims(S)[2];
    vector[dimension] output;
    for (i in 1:dimension){
      output[i] = sqrt(1/L) * sin(sqrt(lambda(S, L, j)[i]) * (x[i] + L));
    }
    return prod(output);
  }
  
  real s_inf(
    real alpha, real l, vector w
  ){
    int D = size(w);
    return alpha * (sqrt(2 * pi())^D) * (l ^ D) * exp(-0.5 * (l ^ 2) * w'*w);
  }
  
  real gaussian_process_approximation(
    vector x, matrix S, real L, real alpha, real l, vector betas
  ){
    int m = dims(S)[1];
    vector[m] output;
    for (i in 1:m){
      output[i] = (s_inf(alpha, l, sqrt(lambda(S, L, i))) ^ 0.5) * phi(x, S, L, i) * betas[i];
    }
    return sum(output);
  }
}
data {
  int approximation_type;
  int K;
  int approximation_type_pow_K;
  
  int time_type;
  
  int N;
  
  matrix[approximation_type_pow_K, K] S;
  
  array[N] vector[K] x;
  array[N] real y;
}
parameters {
  real<lower=0> alpha;
  real<lower=0> l;
  vector[approximation_type] betas;
  
  real intercept;
  real<lower=0> sigma;
}
transformed parameters {
  vector[time_type] f_trend;
  
  for (i in 1:time_type){
    f_trend[i] = gaussian_process_approximation(x[i], S, 1.5, alpha, l, betas);
  }
  
  f_trend = f_trend - mean(f_trend);
}
model {
  alpha ~ gamma(0.1, 0.1);
  l ~ gamma(0.1, 0.1);
  betas ~ normal(0, 1);
  
  intercept ~ normal(0, 10);
  sigma ~ inv_gamma(0.1, 0.1);
  
  y ~ normal(intercept + f_trend, sigma);
}

モデル推定

さて、ここではまず二つのモデルをコンパイルして、

m_org_init <- cmdstanr::cmdstan_model("gp_org.stan")

m_approx_init <- cmdstanr::cmdstan_model("gp_approx.stan")

ここでは、ガウス過程と、K = 5 および K = 20 の二種類の近似法による推定結果を比較します。

> m_org_summary <- m_org_init$variational(
     seed = 123,
     data = list(
         time_type = nrow(gp_df),
         time_seq = gp_df$x_std,
         
         N = nrow(gp_df),
         y = gp_df$y
     )
 )$
     summary()
------------------------------------------------------------ 
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.016977 seconds 
1000 transitions using 10 leapfrog steps per transition would take 169.77 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        -1738.102             1.000            1.000 
   200        -1675.198             0.519            1.000 
   300        -1664.798             0.348            0.038 
   400        -1659.580             0.262            0.038 
   500        -1661.346             0.210            0.006   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  21.5 seconds.

> m_approx_5_summary <- m_approx_init$variational(
     seed = 123,
     data = list(
         approximation_type = 5,
         K = 1,
         approximation_type_pow_K = 5,
         
         N = nrow(gp_df),
         
         S = matrix(1:5),
         
         time_type = nrow(gp_df),
         time_seq = gp_df$x_std,
         
         x = matrix(gp_df$x_std),
         y = gp_df$y
     )
 )$
     summary()
------------------------------------------------------------ 
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.002228 seconds 
1000 transitions using 10 leapfrog steps per transition would take 22.28 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        -7384.488             1.000            1.000 
   200        -3420.027             1.080            1.159 
   300        -2361.051             0.869            1.000 
   400        -2044.689             0.691            1.000 
   500        -1976.737             0.559            0.449 
   600        -1959.021             0.468            0.449 
   700        -1953.099             0.401            0.155 
   800        -1948.132             0.351            0.155 
   900        -1938.774             0.313            0.034 
  1000        -1916.906             0.283            0.034 
  1100        -1904.539             0.183            0.011 
  1200        -1886.912             0.068            0.009   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  11.3 seconds.

> m_approx_20_summary <- m_approx_init$variational(
     seed = 123,
     data = list(
         approximation_type = 20,
         K = 1,
         approximation_type_pow_K = 20,
         
         N = nrow(gp_df),
         
         S = matrix(1:20),
         
         time_type = nrow(gp_df),
         time_seq = gp_df$x_std,
         
         x = matrix(gp_df$x_std),
         y = gp_df$y
     )
 )$
     summary()
------------------------------------------------------------ 
EXPERIMENTAL ALGORITHM: 
  This procedure has not been thoroughly tested and may be unstable 
  or buggy. The interface is subject to change. 
------------------------------------------------------------ 
Rejecting initial value: 
  Gradient evaluated at the initial value is not finite. 
  Stan can't start sampling from this initial value. 
Gradient evaluation took 0.007288 seconds 
1000 transitions using 10 leapfrog steps per transition would take 72.88 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        -1885.455             1.000            1.000 
   200        -1768.022             0.533            1.000 
   300        -1768.368             0.356            0.066 
   400        -1765.720             0.267            0.066 
   500        -1766.652             0.214            0.001   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  26.4 seconds.

所要時間は、ガウス過程で 21.5 秒、K = 5 の近似法で 11.3 秒、K = 20 の近似法で 26.4 秒となりました。近似法なのに逆に所要時間が長くなるのは少し意外ですが、このような現象もあることを事実として認識しておくことが重要です。基底関数の数が増えると計算量も増えるため、場合によっては元のガウス過程よりも時間がかかってしまう可能性があります。

性能比較

さて、性能比較に入りましょう！

g_answer <- gp_df |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = x, y = y)) + 
  ggplot2::labs(
    title = "データ"
  ) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

g_org <- m_org_summary |>
  dplyr::filter(stringr::str_detect(variable, "f_trend")) |>
  dplyr::bind_cols(
    x = gp_df$x
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = x, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = x, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::labs(
    title = "ガウス過程"
  ) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")


g_approx_5 <- m_approx_5_summary |>
  dplyr::filter(stringr::str_detect(variable, "f_trend")) |>
  dplyr::bind_cols(
    x = gp_df$x
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = x, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = x, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::labs(
    title = "ガウス過程近似 K = 5"
  ) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

g_approx_20 <- m_approx_20_summary |>
  dplyr::filter(stringr::str_detect(variable, "f_trend")) |>
  dplyr::bind_cols(
    x = gp_df$x
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = x, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = x, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::labs(
    title = "ガウス過程近似 K = 20"
  ) + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")


gridExtra::grid.arrange(g_answer, g_org, g_approx_5, g_approx_20, ncol = 1)

あれ、、、、、、なんと、K = 20 の近似法のほうが、元のランダムウォークに近いですねwww

定性的な評価で恐縮ですが、性能としては K = 20 > ガウス過程 > K = 5 となるのではないかと思われます。つまり、近似法でありながら逆説的に、本物のガウス過程よりも性能が高くなる場合も実際には起こり得る現象です。

結論

いかがでしたか？本記事の性能比較を通じて、まず基底関数の数 K を増やすと性能が向上することを確認しました。ただし、逆説的に、基底関数の数を増やすことで元のガウス過程を上回る性能が出る場合もあります。今回紹介した近似法はまだ比較的新しい手法であるため、より精緻な数理的解析やシミュレーションを通じて、その強みと弱みを明確に可視化する必要があると感じました。

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