この記事は Qiita Advent Calendar 2025 - 時系列データの13日目の記事です。

はじめに

パターン抽出方法を元にした周期性解析を発案し特許を取得したので、実際どんな時に使えるのか、例として気温の時系列データを使って解析を行ってみたいと思います。
かなりざっくりとした使い方を掲載しています。
手法に関しては、数式を書いた記事と、図で説明した記事があるので読んでください。

忙しい人に向けて

新手法の利点
新手法は、一瞬の値の大きいノイズデータ(インパルスノイズ)が入っていても周期性を安定して計算する事が可能です。

課題

気象庁のデータ(ダウンロード先リンク)を使って、東京と釧路市阿寒湖畔の一日の気温差がどの程度違うのか定量評価したいと思います。
ただ、手法の説明をするだけなので深い意味はないです。

最低気温と最高気温の差の平均を比較する
新手法(bedcmm周期性解析を用いた特徴量)を比較する
フーリエ変換で振幅を比較する

正常な場合の計算

前処理

新手法が負の値が絡むと計算が大変なので、絶対温度にしています。

import numpy as np
import pandas as pd
from bedcmm import pattern as bedcmm # Githubで公開しているリポジトリがあります。
import matplotlib.pyplot as plt
import statsmodels.api as sm

# load data
df_tokyo = pd.read_csv('sample_data/data_東京_utf-8.csv',skiprows=4)
df_akankohan = pd.read_csv('sample_data/data_釧路市阿寒湖畔_utf-8.csv',skiprows=4)
# 欠損値補完
df_tokyo.iloc[:,1] = df_tokyo.iloc[:,1].interpolate(method='linear')
df_akankohan.iloc[:,1] = df_akankohan.iloc[:,1].interpolate(method='linear')
# 絶対温度化
K_cont = 273.15
np_temp_tokyo = np.array(df_tokyo.iloc[:,1]) + K_cont
np_temp_akankohan = np.array(df_akankohan.iloc[:,1]) + K_cont

気温の生データ

こんな感じのデータです。見るだけで、阿寒湖畔の方が東京より、気温差が大きそうですね。

1. 最低気温と最高気温の差の平均を比較する

普通に思いつくのはこれでしょう。最高気温と最低気温の差を計算する。
これを基準とします。

max_temp_day_tokyo = np.array([np.max(np_temp_tokyo[i*24:i*24+24]) for i in range(int(len(np_temp_tokyo)/24))])
min_temp_day_tokyo = np.array([np.min(np_temp_tokyo[i*24:i*24+24]) for i in range(int(len(np_temp_tokyo)/24))])
max_temp_day_akankohan = np.array([np.max(np_temp_akankohan[i*24:i*24+24]) for i in range(int(len(np_temp_akankohan)/24))])
min_temp_day_akankohan = np.array([np.min(np_temp_akankohan[i*24:i*24+24]) for i in range(int(len(np_temp_akankohan)/24))])
diff_temp_tokyo = max_temp_day_tokyo - min_temp_day_tokyo
diff_temp_akankohan = max_temp_day_akankohan - min_temp_day_akankohan
print(f"tokyo min max diff maen:{np.mean(diff_temp_tokyo)}")
print(f"akankohan min max diff mean:{np.mean(diff_temp_akankohan)}")
print(f"akankohan/tokyo diff ratio:{np.mean(diff_temp_akankohan)/np.mean(diff_temp_tokyo)}")

出力

tokyo min max diff maen:7.508743169398909
akankohan min max diff mean:10.229234972677595
akankohan/tokyo diff ratio:1.3623098755549083

阿寒湖畔と東京を比較すると、1.36倍気温差があるようです。

2. 新手法(bedcmm周期性解析を用いた特徴量)を比較する

新手法で周期性の分布を取得して、決まった範囲の周期性の標準偏差を計算する事によって、気温差を計算します。
ちなみに同じく周期性を出す手法のacfでも周期性の分布の標準偏差を見てみました。

max_lag = 97
min_lag = 12
periodicity_tokyo = bedcmm.periodicity_1d(np_temp_tokyo)
periodicity_akankohan = bedcmm.periodicity_1d(np_temp_akankohan)
acf_tokyo = sm.tsa.stattools.acf(np_temp_tokyo, nlags=max_lag)
acf_akankohan = sm.tsa.stattools.acf(np_temp_akankohan, nlags=max_lag)
print(f"tokyo periodicity std :{np.std(periodicity_tokyo[min_lag:max_lag])}")
print(f"akankohan periodicity std :{np.std(periodicity_akankohan[min_lag:max_lag])}")
print(f"akankohan/tokyo periodicyty ratio :{np.std(periodicity_akankohan[min_lag:max_lag])/np.std(periodicity_tokyo[min_lag:max_lag])}")
print(f"tokyo acf std :{np.std(acf_tokyo[min_lag:max_lag])}")
print(f"akankohan acf std :{np.std(acf_akankohan[min_lag:max_lag])}")
print(f"akankohan/tokyo acf ratio :{np.std(acf_tokyo[min_lag:max_lag])/np.std(acf_akankohan[min_lag:max_lag])}")

