kernel = GPy.kern.Bias(1, variance=30) + \
         GPy.kern.Matern32(1, lengthscale=20, variance=0.1) + \
         GPy.kern.PeriodicMatern52(1, variance=1, lengthscale=2, period=1)
    
model = GPy.models.GPRegression(x, np.log(y), kernel)
model["sum.bias.variance"].constrain_fixed(30)
model["sum.periodic_Matern52.period"].constrain_fixed(1)

Test区間で推定値は実際よりも下にいきますが、これはバイアス項をかなり大きくしているためです。ただし、95％内には入っており、少なくともCI区間は感覚的にはいい感じと思ってます。

まだ腑に落ちていない点としては、モデルの決定方法。例えばconstrain_fixedせずにパラメータをoptimizeすると、bias.varianceはもっと小さくなります。それがそんなに悪いか？と言われるとちょっと気持ち悪い程度でうまく答えれられないのですが、理由がわからないのが気持ち悪い点でしょうか。

線形に増加するのだから、kern.Liniar()を使えば？とも思いますが、線形の仮定がめちゃめちゃつよく、外挿の際に分散が減少しないのが気持ち悪い。

他にも複数のkern.ExpQuadCosine(1)の重ね合わせで表現する際に、その数をAICで決定しようとすると、手元の環境では7つになりました。一方で、せいぜい3つぐらいから感覚的にはほとんどかわらないので、モデル選択が難しいですね。

RNN

やっとNNです。当然、RNNの実装例は多く、以下などがあります。

https://www.kaggle.com/code/xcwang21/flights-passenger-dataset-with-rnn-lstm-prophet

https://github.com/abhiupes01/Multi_Step_Time_Series/blob/master/RNN_flight_passanger.ipynb

が、1個目だと一見よく推定できているように見えますが、2個目で指摘しているように、Multi_stepで推定したときに全く推定できていません。直前のデータをそのまま返すだけになってますね。1個目でうまくできてそうに見えるのも、単純に入力を1ステップシフトしただけなんでしょう。

学習がうまくいっていないのはデータ数のせいなのか。そもそもAirplane passengersが高々120程度のデータ量で、RNNを使うのがどうかしてるのか。合成データをつくるのも一案ですが、とりあえず保留しています。

BNN(AR)

次にベイズ深層学習です。

ベイズ深層学習 (機械学習プロフェッショナルシリーズ)

作者:須山敦志
講談社

Amazon

実はこの本は3年前には購入していたのですが、なかなか難しく長い間放置されていました。

本来ならSimpleRNNか、LSTMを使いたかったのですが、先のとおり通常のRNNでうまく行っていないので、単純なAR（自己回帰モデル）を用います。 AR（auto regression）はその名のとおり、過去の自身の線形和で表されるので、単純なLinear層のみで表せます。例えば、過去12個まで見るなら、+bias項の計13個のパラメータと、RNNよりも圧倒的にパラメータ数が少なくて済みます。

BNNに話を戻して、本はとりあえず置いたままにして、結局他のブログ等を参考にBNNを理解しようとしました。そして今回やっと気づいたこと。

パラメータの事後分布は自由ではなく、基本的には正規分布とするよう。つまり変分ベイズ。パラメータ数は単純な2倍ですむ。
BNNは推定＋分散推定ではない => サンプルを出力するので、それから信用区間を算出（たぶんMCMCと同じような感じ）
パラメータの事前分布は一様ではなく、これも正規分布が基本のよう。

今回のような単純なNNなら分散の出力が可能だろうが、例えば活性化関数が入ったらそれだけで困難な気がする。なので思ったよりはBNN難しくないかも。

そもそもARのフィットがよかったのもあるが、信用区間にもそれっぽさがあります。

なお、BNNで検索すると、pyroとやらがたくさん出てきます。

Simple Bayesian Neural Network in Pyro | Kaggle

が、以下のように単純にpytorchのみで記述できそうで、なにが簡略化できるのかよくわかっていません。要勉強。

code example

