はじめに
こんにちは!この記事では、データ分析初心者の私が、Python
と時系列予測ライブラリProphet
を使って、野沢温泉(長野県)の過去の気象データから未来の積雪量を予測するプロジェクトに挑戦した記録をまとめます。単純な予測で終わらず、特徴量を追加してモデルを改善し、評価指標で客観的に評価するまでの一連のプロセスを体験しました。
積雪量を分析してみようと思った背景
今年の夏も非常に暑い日が続き、各所で最高気温を更新しました。地球温暖化という言葉を初めて耳にした約35年前から「Save The Snow、Save The Winter!」という合言葉の元、スキー場のクリーン活動を行い続けている昔からの馴染みに敬意を払い、少しでも今後の活動の後押しになれば!!そんな思いから、今回その活動拠点でウィンタースポーツの聖地として名高い 野沢温泉リゾート の積雪予測を行い、この活動の活性化に繋がればと思い至った次第です。また、個人的には、良い結果・悪い結果を全て受け入れて今後のウインタースポーツとの向き合い方、楽しみ方の一つの指標になれば良いなぁと考えております。
今回の目的
- 目的: 野沢温泉の地域活性化のため、各種イベントや宿泊施設の経営企画などに役立てられる指標として、2025~2026年の冬シーズンの野沢温泉における積雪量を予測する
- 使用言語: Python
-
主要ライブラリ:
pandas
,matplotlib
,seaborn
,prophet
-
データセット: 気象庁ホームページより野沢温泉日別気象データを収集 (詳細は次の通り)
- 期間:1984年1月~2025年7月まで
- 項目:年月・最深積雪(cm)・平均気温(℃)・日最高気温の平均(℃)・日最低気温の平均(℃)・最高気温(℃)・最低気温(℃)・日平均気温0℃未満日数(日)・日最高気温0℃未満日数(日)・日最低気温0℃未満日数(日)・降雪量合計(cm)・降雪量日合計最大(cm)・降雪量日合計3cm以上日数(日)・降水量の合計(mm)
- 環境: Windows + Jupyter Notebook
私的な分析ではありますが、この記事が以下のような方の参考になれば幸いです。
- データ分析や機械学習に興味があるけど、何から始めたらいいか分からない方
- Pythonやpandasの基本的な使い方を知りたい方
-
Prophet
を使った時系列予測の具体的な流れを知りたい方 - モデルの評価や改善のプロセスに興味がある方
それでは早速、始めてまいります。
ステップ1:データの読み込みと前処理
まずは、データをPythonで扱えるように準備します。
1-1. データの読み込みと確認
pandas
ライブラリを使ってCSVファイルを読み込み、head()
で最初の5行を表示して、データの中身を確認します。同時にデータ型も確認します。
# ライブラリをインポートして準備
import numpy as np
import pandas as pd
# CSVファイルを読み込みます
df = pd.read_csv('nozawa_data.csv')
# 最初の5行を表示して、データの中身を確認します
df.head()
# データの型と欠損値を確認します
df.info()
# もう少し詳しく欠損値を確認
df.isnull().sum()
# データの基本統計量を確認します
df.describe()
info()
でデータ型を確認したところ、日付(年月
)がただの文字列(object
)になっており、isnull().sum()
で欠損値を確認すると、積雪量や降雪量の列に135~136あることが分かりました。
【私の考察】
欠損値は、別途Excelで確認、夏期間のものだと判明したので
0
で補完の必要性がある。
年月列は、時系列分析のために正しい日付型(datetime
)に変換する必要がある。
df.describe()
で基本統計量を見ると面白いデータが並んでますね。日平均気温0℃未満日数と降雪量日合計3cm以上日、更に、日最低気温0℃未満日数がすべて共に31日となってますね。気象データが欲しくなる(またの機会に…)。
この方針に基づき、次の前処理を行います。
1-2. データ型の変換と欠損値の処理
pandasのto_datetimeを使って、日付の列を正しい日付型に変換します。
fillna(0)という命令を使って、すべての欠損値を0で埋めます。
# '年月'列を日付型に変換します
df['年月'] = pd.to_datetime(df['年月'], format='%b-%y')
# 欠損値を0で埋めます
df.fillna(0, inplace=True)
# 処理が正しく行われたか、再度確認します
df.info()
【ポイント】
-
pd.to_datetime()
で日付型に変換、format='%b-%y'
とすることで、現在の形式を指定し確実に変換することで、後に年や月を簡単に取り出したり、期間での絞り込みが楽になります。 -
fillna(0, inplace=True)
で、データフレーム内の全ての欠損値を0に置き換えています。
ステップ2:データの可視化と特徴量の検討
次に、データにどのような傾向があるか、グラフにして探ります。
2-1. ラインプロットで時間の経過と積雪量の変化を可視化
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import japanize_matplotlib
import seaborn as sns
# seabornで時間の経過と最深積雪、平均気温、最低気温の変化を確認
time = df[df['年月'] > '2014-12-01']
# 1つ目のグラフを描画(最深積雪(cm))
fig, ax = plt.subplots(figsize=(10, 4))
sns.lineplot(x='年月', y='最深積雪(cm)', data=time, color='blue')
plt.title('最深積雪(cm)/過去10年')
# 11月と4月だけ日付を出します
ax.xaxis.set_major_locator(mdates.MonthLocator(bymonth=[4, 11]))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax.tick_params(axis='x', rotation=45)
plt.ylabel('最深積雪')
plt.show()
# 2つ目のグラフを描画(平均気温(℃))
plt.figure(figsize=(10, 4))
sns.barplot(x='年月', y='平均気温(℃)', data=time, color='red')
plt.title('平均気温(℃)/年')
# 11月と4月だけ日付を出します
xticks = []
xlabels = []
for i, d in enumerate(time['年月']):
if d.month in [4, 11]: # 4月と11月のみ
xticks.append(i)
xlabels.append(d.strftime('%Y-%m'))
plt.xticks(xticks, xlabels, rotation=45)
plt.ylabel('平均気温')
plt.show()
# 3つ目のグラフを描画(最低気温(℃))
plt.figure(figsize=(10, 4))
sns.barplot(x='年月', y='最低気温(℃)', data=time, color='blue')
plt.title('最低気温(℃)/年')
# 11月と4月だけ日付を出します
xticks = []
xlabels = []
for i, d in enumerate(time['年月']):
if d.month in [4, 11]: # 4月と11月のみ
xticks.append(i)
xlabels.append(d.strftime('%Y-%m'))
plt.xticks(xticks, xlabels, rotation=45)
plt.ylabel('最低気温')
plt.show()
# 4つ目のグラフを描画(降雪量合計(cm))
plt.figure(figsize=(12, 4))
sns.barplot(x='年月', y='降雪量合計(cm)', data=time, color='blue')
plt.title('降雪量合計(cm)/年')
# 11月と4月だけ日付を出します
xticks = []
xlabels = []
for i, d in enumerate(time['年月']):
if d.month in [4, 11]: # 4月と11月のみ
xticks.append(i)
xlabels.append(d.strftime('%Y-%m'))
plt.xticks(xticks, xlabels, rotation=45)
plt.ylabel('降雪量合計(cm)')
plt.show()
その年によって積雪量にもバラツキがあり、予測が難しそうですね。
グラフの山と谷はそれぞれ、冬季と夏季だということが分かります。
【私の考察】
シーズンによっての降雪量の変化が気になる!積雪に大きな影響を及ぼすであろう特徴量なので冬季シーズンに限定した方が良いかもしれない。
2-2. barplotで冬季の降雪量(cm)を可視化
ということで、冬季に限定(11月~4月間)して降雪量合計(cm)の合計を可視化してみたいと思います。
# 年と月を抽出
df['year'] = df['年月'].dt.year
df['month'] = df['年月'].dt.month
# 積雪シーズンを定義する'season_year'列を作成
# 11月と12月は同じ年、1月から4月は前年の年として扱う
df['season_year'] = df['year'].where(df['month'] >= 11, df['year'] - 1)
# 11月~4月の期間のデータのみを抽出
snow_season_df = df[df['month'].isin([11, 12, 1, 2, 3, 4])]
# ★ ここで2014年以降のデータのみに絞り込みます
filtered_df = snow_season_df[snow_season_df['season_year'] > 2004]
# 'season_year'ごとに積雪量を合計
yearly_total_snowfall = filtered_df.groupby('season_year')['降雪量合計(cm)'].sum().reset_index()
# バーグラフを作成
plt.figure(figsize=(10, 6))
bars = plt.bar(yearly_total_snowfall['season_year'], yearly_total_snowfall['降雪量合計(cm)'], color='skyblue')
# 各バーの上に積雪合計量を表示
for bar in bars:
yval = bar.get_height()
plt.text(bar.get_x() + bar.get_width()/2.0, yval, int(yval), va='bottom', ha='center', fontsize=12)
# グラフのタイトルとラベルを設定
plt.title('年ごとの降雪量合計(cm)(11月~翌年4月)', fontsize=16)
plt.xlabel('年', fontsize=12)
plt.ylabel('降雪量合計(cm)', fontsize=12)
# X軸の目盛りをシーズン年に設定
plt.xticks(yearly_total_snowfall['season_year'])
# グリッド線を追加
plt.grid(axis='y', linestyle='--', alpha=0.7)
年ごとの降雪量にも差がありますね。数字も加えるとその差が良く分かります。この降雪量を積雪につなげるのがズバリ「気温」ですね。気温が低いほど積雪に繋がるのでしょう。
【私の考察】
気温と降雪の関係も可視化してみよう。気温は低いときに、積雪に繋がるだろう。
2-3. 気温と降雪量のサブプロットの作成
関連性のありそうな、気温と降雪量の関係をグラフで可視化
# 1. 2014年以降のデータに絞り込みます
df_last10 = df[df['年月'].dt.year >= 2020].copy()
# 2. グラフの描画領域を準備します
fig, ax1 = plt.subplots(figsize=(15, 7))
# 3. 1つ目の軸(左側)に積雪量のグラフを描画します
ax1.bar(df_last10['年月'], df_last10['最深積雪(cm)'], color='blue', width=10, label='最深積雪(cm)')
ax1.set_ylabel('最深積雪(cm)')
ax1.set_xlabel('年月')
# 4. 2つ目の軸(右側)を作成します
ax2 = ax1.twinx()
# 5. 2つ目の軸に平均気温のグラフを描画します
ax2.plot(df_last10['年月'], df_last10['平均気温(℃)'], color='r', marker='o', linestyle='-', label='平均気温(℃)')
ax2.set_ylabel('平均気温 (℃)')
# 6. グラフの見た目を整えます
plt.title('過去10年間の積雪量と平均気温(野沢温泉)')
fig.legend(loc='upper left', bbox_to_anchor=(0.1, 0.9))
plt.grid(True)
# x軸のフォーマットを年月が見やすいように調整
ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
# 7. グラフを表示します
plt.show()
やはり負の相関関係ですね。(目的が積雪量の予測ということで、冬季の変化を確認しやすいように過去5年に絞った表を描画させています。
【私の考察】
平均気温が低い12月~4月にかけての「最深積雪」が、主な予測の目的変数となる。
反対に目標変数に対して積雪の無い時期が(5月~10月)と長く実績が図れないため時期を絞る必要がありそうだ。
他の特徴量との関係も気になるので、ヒートマップで相関関係の全体像を把握したい。
2-4. 相関ヒートマップの作成
どの特徴量が積雪量と関係が深いかを探るため、ヒートマップを作成しました。
# 相関係数を計算します
correlation_matrix = df.corr()
# seabornを使ってヒートマップを描画します
plt.figure(figsize=(12, 10))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt='.2f')
plt.title('特徴量の相関ヒートマップ')
plt.show()
「平均気温(℃)」を予測モデル構築時の一番の有力特徴量だと思い込んでいましたが、「日最高気温の平均(℃)」や 「最高気温(℃)」の方が、より強い負の相関で、「日最低気温0℃未満日数(日)」、「降雪量日合計3cm以上日数(日)」 が、強い正の相関があることが分りました。
【私の考察】
多重共線性の高い特徴量が多い場合、モデルに不要なノイズが加わり、予測精度が低下する可能性があります。
そこで今回は、まずはシンプルに「日最高気温の平均(℃)」と「降雪量日合計3cm以上日数(日)」の2つを有力な特徴量として選び、予測モデルを構築します。更に「気温関連」と「降雪関連」に分類した中から特徴量の追加や削除をして、モデルの精度を確認するのが良さそう。
【ポイント】
- ヒートマップを見ると、「最深積雪(cm)」は「日最高気温の平均(℃)」と強い負の相関(色が濃い青)があることが一目で分かりました。これは「日最高気温が低いほど積雪が多くなる」という結果が可視化により明確になったことで、直感(平均気温(℃))とは違い、予測モデル構築の重要な手がかりとなります。
- また、似たような特徴量が多く、多重共線性(ノイズの原因)が考えられるので、予測モデル構築の使用する特徴量の選定には注意を払う必要があります。
ステップ3:Prophetによる予測モデル構築 (V1)
いよいよProphet
を使ってモデルを構築します。
上記の考察を基にデータフレームから有力特徴量を絞り込みモデル構築の準備を行います。
3-1. データの準備とモデル学習
Prophetが扱える形式(日付列ds
、予測値列y
)にデータフレームを整え、学習用(train
)と評価用(test
)に分割します。
-
学習(train)期間(
ds
): 1984-01 ~ 2023-12 -
評価(test)期間(
ds
): 2024-01 ~ 2025-03 - 学習期間: 11月 ~ 04月に限定(※降雪の無い時期は除外)⇒考察より
-
予測対象(
y
):最深積雪(cm)
-
特徴量(regressor):
日最高気温の平均(℃)
,降雪量日合計3cm以上日数(日)
⇒考察より
from prophet import Prophet
train
/test
の期間と使用する特徴量を抽出
# 1. Prophet用に列名を変更し、必要な列を抽出
df_prophet = df[['年月', '最深積雪(cm)', '日最高気温の平均(℃)', '降雪量日合計3cm以上日数(日)']].copy()
df_prophet.rename(columns={'年月': 'ds', '最深積雪(cm)': 'y'}, inplace=True)
# 2. データを冬季 (11月, 12月, 1月, 2月, 3月, 4月) のみに絞り込む
df_winter = df_prophet[df_prophet['ds'].dt.month.isin([11, 12, 1, 2, 3, 4])]
# 3. データを学習用と評価用に分割
train = df_winter[df_winter['ds'] <= '2023-12-31']
test = df_winter[df_winter['ds'] > '2023-12-31']
指定の通りtrain
/test
に分割されたか確認
# 4. 念のため分割された中身を確認
print("--学習データ(最後の5行)")
print(train.tail()) # trainの終わり部分を確認
print("\n --評価用データ(最後の5行)") # >綺麗に分かれているか確認
print(test.head()) # testの開始部分を確認
対象の月に絞られていて、綺麗に分かれていることが確認出来た
# 3. Prophetモデルの初期化と学習
model = Prophet() # 特別な設定は無くFull Autoで行こうと思います。
# 特徴量を追加
# add_regressorという命令で、予測に使う追加の特徴量を設定します
model.add_regressor('日最高気温の平均(℃)')
model.add_regressor('降雪量日合計3cm以上日数(日)')
# trainに対し作成したモデルをfitさせて学習を実施します。
model.fit(train)
3-2. モデルの評価 (グラフ)
学習済みモデルで評価期間の予測を行い、実績値と予測値をグラフで比較します。
ここでは、データのふるまいを再現できているかを可視化するグラフの表示を行います。
regs = ['日最高気温の平均(℃)', '降雪量日合計3cm以上日数(日)']
# 予測を実行
forecast_test_only = model.predict(test[['ds'] + regs]).sort_values('ds')
# ☆★☆
# データセット(実績)のふるまいを捉えられているか可視化用
# 履歴+未来を縦に結合して predict(プロット用)
future_all = pd.concat([train[['ds'] + regs], test[['ds'] + regs]], ignore_index=True)
future_all = future_all.drop_duplicates('ds').sort_values('ds')
# 学習データも入れた内容で予測を実行
forecast_all = model.predict(future_all).sort_values('ds')
# 予測が変わっていないか念のため検証(最大誤差を表示)※カンニング対策
fa_test_part = forecast_all[forecast_all['ds'].isin(test['ds'])].sort_values('ds')
max_abs_diff = np.max(np.abs(fa_test_part['yhat'].to_numpy() - forecast_test_only['yhat'].to_numpy()))
# カンニングになっていないかを確認(誤差を確認)
print("max |Δyhat| on test:", max_abs_diff)
実際のデータにどれくらい適合しているかを可視化しているため、実績を含んだ検証用のデータに対しても予測を実施しています。カンニングになっていないかも確認しています。
マイナスの予測値を0に補正
# マイナスの予測値を0に補正
# forecast_test_only
forecast_test_only['yhat'] = forecast_test_only['yhat'].clip(lower=0)
forecast_test_only['yhat_lower'] = forecast_test_only['yhat_lower'].clip(lower=0)
forecast_test_only['yhat_upper'] = forecast_test_only['yhat_upper'].clip(lower=0)
# ☆★☆
# forecast_all
forecast_all['yhat'] = forecast_all['yhat'].clip(lower=0)
forecast_all['yhat_lower'] = forecast_all['yhat_lower'].clip(lower=0)
forecast_all['yhat_upper'] = forecast_all['yhat_upper'].clip(lower=0)
範囲を指定してグラフで確認
# 予測グラフ描画
start = pd.Timestamp('2015-07-31')
end = pd.Timestamp('2025-07-31')
fig_forecast = model.plot(forecast_all)
plt.xlim(start, end)
plt.show()
【ポイント】
-
clip(lower=0)
: 積雪0cm時の予測でマイナスになる事も有るのでlower=0
としてマイナス(ー)の予測値を0に留めておきます。
予測モデル比較用のグラフも作成しておきます。
# 予測結果と実測値を1つのグラフに描画して比較
fig, ax = plt.subplots(figsize=(15, 7))
# 実際のプロット(testデータのy)
ax.plot(test['ds'], test['y'], 'r.-', label='実績値(最深積雪)')
# 予測値のプロット(forecast_testのyhat)
# yhatが予測値、yhat_lowerとyhat_upperが予測の信頼区間(予測の振れ幅)
ax.plot(forecast_test['ds'], forecast_test['yhat'], 'b.-', label='予測値')
ax.fill_between(forecast_test['ds'], forecast_test['yhat_lower'],
forecast_test['yhat_upper'], color='blue', alpha=0.2, label='予測の信頼区間')
#グラフの見た目を整える
ax.set_title('積雪量の予測結果と実績値の比較(2024/01 - 2024/03)')
ax.set_xlabel('年月')
ax.set_ylabel('最深積雪(cm)')
ax.legend()
plt.grid(True)
plt.show()
【私の考察】
特徴量の情報がもう少しあると良いモデルが作れるかもしれない。特徴量を徐々に追加して予測モデルの変動を確認する必要がある。
【ポイント】
-
yhat
,yhat_lower
,yhat_upper
に負の値が含まれるため、clip(lower=0)
で0未満を補正。特にシーズン前半の11月と、後半の4月に積雪が0の年もあるため、この補正は実用上問題ないと判断できます。 -
ax.fill_between(...)
:yhat_lower
とyhat_upper
の間を塗りつぶす命令です。これにより、予測の振れ幅が視覚的に分か りやすくなります。alpha=0.2
は、色の透明度を20%に設定するオプションです。
ステップ4:モデルの評価と改善
グラフでの評価を踏まえ、特徴量を追加してモデルの改善を試してみる。
4-1. 評価指標(MAE, RMSE)の導入
改善度合いを客観的に測るため、数値による評価指標を導入しました。
- MAE (平均絶対誤差): 予測が平均してどれくらいズレているか (cm)
- RMSE (二乗平均平方根誤差): 大きなズレをより重視する指標 (cm)
from sklearn.metrics import mean_absolute_error, mean_squared_error
# 1. 最初のモデルの評価指標を計算
mae_v1 = mean_absolute_error(test['y'], forecast_test['yhat'])
rmse_v1 = np.sqrt(mean_squared_error(test['y'], forecast_test['yhat']))
print("--- 最初のモデルの評価指標 ---")
print(f"MAE (平均絶対誤差): {mae_v1:.2f} cm")
print(f"RMSE (二乗平均平方根誤差): {rmse_v1:.2f} cm")
print("-" * 30)
どちらも小さいほど良いモデル。(最後にまとめて比較します。)
【ポイント】
- どちらも値が小さいほど良いモデルといえます。
- 夏期(降雪0)を含めないことで、誤差が少ない期間によるスコアの“かさ上げ”を防ぎ、より厳しい評価になっています。
4-2. 予測結果のテーブルも確認
# どのようなデータが予測されているか、先頭行を確認
print("--予測結果データ(先頭5行)--")
print(forecast_test[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].head())
【ポイント】
-
model.predict(test)
:学習済みの予測モデルに、未来の期間と、その期間の特徴量(気温など)が含まれたtest
データフレー ムを渡すことで、予測を計算させます。(ステップ3-2参照) -
forecast_test:predict
メソッドが返す結果は、新しいデータフレームです。この中にはたくさんの情報が含まれていま すが、特に重要なのは以下の3つです。-
ds
:予測した日付 -
yhat
:予測値 (y-hatと読みます) -
yhat_lower
,yhat_upper
:予測の信頼区間。「この範囲に80%(デフォルト設定)の確率で収まるだろう」という、予測の上下限の範囲を表します。
-
ステップ5':特徴量を追加してモデルを改善する
改善方針 今度は「日最低気温0℃未満日数(日)」を特徴量として加え、再度モデルの精度向上を試みます。
train
/test
の期間と使用する特徴量を抽出
# 1. Prophetで使うためのデータフレームを再作成(新しい特徴量を追加)
df_prophet_v2 = df[['年月', '最深積雪(cm)','日最高気温の平均(℃)', '降雪量日合計3cm以上日数(日)',
'日最低気温0℃未満日数(日)']].copy()
df_prophet_v2.rename(columns={'年月': 'ds', '最深積雪(cm)': 'y'}, inplace=True)
train
/test
の期間を指定
# 2. データを冬季 (11月, 12月, 1月, 2月, 3月, 4月) のみに絞り込む
df_winter_v2 = df_prophet_v2[df_prophet_v2['ds'].dt.month.isin([11, 12, 1, 2, 3, 4])]
# 3. データを学習用(train_v3)と評価用(test_v3)に再分割
train_v2 = df_winter_v2[df_winter_v2['ds'] <= '2023-12-31']
test_v2 = df_winter_v2[df_winter_v2['ds'] > '2023-12-31']
指定の通りtrain
/test
に分割されたか確認
# 4. 念のため分割された中身を確認
print("--学習データ(最後の5行)")
print(train_v2.tail()) # trainの終わり部分を確認
print("\n --評価用データ(最後の5行)") # >綺麗に分かれているか確認
print(test_v2.head()) # testの開始部分を確認
新たな特徴量で予測モデルを作成
# 4. 新しいProphetモデルを初期化
model_v2 = Prophet()
追加で使用する特徴量は「日最低気温0℃未満日数(日)
」
# 5. 3つの特徴量を追加
model_v2.add_regressor('日最高気温の平均(℃)')
model_v2.add_regressor('降雪量日合計3cm以上日数(日)')
model_v2.add_regressor('日最低気温0℃未満日数(日)') # 今度はこちらの特徴量を追加してみる
予測を実施
model_v2.fit(train_v2)
検証用データを作成
regs_v2 = ['日最高気温の平均(℃)', '降雪量日合計3cm以上日数(日)','日最低気温0℃未満日数(日)']
# 予測を実行
forecast_test_only_v2 = model_v2.predict(test_v2[['ds'] + regs_v2]).sort_values('ds')
# ☆★☆
# データセット(実績)のふるまいを捉えられているか可視化用
# 履歴+未来を縦に結合して predict(プロット用)
future_all_v2 = pd.concat([train_v2[['ds'] + regs_v2], test_v2[['ds'] + regs_v2]], ignore_index=True)
future_all_v2 = future_all_v2.drop_duplicates('ds').sort_values('ds')
# 学習データも入れた内容で予測を実行
forecast_all_v2 = model_v2.predict(future_all_v2).sort_values('ds')
# 予測が変わっていないか念のため検証(最大誤差を表示) ※カンニング対策
fa_test_part_v2 = forecast_all_v2[forecast_all_v2['ds'].isin(test_v2['ds'])].sort_values('ds')
max_abs_diff_v2 = np.max(np.abs(fa_test_part_v2['yhat'].to_numpy() - forecast_test_only_v2['yhat'].to_numpy()))
# カンニングになっていないかを確認(誤差を確認)
print("max |Δyhat| on test:", max_abs_diff_v2)
マイナスの予測値を0に補正
# マイナスの予測値を0に補正
# forecast_test_only
forecast_test_only_v2['yhat'] = forecast_test_only_v2['yhat'].clip(lower=0)
forecast_test_only_v2['yhat_lower'] = forecast_test_only_v2['yhat_lower'].clip(lower=0)
forecast_test_only_v2['yhat_upper'] = forecast_test_only_v2['yhat_upper'].clip(lower=0)
# forecast_all
forecast_all_v2['yhat'] = forecast_all_v2['yhat'].clip(lower=0)
forecast_all_v2['yhat_lower'] = forecast_all_v2['yhat_lower'].clip(lower=0)
forecast_all_v2['yhat_upper'] = forecast_all_v2['yhat_upper'].clip(lower=0)
期間を指定して可視化(ふるまいを確認)
# 予測グラフ描画
start = pd.Timestamp('2015-07-31')
end = pd.Timestamp('2025-07-31')
fig_forecast = model.plot(forecast_all_v2)
plt.xlim(start, end)
plt.show()
予測モデル比較用のグラフ(前のモデルと比較)
# 7. 改善されたかグラフで比較
fig, ax = plt.subplots(figsize=(15, 7))
# 実績値
ax.plot(test_v2['ds'], test_v2['y'], 'r.-', label='実績値')
# 最初のモデルの予測値 (比較のため)
ax.plot(forecast_test_only['ds'], forecast_test_only['yhat'],
'g:', label='予測値 (降雪量日合計3cm以上)', alpha=0.7)
# 今回のモデルの予測値
ax.plot(forecast_test_only_v2['ds'],
forecast_test_only_v2['yhat'],
'b.-', label='予測値 (降雪量日合計3cm以上 + 日最低気温0℃未満日数モデル)')
ax.fill_between(forecast_test_only_v2['ds'],
forecast_test_only_v2['yhat_lower'],
forecast_test_only_v2['yhat_upper'],
color='blue', alpha=0.2, label='降雪量日合計3cm以上 + 日最低気温0℃未満日数モデルの信頼区間')
# グラフの見た目を整える
ax.set_title('新旧モデルの予測結果比較 (2024/01 - 2025/03)')
ax.set_xlabel('年月')
ax.set_ylabel('最深積雪 (cm)')
ax.legend()
plt.grid(True)
plt.show()
【私の考察】
上限と下限は僅かに良くなった気がするが、グラフだけでは判断できないので、後でMAE値とRMSE値のスコアで数値的に判断することとします。
ステップ5'':特徴量を「降水量」に変更してモデルを改善する
改善方針 最後は更に「降水量の合計(mm)」を特徴量として加え、再度モデルの精度向上を試みます。
train
/test
の期間と使用する特徴量を抽出
# 1. Prophetで使うためのデータフレームを再作成(新しい特徴量を追加)
df_prophet_v3 = df[['年月', '最深積雪(cm)',
'日最高気温の平均(℃)', '降雪量日合計3cm以上日数(日)',
'降水量の合計(mm)']].copy()
df_prophet_v3.rename(columns={'年月': 'ds', '最深積雪(cm)': 'y'}, inplace=True)
train
/test
の期間を指定
# 2. データを冬季 (11月, 12月, 1月, 2月, 3月, 4月) のみに絞り込む
df_winter_v3 = df_prophet_v3[df_prophet_v3['ds'].dt.month.isin([11, 12, 1, 2, 3, 4])]
# 3. データを学習用(train_v3)と評価用(test_v3)に再分割
train_v3 = df_winter_v3[df_winter_v3['ds'] <= '2023-12-31']
test_v3 = df_winter_v3[df_winter_v3['ds'] > '2023-12-31']
指定の通りtrain
/test
に分割されたか確認
# 4. 念のため分割された中身を確認
print("--学習データ(最後の5行)")
print(train_v3.tail()) # trainの終わり部分を確認
print("\n --評価用データ(最後の5行)") # >綺麗に分かれているか確認
print(test_v3.head()) # testの開始部分を確認
新たな予測モデルを作成
# 4. 新しいProphetモデルを初期化
model_v3 = Prophet()
使用する特徴量を「降水量の合計(mm)
」に変更
# 5. 3つの特徴量を追加
model_v3.add_regressor('日最高気温の平均(℃)')
model_v3.add_regressor('降雪量日合計3cm以上日数(日)')
model_v3.add_regressor('降水量の合計(mm)') # 今度はこちらの特徴量
予測を実行
model_v3.fit(train_v3)
検証用データを作成
regs_v3 = ['日最高気温の平均(℃)', '降雪量日合計3cm以上日数(日)','降水量の合計(mm)']
# 予測を実行
forecast_test_only_v3 = model_v3.predict(test_v3[['ds'] + regs_v3]).sort_values('ds')
# ☆★☆
# データセット(実績)のふるまいを捉えられているか可視化用
future_all_v3 = pd.concat([train_v3[['ds'] + regs_v3], test_v3[['ds'] + regs_v3]], ignore_index=True)
future_all_v3 = future_all_v3.drop_duplicates('ds').sort_values('ds')
# 学習データも入れた内容で予測を実行
forecast_all_v3 = model_v3.predict(future_all_v3).sort_values('ds')
# 予測が変わっていないか念のため検証(最大誤差を表示) ※カンニング対悪
fa_test_part_v3 = forecast_all_v3[forecast_all_v3['ds'].isin(test_v3['ds'])].sort_values('ds')
max_abs_diff_v3 = np.max(np.abs(fa_test_part_v3['yhat'].to_numpy() - forecast_test_only_v3['yhat'].to_numpy()))
# カンニングになっていないかを確認(誤差を確認)
print("max |Δyhat| on test:", max_abs_diff_v3)
マイナスの予測値を0に補正
# マイナスの予測値を0に補正
# forecast_test_only
forecast_test_only_v3['yhat'] = forecast_test_only_v3['yhat'].clip(lower=0)
forecast_test_only_v3['yhat_lower'] = forecast_test_only_v3['yhat_lower'].clip(lower=0)
forecast_test_only_v3['yhat_upper'] = forecast_test_only_v3['yhat_upper'].clip(lower=0)
# forecast_all
forecast_all_v3['yhat'] = forecast_all_v3['yhat'].clip(lower=0)
forecast_all_v3['yhat_lower'] = forecast_all_v3['yhat_lower'].clip(lower=0)
forecast_all_v3['yhat_upper'] = forecast_all_v3['yhat_upper'].clip(lower=0)
期間を指定して可視化(ふるまいを確認)
# 予測グラフ描画
start = pd.Timestamp('2015-07-31')
end = pd.Timestamp('2025-07-31')
fig_forecast = model.plot(forecast_all_v3)
plt.xlim(start, end)
plt.show()
予測モデル比較用のグラフ(歴代のモデルと比較)
# 7. 改善されたかグラフで比較
fig, ax = plt.subplots(figsize=(15, 7))
# 実績値
ax.plot(test_v3['ds'], test_v3['y'], 'r.-', label='実績値')
# 最初のモデルの予測値 (比較のため)
ax.plot(forecast_test_only['ds'], forecast_test_only['yhat'],
'm:', label='予測値 (降雪量日合計3cm以上)', alpha=0.7)
# V2モデルの予測値 (比較のため)
ax.plot(forecast_test_only_v2['ds'], forecast_test_only_v2['yhat'],
'g:', label='予測値 (降雪量日合計3cm以上 + 日最低気温0℃未満日数)', alpha=0.7)
# 今回のモデルの予測値
ax.plot(forecast_test_only_v3['ds'],
forecast_test_only_v3['yhat'], 'b.-', label='予測値 (降雪量日合計3cm以上 + 降水量モデル)')
ax.fill_between(forecast_test_only_v3['ds'],
forecast_test_only_v3['yhat_lower'],
forecast_test_only_v3['yhat_upper'],
color='blue',
alpha=0.2,
label='降水量モデルの信頼区間'
)
# グラフの見た目を整える
ax.set_title('新旧モデルの予測結果比較 (2024/01 - 2025/03)')
ax.set_xlabel('年月')
ax.set_ylabel('最深積雪 (cm)')
ax.legend()
plt.grid(True)
plt.show()
歴代の予測モデルと比較がグラフ化されました。視覚的に捉えることでモデルの判定がし易くなります。
【私の考察】
視覚的には2番目のモデル「日最低気温0℃未満日数」を追加したモデルが全体的にフィットしている様な気がしますが、判断が難しいのでスコアを確認しよう。
ステップ6''':評価指標でモデルを数値的に比較する
scikit-learnのライブラリの機能を使って、最初のモデルと、その後行った2つのモデル(V2とV3)のMAEとRMSEをそれぞれ計算し、比較します。
# 1. 最初のモデルの評価指標を計算
mae_v1 = mean_absolute_error(test['y'], forecast_test_only['yhat'])
rmse_v1 = np.sqrt(mean_squared_error(test['y'], forecast_test_only['yhat']))
print("--- 最初のモデルの評価指標 ---")
print(f"MAE (平均絶対誤差): {mae_v1:.2f} cm")
print(f"RMSE (二乗平均平方根誤差): {rmse_v1:.2f} cm")
print("-" * 30)
# 2. 日最低気温0℃未満日数(日)を加えたモデルの評価指標を計算
mae_v2 = mean_absolute_error(test_v2['y'], forecast_test_only_v2['yhat'])
rmse_v2 = np.sqrt(mean_squared_error(test_v2['y'], forecast_test_only_v2['yhat']))
print("--- 日最低気温0℃未満日数(日)モデルの評価指標 ---")
print(f"MAE (平均絶対誤差): {mae_v2:.2f} cm")
print(f"RMSE (二乗平均平方根誤差): {rmse_v2:.2f} cm")
print("-" * 30)
# 3. 降水量を加えたモデルの評価指標を計算
mae_v3 = mean_absolute_error(test_v3['y'], forecast_test_only_v3['yhat'])
rmse_v3 = np.sqrt(mean_squared_error(test_v3['y'], forecast_test_only_v3['yhat']))
print("--- 降水量モデルの評価指標 ---")
print(f"MAE (平均絶対誤差): {mae_v3:.2f} cm")
print(f"RMSE (二乗平均平方根誤差): {rmse_v3:.2f} cm")
print("-" * 30)
数値的に判断しやすくなりました。全体的には2番目のモデルが良さそうです。
【私の考察】
モデルの比較グラフからの考察どおり2番目の「日最低気温0℃未満日数(日)」が一番よいモデルと決定!
このモデルを使って未来予測をしてみよう。
ポイント
-
from sklearn.metrics import ...
:scikit-learn
ライブラリのmetrics
(評価指標)モジュールから、必要な関数をインポートしています。 -
mean_absolute_error(y_true, y_pred)
:第一引数に実績値(test['y']
)、第二引数に予測値(forecast_test['yhat']
) を渡すことで、MAEを計算してくれます。 -
mean_squared_error(...)
:こちらはMSE(平均二乗誤差)を計算します。 -
np.sqrt(...)
:numpy
ライブラリの平方根(square root
)を計算する関数です。MSEの結果にnp.sqrt
を適用することで、RMSEを算出しています。 - f"文字列: {変数:.2f}":これは
f-string
というPython
の便利な文字列フォーマット機能です。{変数:.2f}と書くことで、 変数の値を小数点以下2桁で表示するように整形しています。
最終ステップ':V2モデルによる未来予測
V2モデルの特徴量を使いtrain
+test
の全期間で学習させ予測モデルを完成させます。その後、未来(2025年11月~2026年4月)の未来予測を行います。
期間を分けず全期間を指定し、最良モデルV2で使用した特徴量を抽出
# 1. Prophetで使うためのデータフレームを再作成(新しい特徴量を追加)
df_prophet_final = df[['年月', '最深積雪(cm)','日最高気温の平均(℃)', '降雪量日合計3cm以上日数(日)', '日最高気温0℃未満日数(日)']].copy()
df_prophet_final.rename(columns={'年月': 'ds', '最深積雪(cm)': 'y'}, inplace=True)
指定の通りの期間が対象になったか確認
# 2. データを冬季 (11月, 12月, 1月, 2月, 3月, 4月) のみに絞り込む
df_winter_final = df_prophet_final[df_prophet_final['ds'].dt.month.isin([11, 12, 1, 2, 3, 4])]
# 3. 全データでを使用して予測モデルを構築(train_v2)
train_final = df_winter_final
print(train_final['ds'].min())
print(train_final['ds'].max())
予測モデルを作成
# 4. 新しいProphetモデルを初期化
model_final = Prophet()
最良モデルV2で使用した特徴量「日最低気温0℃未満日数(日)
」でセット
# 5. 3つの特徴量を追加
model_final.add_regressor('日最高気温の平均(℃)')
model_final.add_regressor('降雪量日合計3cm以上日数(日)')
model_final.add_regressor('日最高気温0℃未満日数(日)') # 最良モデルV2で追加した特徴量
予測を実施
model_final.fit(train_final)
最終ステップ'':いよいよ未来予測
予測したい未来のデータフレームを作成します。
期間:2025年11月~2026年4月の冬期シーズン
# 1. 予測したい未来の期間のデータフレームを作成
future_dates_final = pd.date_range(start='2025-11-01', end='2026-04-30', freq='MS')
future_final = pd.DataFrame({'ds': future_dates_final})
期間指定した予測用データフレームに予測モデルと同じにします。
特徴量 | 値 |
---|---|
日最高気温の平均(℃): | mean (平均値) |
降雪量日合計3cm以上日数(日): | mean (平均値) |
日最高気温0℃未満日数(日): | mean (平均値) |
月次集計し平均値(mean
)をセットします。
# 2. 未来の気温・降水量を「過去の各月の平均値」で代用する
# 月(month)ごとに各特徴量の平均値を計算
monthly_avg_final = df_prophet_final.groupby(
df_prophet_final['ds'].dt.month)[['日最高気温の平均(℃)',
'降雪量日合計3cm以上日数(日)', '日最高気温0℃未満日数(日)']].mean()
monthly_avg_final.rename_axis('month', inplace=True)
# future_finalの各行に、月ごとの平均値を割り当てる
future_final['month'] = future_final['ds'].dt.month
future_final = pd.merge(future_final, monthly_avg_final, on='month', how='left')
future_final.drop(columns=['month'], inplace=True)
future_final
準備した未来予測用データフレームで予測を実施
# 3. 準備したfuture_finalデータフレームで予測を実行
regs_final = ['日最高気温の平均(℃)', '降雪量日合計3cm以上日数(日)','日最高気温0℃未満日数(日)']
# 予測を実行
forecast_test_only_final = model_final.predict(future_final[['ds'] + regs_final]).sort_values('ds')
# ☆★☆
# データセット(実績)のふるまいを捉えられているか可視化用
future_all_final = pd.concat([train_final[['ds'] + regs_final], future_final[['ds'] + regs_final]], ignore_index=True)
future_all_final = future_all_final.drop_duplicates('ds').sort_values('ds')
# 学習データも入れた内容で予測を実行
forecast_all_final = model_final.predict(future_all_final).sort_values('ds')
# 予測が変わっていないか念のため検証(最大誤差を表示) ※カンニング対悪
fa_test_part_final = forecast_all_final[forecast_all_final['ds'].isin(future_final['ds'])].sort_values('ds')
max_abs_diff_final = np.max(np.abs(fa_test_part_final['yhat'].to_numpy() - forecast_test_only_final['yhat'].to_numpy()))
# カンニングになっていないかを確認(誤差を確認)
print("max |Δyhat| on test:", max_abs_diff_final)
マイナスの予測値を0に補正
# マイナスの予測値を0に補正
# forecast_test_only
forecast_test_only_final['yhat'] = forecast_test_only_final['yhat'].clip(lower=0)
forecast_test_only_final['yhat_lower'] = forecast_test_only_final['yhat_lower'].clip(lower=0)
forecast_test_only_final['yhat_upper'] = forecast_test_only_final['yhat_upper'].clip(lower=0)
# forecast_all
forecast_all_final['yhat'] = forecast_all_final['yhat'].clip(lower=0)
forecast_all_final['yhat_lower'] = forecast_all_final['yhat_lower'].clip(lower=0)
forecast_all_final['yhat_upper'] = forecast_all_final['yhat_upper'].clip(lower=0)
最終モデルのフィット感を確認
# 予測グラフ描画
start = pd.Timestamp('2015-07-31')
end = pd.Timestamp('2025-07-31')
fig_forecast = model.plot(forecast_all_v3)
plt.xlim(start, end)
plt.show()
実績も予測の範囲内に"まずます"収まっており良いモデルかと思います。
予測部分を抽出して拡大グラフを描画
# 4. 最終的な予測結果をグラフで可視化
fig, ax = plt.subplots(figsize=(10, 6))
# 予測値のプロット
ax.plot(forecast_future_final['ds'], forecast_future_final['yhat'], 'b.-', label='予測値 (最深積雪)')
# 予測の信頼区間
ax.fill_between(forecast_future_final['ds'], forecast_future_final['yhat_lower'], forecast_future_final['yhat_upper'], color='blue', alpha=0.2, label='予測の信頼区間')
# グラフの見た目を整える
ax.set_title('【最終予測】積雪量の未来予測 (2025/11 - 2026/04) - Finalモデル')
ax.set_xlabel('年月')
ax.set_ylabel('最深積雪 (cm)')
ax.legend()
plt.grid(True)
plt.show()
未来予測のテーブルを取得します。
# 5. 最終的な予測結果のテーブルも確認
print("---【最終予測】未来の予測結果 (Finalモデル) ---")
print(forecast_future_final[['ds', 'yhat', 'yhat_lower', 'yhat_upper']])
ステップ+α 数値入りのグラフで描画
数値入りの年別比較データを作成します。
snow_season_df
は再利用でも追いですが、念のためはじめから…
# 11月~4月の期間のデータのみを抽出
snow_season_df = df[df['month'].isin([11, 12, 1, 2, 3, 4])]
# ※※※コピペしてるのでココのシーズン絞り込みは必要ないです※※※
# ★ ここで2014年以降のデータのみに絞り込みます
#winter_final = snow_season_df[snow_season_df['season_year'] > 2014]
# 直近5シーズン(実績のみで抽出)
last5 = (snow_season_df['season_year'].dropna().astype(int).sort_values().unique())[-4:]
winter_final = snow_season_df[snow_season_df['season_year'].isin(last5)]
# 月別に集計(合計。平均/最大にしたい場合は mean()/max() に変更)
monthly_final = (winter_final
.groupby(['season_year','month'], as_index=False)['最深積雪(cm)']
.sum())
monthly_final['label'] = monthly_final['season_year'].astype(str) + '-' + (monthly_final['season_year']+1).astype(str) + '冬'
# ===== 予測(2025-2026シーズン)=====
future_fin = forecast_test_only_final.copy()
future_fin['year'] = future_fin['ds'].dt.year
future_fin['month'] = future_fin['ds'].dt.month
future_fin['season_year'] = np.where(future_fin['month'] >= 11, future_fin['year'], future_fin['year'] - 1)
future_fin = future_fin[future_fin['season_year'] == 2025].copy() # 2025-2026シーズン
future_fin['最深積雪(cm)'] = future_fin['yhat']
future_monthly_fin = (future_fin.groupby(['season_year','month'], as_index=False)['最深積雪(cm)']
.sum())
future_monthly_fin['label'] = '2025-2026冬(予測)'
# ===== 結合 → ピボット =====
all_monthly_fin = pd.concat([monthly_final, future_monthly_fin], ignore_index=True)
# 月の表示順を固定
month_order_fin = [11,12,1,2,3,4]
all_monthly_fin['month'] = pd.Categorical(all_monthly_fin['month'], categories=month_order_fin, ordered=True)
# 行=月, 列=シーズン, 値=最深積雪
pivot_fin = (all_monthly_fin
.pivot_table(index='month', columns='label', values='最深積雪(cm)', aggfunc='sum')
.reindex(index=month_order_fin))
# ===== 描画(カラフルなグループ化棒)=====
plt.figure(figsize=(14,6))
x = np.arange(len(pivot_fin.index))
cols = pivot_fin.columns.tolist()
n = len(cols)
bar_w = 0.8 / max(n,1)
for i, c in enumerate(cols):
y = pivot_fin[c].values
# グループ位置をずらして横並び
xpos = x + (i - n/2 + 0.5)*bar_w
plt.bar(xpos, y, width=bar_w, label=c)
# 値表示(NaNはスキップ)
for xi, yi in zip(xpos, y):
if pd.notna(yi):
plt.text(xi, yi, int(round(yi)), ha='center', va='bottom', fontsize=9)
plt.xticks(x, ['11月','12月','1月','2月','3月','4月'])
plt.title('直近5シーズンと未来予測「2025-2026シーズン」(月別最深積雪)', fontsize=16)
plt.xlabel('月'); plt.ylabel('最深積雪(cm)')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.legend(title='シーズン', ncol=3)
plt.tight_layout()
plt.show()
野沢温泉積雪量予測レポート
【予測概要】
- 予測期間:2025年12月~2026年4月
- 分析手法:Prophet時系列予測モデル + 多重特徴量分析
- データ期間:1984年~2025年(41年間)
- 予測精度:MAE 39.96cm, RMSE 48.37cm
【2025-2026年積雪量予測】
予測時期 | 予測積雪量 | 予測範囲 |
---|---|---|
2025年12月: | 76cm | 23cm-130cm |
2026年01月: | 158cm | 98cm-209cm |
2026年02月: | 183cm | 130cm-237cm |
2026年03月: | 150cm | 96cm-204cm |
2026年04月: | 57cm | 6cm-111cm |
【活用方法】
- スキー場運営計画の参考資料として
- 観光業界の需要予測の基礎データとして
- 除雪計画や設備準備の目安として
- 気象データと組み合わせた総合判断に活用
【重要な注意事項】
予測の限界:
- この予測は過去41年間のデータパターンに基づく統計的推定です。
- 異常気象や急激な気候変動は完全には反映されていません。
- 実際の気象条件により大きく変動する可能性があります。
- 予測区間は統計的な不確実性を示しており、参考値として活用してください。
【主要特徴量の重要度】
特徴量 | 相関 |
---|---|
降雪量合計3cm以上日数(日) | 0.86 |
降雪量合計(cm) | 0.86 |
日最高気温の平均(℃) | -0.79 |
平均気温(℃) | -0.78 |
【分析手法】
-
Prophet時系列予測モデル
Facebookが開発した時系列予測ライブラリを使用。季節性とトレンドを自動学習し、外部変数(気温、降雪量等)を組み込んだ多変量予測を実現。- モデル設定:標準(
Default
) - 使用特徴量の追加:
add_regresor
(使用する特徴量)で追加
- モデル設定:標準(
【データ分割戦略】
- 訓練データ:1984年1月~2023年12月(40年間)
- テストデータ:2024年1月~2025年3月(検証用)
- 最良モデル:1984年1月~2025年7月(最終予測モデル)
- 予測期間:2025年12月~2026年4月(最終目標)
【特徴量】
- 日最高気温の平均(℃)
- 降雪量日合計3cm以上日数(日)
- 日最高気温0℃未満日数(日)
【検証方法】
- MAE(平均絶対誤差):予測値と実測値の差の絶対値の平均
- RMSE(二乗平均平方根誤差):大きな誤差により敏感な指標
- 散布図分析:実測値vs予測値の相関確認
- 時系列プロット:予測トレンドの妥当性確認
まとめ
今回のプロジェクトを通して、データ分析の一連のプロセスを体験することができました。特に、仮説を立てて特徴量を追加し、評価指標で改善を確認しながらの進め方は、非常に実践的で学びが多かったです。
最後までお読みいただき、ありがとうございました!