出力

tokyo periodicity std :0.2938155516187851
akankohan periodicity std :0.35027607630514623
akankohan/tokyo periodicyty ratio :1.192163159421924
tokyo acf std :0.04041222090078935
akankohan acf std :0.038326656597164846
akankohan/tokyo acf ratio :1.0544155031717215

新手法だと、阿寒湖畔と東京の比率が1.2倍程度になりました。少々小さいですが差は出ているようです。
acfでやると1.05倍と少ししかでないので感覚的に差を表現できていないと感じます。
ちなみに、bedcmm周期性解析とacfの結果も載せておきます。
24時間おきに周期性が高くなっている事が確認できます。この標準偏差を計算しています。

フーリエ変換で振幅を比較する

フーリエ変換で24[/hour]の振幅を比較してやる事もできます。
振幅とってくるあたりの書き方がちょっと書き方雑ですが。

tokyo_temp_fft = np.abs(np.fft.fft(np_temp_tokyo))[:int(len(np_temp_tokyo)/2)]/len(np_temp_tokyo)*2
akankohan_temp_fft = np.abs(np.fft.fft(np_temp_akankohan))[:int(len(np_temp_akankohan)/2)]/len(np_temp_akankohan)*2
temp_freq = np.fft.fftfreq(len(np_temp_tokyo))[:int(len(np_temp_akankohan)/2)]
day_freq_ind = np.where(temp_freq >= 1/24)[0][0]
print(f"day freq amp of tokyo:{tokyo_temp_fft[day_freq_ind]}")
print(f"day freq amp of akankohan:{akankohan_temp_fft[day_freq_ind]}")
print(f"akankohan/tokyo ratio of day freq amp:{akankohan_temp_fft[day_freq_ind]/tokyo_temp_fft[day_freq_ind]}")

出力

day freq amp of tokyo:2.6999921647291574
day freq amp of akankohan:3.498352907193192
akankohan/tokyo ratio of day freq amp:1.2956900219538674

阿寒湖畔と東京を比較すると、1.29倍となりました。1.のケースと同じような数字で納得感のある数字ではないでしょうか。
ちなみに、フーリエ変換の結果も載せておきます。

課題となるパターン

結局どれでもそれなりの数値が出てくるわけですが、新手法が何に強いかというと一瞬の値の大きいノイズデータ(インパルスノイズ)です。
前処理時に、インパルスノイズを付加して同じ計算を回します。

# load data
df_tokyo = pd.read_csv('sample_data/data_東京_utf-8.csv',skiprows=4)
df_akankohan = pd.read_csv('sample_data/data_釧路市阿寒湖畔_utf-8.csv',skiprows=4)
# 欠損値補完
df_tokyo.iloc[:,1] = df_tokyo.iloc[:,1].interpolate(method='linear')
df_akankohan.iloc[:,1] = df_akankohan.iloc[:,1].interpolate(method='linear')
# 絶対温度化
K_cont = 273.15
np_temp_tokyo = np.array(df_tokyo.iloc[:,1]) + K_cont
np_temp_akankohan = np.array(df_akankohan.iloc[:,1]) + K_cont

# インパルスノイズを付加
np_temp_tokyo[180] = 5000

グラフからとびぬけますが、こんな形になります。

1. 最低気温と最高気温の差の平均を比較する

出力

tokyo min max diff maen:20.397404371584702
akankohan min max diff mean:10.229234972677595
akankohan/tokyo diff ratio:0.5014968956995224

一つ大きな値が入るだけで真逆の結果になります。

2. 新手法(bedcmm周期性解析を用いた特徴量)を比較する

出力

tokyo periodicity std :0.2939404381995401
akankohan periodicity std :0.35027607630514623
akankohan/tokyo periodicyty ratio :1.1916566446273138

ほぼ値は変わりません。突発的なノイズに強い事がわかります。
周期性分布は以下のようになります。

acfは変わってしまいますが、bedcmm周期性解析はほとんど変わりません。

3. フーリエ変換で振幅を比較する

出力

day freq amp of tokyo:3.7184482157661316
day freq amp of akankohan:3.498352907193192
akankohan/tokyo ratio of day freq amp:0.9408099035399388

行けると思う人もいたかもしれませんが、いけません。
インパルスノイズは、全周波数帯域に平均してでるので、値が影響を受けてしまします。
ちなみに、周波数分布は以下のようになります。

まとめ

新手法では、一瞬の値の大きいノイズデータ(インパルスノイズ)がデータに入っても安定して、東京と阿寒湖畔の気温差の違いを計算する事ができました。
計算する際に最小値を利用しているので、インパルスノイズに強い処理となります。
ちょっと雑に書いてしまったので、興味や質問があればコメントください。

Github

新手法のGithub

明日の時系列データ Advent Calendarは @maskot1977 さんの時系列解析で読み解くCOVID-19波の伝播：差分系列と残差成分による比較分析です！