桜の満開を予測する②（LSTMを試す）

import numpy as np
import ast
from datetime import datetime

def extract_features(weather_data):

    daily_data = ast.literal_eval(weather_data)  # 文字列をリストに変換
    print(datetime.utcfromtimestamp(daily_data[0]["dt"]))

    return {
        'date': datetime.utcfromtimestamp(daily_data[0]["dt"]),
        'temp': daily_data[0]["temp"],
        'feels_like': daily_data[0]["feels_like"],
        'pressure': daily_data[0]["pressure"],
        'humidity': daily_data[0]["humidity"],
        'dew_point': daily_data[0]["dew_point"],
        'dew_point': daily_data[0]["dew_point"],
        'clouds': daily_data[0]["clouds"],
        'wind_speed': daily_data[0]["wind_speed"],
        'wind_deg': daily_data[0]["wind_deg"],
    }

features_df = combined_df['data'].apply(extract_features)
features_df = pd.DataFrame(features_df.tolist())

import pandas as pd
import numpy as np
from datetime import datetime
import torch
import torch.nn as nn
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler

# 前処理
# 年・月・日を追加し、1月〜4月に絞る
features_df['date'] = pd.to_datetime(features_df['date'])
features_df['year'] = features_df['date'].dt.year
features_df['month'] = features_df['date'].dt.month
features_df['day'] = features_df['date'].dt.day
features_df = features_df[features_df['month'].between(1, 4)]

# 特徴量を標準化
feature_cols = [col for col in features_df.columns if col not in ['date', 'year', 'month', 'day']]
scaler = StandardScaler()
features_df[feature_cols] = scaler.fit_transform(features_df[feature_cols])

# 満開日（ターゲット）と年情報の取得
tokyo_df['満開の日'] = pd.to_datetime(tokyo_df['満開の日'])
tokyo_df['year'] = tokyo_df['満開の日'].dt.year
tokyo_df['day_of_year'] = tokyo_df['満開の日'].dt.dayofyear
y_map = tokyo_df.set_index('year')['day_of_year'].to_dict()

# 年ごとの時系列データを生成
SEQUENCE_LENGTH = 110

X_seq, y_seq, years_seq = [], [], []

for year, group in features_df.groupby('year'):
    group_sorted = group.sort_values('date')
    if len(group_sorted) >= SEQUENCE_LENGTH and year in y_map:
        values = group_sorted[feature_cols].values
        X_seq.append(values[:SEQUENCE_LENGTH])  # ← 最初の90日間だけを使う
        y_seq.append(y_map[year])  # 満開日
        years_seq.append(year)
        
X_seq = np.array(X_seq)
y_seq = np.array(y_seq)

class BlossomDataset(torch.utils.data.Dataset):
    def __init__(self, X, y, years):
        self.X = torch.tensor(X, dtype=torch.float32)
        self.y = torch.tensor(y, dtype=torch.float32).unsqueeze(1)
        self.years = torch.tensor(years, dtype=torch.int)

    def __len__(self):
        return len(self.X)

    def __getitem__(self, idx):
        return self.X[idx], self.y[idx], self.years[idx]

dataset = BlossomDataset(X_seq, y_seq, years_seq)

# 年で train/test 分割
train_idx = [i for i, y in enumerate(years_seq) if y < 2020]
test_idx = [i for i, y in enumerate(years_seq) if y >= 2020]

train_dataset = torch.utils.data.Subset(dataset, train_idx)
test_dataset = torch.utils.data.Subset(dataset, test_idx)

train_loader = torch.utils.data.DataLoader(train_dataset, batch_size=4, shuffle=True)
test_loader = torch.utils.data.DataLoader(test_dataset, batch_size=4, shuffle=False)

# LSTMモデル
class LSTMModel(nn.Module):
    def __init__(self, input_size, hidden_size=64):
        super().__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        _, (hn, _) = self.lstm(x)
        return self.fc(hn[-1])

model = LSTMModel(input_size=len(feature_cols))
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)

early_stop_threshold = 0.05
for epoch in range(1000):
    model.train()
    for X_batch, y_batch, _ in train_loader:  # 年は無視
        optimizer.zero_grad()
        output = model(X_batch)
        loss = criterion(output, y_batch)
        loss.backward()
        optimizer.step()
    if epoch % 50 == 0:
        print(f"Epoch {epoch}, Loss: {loss.item():.4f}")
    # ある損失以下なら中断
    if loss.item() < early_stop_threshold:
        print(f"Early stopping at epoch {epoch}, loss: {loss.item():.4f}")
        break

# 評価
model.eval()

y_true, y_pred, y_years = [], [], []

with torch.no_grad():
    for X_batch, y_batch, year_batch in test_loader:
        pred = model(X_batch)
        y_true.extend(y_batch.numpy().flatten())
        y_pred.extend(pred.numpy().flatten())
        y_years.extend(year_batch.numpy().flatten())

# 結果表示
for year, doy in zip(y_years, y_pred):
    predicted_date = datetime(year, 1, 1) + pd.to_timedelta(doy - 1, unit='D')
    print(f"{year}年の予測満開日: {predicted_date.strftime('%Y-%m-%d')}")

Epoch 0, Loss: 8122.0400
Epoch 50, Loss: 3703.0417
Epoch 100, Loss: 1634.8743
Epoch 150, Loss: 432.9744
Epoch 200, Loss: 30.7069
Epoch 250, Loss: 8.4507
Epoch 300, Loss: 106.8499
Epoch 350, Loss: 26.9567
Epoch 400, Loss: 6.4260
Epoch 450, Loss: 108.6270
Epoch 500, Loss: 47.5118
Epoch 550, Loss: 59.1875
Early stopping at epoch 587, loss: 0.0033

2020年の予測満開日: 2020-04-01
2021年の予測満開日: 2021-04-02
2022年の予測満開日: 2022-04-02
2023年の予測満開日: 2023-04-02
2024年の予測満開日: 2024-04-01

2020-03-22
2021-03-22
2022-03-27
2023-03-22
2024-04-04