class BayesianLinear(nn.Module):
    def __init__(self, in_features, out_features):
        super().__init__()
        self.mu = nn.Parameter(torch.Tensor(out_features, in_features))
        self.rho = nn.Parameter(torch.Tensor(out_features, in_features))
        self.bias_mu = nn.Parameter(torch.Tensor(out_features))
        self.bias_rho = nn.Parameter(torch.Tensor(out_features))
        self.reset_parameters()

    def reset_parameters(self):
        nn.init.kaiming_uniform_(self.mu)
        nn.init.constant_(self.rho, -3.0) # 初期値は小さめに設定

    def forward(self, x):
        sigma = torch.log1p(torch.exp(self.rho))
        weight = self.mu + sigma * torch.randn_like(self.mu) # 重みをサンプリング
        bias = self.bias_mu + torch.log1p(torch.exp(self.bias_rho)) * torch.randn_like(self.bias_mu) # biasをサンプリング
        return nn.functional.linear(x, weight, bias)

    def kl_loss(self,sig=3):
        mu_prior = torch.zeros_like(self.mu)
        sigma_prior = sig*torch.ones_like(self.mu)
        p = dist.Normal(mu_prior, sigma_prior)
        q = dist.Normal(self.mu, torch.log1p(torch.exp(self.rho)))
        kl_weight = dist.kl_divergence(q, p).sum()

        mu_prior = torch.zeros_like(self.bias_mu)
        sigma_prior = sig*torch.ones_like(self.bias_mu)
        p = dist.Normal(mu_prior, sigma_prior)
        q = dist.Normal(self.bias_mu, torch.log1p(torch.exp(self.bias_rho)))
        kl_bias = dist.kl_divergence(q, p).sum()
        return kl_weight + kl_bias

LLM4TS

zenn.dev

最近急激に身近になったLLMですが、これを時系列データに適用しようというもの。「いや、そんな簡単にできねーだろ」と思ったのですが、強化学習への適用例もあるので*1、なんでもありなのかな。

それでpromptのみでできるものもあるらしいので、試しに。

Prompt1

Here is some values for each month. Without using codes, Estimate the value over the 48 months following this.

37,  39,  44,  43,  40,  45,  49,  49,  45,  39,  34,  39,
38,  42,  47,  45,  41,  49,  56,  56,  52,  44,  38,  46,
...略
80,  77,  89,  89,  90, 105, 121, 115, 104,  91,  79,  92,
94,  92, 105, 104, 106, 124, 137, 135, 118, 102,  90, 102

Output

To estimate the values over the next 48 months based on the given data, I’ll assume the following:

Seasonality: There is a repeating seasonal pattern in the data, evident from the oscillations over months.
Trend: There is an upward trend over time, as the values generally increase across the years.
Here’s how the estimate can be constructed:

Decompose the data into trend and seasonality:

The trend increases roughly linearly, with some fluctuations. From the first year (37–49 range) to the last year (92–137 range), the trend grows by approximately 10–15 units per year.
The seasonality appears to repeat every 12 months with highs around mid-year and lows toward the end or start of the year.
Project the trend:

If the current year ends around a trend of ~120, the next year might start around 125 and increase to ~135 by year-end. This continues with an approximate 10–15 unit annual increase.
Add seasonal adjustments:

Seasonal variations are around ±10 to ±15 units from the trend.
Estimated Values (48 Months)
Year 1 (Months 1–12):
Trend grows from 125 to 135, adjusted seasonally:
126, 128, 134, 133, 130, 138, 144, 144, 139, 133, 125, 132

Year 2 (Months 13–24):
Trend grows from 135 to 145, adjusted seasonally:
136, 138, 144, 143, 140, 148, 154, 154, 149, 143, 135, 142

Year 3 (Months 25–36):
Trend grows from 145 to 155, adjusted seasonally:
146, 148, 154, 153, 150, 158, 164, 164, 159, 153, 145, 152

Year 4 (Months 37–48):
Trend grows from 155 to 165, adjusted seasonally:
156, 158, 164, 163, 160, 168, 174, 174, 169, 163, 155, 162

This projection reflects the combination of a consistent upward trend and repeating seasonal patterns over the next 48 months.

Prompt2

Estimate with additional 95% confidence intervals.

Output2

To estimate the 95% confidence intervals, we account for variability in the data around the trend and seasonality. Based on the provided values, here are the steps:

Steps for Adding Confidence Intervals
Calculate Variability:

The given data shows variability around the trend and seasonality. The deviation from the trend appears to range between ±5 and ±15 units.
Assume Normal Distribution:

A 95% confidence interval implies the range within approximately ±2 standard deviations from the mean. For simplicity, we can assume the variability (standard deviation) is about ±10 units, based on the observed fluctuations.
Confidence Interval Range:

For each estimated value, add and subtract 20 units (2 × 10 = 20) to establish the 95% confidence bounds.
Revised Estimates with Confidence Intervals
Year 1 (Months 1–12):

Month   Estimated Value Confidence Interval (±20)
1   126 106 to 146
2   128 108 to 148
3   134 114 to 154
4   133 113 to 153
5   130 110 to 150
6   138 118 to 158
7   144 124 to 164
8   144 124 to 164
9   139 119 to 159
10  133 113 to 153
11  125 105 to 145
12  132 112 to 152
Year 2 (Months 13–24):

Month   Estimated Value Confidence Interval (±20)
1   136 116 to 156
2   138 118 to 158
3   144 124 to 164
4   143 123 to 163
5   140 120 to 160
6   148 128 to 168
7   154 134 to 174
8   154 134 to 174
9   149 129 to 169
10  143 123 to 163
11  135 115 to 155
12  142 122 to 162
Year 3 (Months 25–36):

Month   Estimated Value Confidence Interval (±20)
1   146 126 to 166
2   148 128 to 168
3   154 134 to 174
4   153 133 to 173
5   150 130 to 170
6   158 138 to 178
7   164 144 to 184
8   164 144 to 184
9   159 139 to 179
10  153 133 to 173
11  145 125 to 165
12  152 132 to 172
Year 4 (Months 37–48):

Month   Estimated Value Confidence Interval (±20)
1   156 136 to 176
2   158 138 to 178
3   164 144 to 184
4   163 143 to 183
5   160 140 to 180
6   168 148 to 188
7   174 154 to 194
8   174 154 to 194
9   169 149 to 189
10  163 143 to 183
11  155 135 to 175
12  162 142 to 182
Summary
The estimated values include confidence intervals of ±20 units to account for seasonal and random variations, providing a more realistic range for the projections.

なお、当初Airplane passengersと具体的に書いていたところ、年は書いてないにもかかわらず、それを補ってアウトプットしだしたので、それらの文言を除いて、かつ具体的な数字からも類推されるかと思ったので元の数字を3で割った整数にして、入力しました。

また、まるで"Step by step"と入力したかのように、分析して結果ARIMAがよいと判断して、コードを書きだすことが多かったので、"Without code"を追加しています。とくに最近Reasoningがうまくなっているので、もしかするとGPT3.0などと古いほうがうまくいくのでは？と思っています。*2

何回か試して奇跡的にそれっぽい（つまり適当に）推定してきましたが、結果は推定については悪くはないがよくもない、推定区間は単純な±一定の適当な感じでした。

紹介例は自然言語も兼ねた適用例な一方、今回のは本当に単純なTimeseriesデータだったので、かえってうまくいかないのか。

もう少し実施例が広まったら、再度試してみたい。

FM4TS

Foundation Model for TimeSeries data.

arxiv.org

上では、GoogleのTimesFMを他の手法（ARIMAやGP)と比較しています。

Airplane passengersでは両者にはおよばないものの、他の多くの課題で1-3番目くらいのできで汎用性がある、とのこと。

前に示したとおり、RNNのうまくいかない原因の一つがデータ不足だとしたら、Pre-trainedのFoundation modelは解決策になるのも納得できるところ。

近いうちにやりたい。

まとめ

時系列データに対して複数の対応方法がある。従来法としてARIMAやGPがあり、いまでもそれらはトップクラスの性能だが、適切なモデル選択が必要。一方、自然言語につづき、Foundation modelができるとその適用が期待できるし、そもそもそういった自由度の高さがNNの魅力のひとつだと思う。

ついでに推定の際に、その推定範囲も示すことが実用上必要なことが多々ある。そのうえでBNNは強力な手法になりそう。Foundation modelも最後の出力層か＋αのみベイズ的にすれば変わらずできるのだろうか？

*1:https://qiita.com/pocokhc/items/44a8643e9570efb61a48

*2:が、とりあえずチャットベースですぐに確認できる、ChatGPT, Gemini, Claude, Mistralで試したところ、どれも明確な違いはなさそう

NezumiNoKuni’s blog

日々気になったことをメモしてます。ポケモン、データ分析、育児他を予定してます。

時系列データ分析（GP,RNN,BNN(AR),LLM4TS,FM4TSなどなど）

概要

データ

ARIMA

GP (Gaussian Process)

RNN

BNN(AR)

LLM4TS

FM4TS

まとめ