日本の大気中の二酸化炭素濃度予測
はじめに
こんにちは、今回アイデミーで学んだ事をもとに日本の大気中の二酸化炭素濃度の予測モデルを作成し、こちらにまとめます。
私は社会を取り巻く環境が変化する中で自身でも知識を身に着けたいと考え、アイデミーのデータ分析講座を受講しました。学習を始める前はプログラミングの基礎知識等もなく、全くの素人です。
本記事の概要
プログラミング初心者だけど
将来の地球環境はどのようになるの？大気中の二酸化炭素濃度を予測をしてみたいなんて思っている方もいると思います。
今回は、季節ARIMA（SARIMA:Seasonal Autoregressive Integrated Moving Average）モデルを使って予測をしてみました。
【読んで欲しい方】
・機械学習に興味がある人
・データサイエンスに興味がある人
実行環境
・Google Colaboratory
・Google Drive
データ準備
①データ調査
〈大気中の二酸化炭素濃度〉
ダウンロードサイト
気象庁 大気・海洋環境観測年報 データダウンロードhttps://www.data.jma.go.jp/gmd/env/data/report/data/download/atm_bg_j.html
大気中の温室効果ガス 観測種目 大気二酸化炭素
<< 大気バックグランド汚染観測地点一覧 >>
*********************************************************************************
地点 WMO地点番号 緯度 経度 高度
綾里 47513 39゜02'N 141゜49'E 260 m
南鳥島 47991 24゜17'N 153゜59'E 7 m
与那国島 47912 24゜28'N 123゜01'E 30 m
*********************************************************************************
・使用データ 月別値
・期間 1987年1月～2022年2月
②データの前処理
３地点のダウンロードデータを統合し、csvデータに変換し、Googleドライブへ格納
〈データの読み込み〉
pandasのpd.read_csv()メソッドで文字コードも指定して読み込み
データの可視化
from google.colab import drive
drive.mount('/content/drive')import pandas as pd
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import statsmodels.api as sm
import matplotlib.cm as cm
import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_errordata_df = pd.read_csv('/content/drive/MyDrive/aidemy/CO2monthly.csv', encoding='shift-jis')
#データの形を表示 →420行、5列のデータがある
print(f'Data_df_shape : {data_df.shape}\n')
# 各columnのデータ型を表示→整数(int), 小数(float), 文字(object)があります。
print(f'{data_df.dtypes} \n')
#先頭5つを可視化
display(data_df.head())
予測モデルとして使用できよう、前処理を行う。
①欠損値の代入
②外れ値の処理
〈①欠損値の代入〉
読み込んだデータに、"NaN"が連続する箇所がある。
これは、三地点の観測開始年月が異なる為、データが存在していない等と考えられる。それ以外の欠損値には、三地点の地理的関係から判断し、データが存在している他地点の平均値を代入する。
#欠損値はあるかどうか 、いくつあるか→綾里はないが、南鳥島に73、与那国に120ある
data_df.isnull().sum()
観測地点
綾里、南鳥島、与那国のデータから、欠損値の数によって代入するデータを決定し、入力する。
# 右から3列を取得
last_three_columns = data_df.iloc[:, -3:]
# 各行に対して処理を行う
for index, row in last_three_columns.iterrows():
# 3つのデータのうち1つがNaNの場合
if row.count() == 2:
# NaNでない数値を取得して平均値とする
mean_value = row[~row.isna()].mean()
last_three_columns.iloc[index] = row.fillna(mean_value)
# 3つのデータのうち2つがNaNの場合
elif row.count() == 1:
# NaNでない数値をそのまま代入する
last_three_columns.iloc[index] = row.fillna(row.dropna().iloc[0])
# 元のDataFrameに変更を反映する
data_df.iloc[:, -3:] = last_three_columns
# 変更を保存
data_df.to_csv('data_processed.csv', index=False)
print(data_df)
その後、観測年・観測月のデータ列を統合し削除。
#観測年と観測月を結合して削除
data_df['Date'] = data_df['観測年'].astype(int).astype(str) + '-' + data_df['観測月'].astype(int).astype(str).str.zfill(2)
data_df.insert(0, 'Date', data_df.pop('Date'))
# 結果を表示
print(data_df)
〈②外れ値の処理〉
それぞれの地点で外れ値が以下の通りある。
# データ列を選択
data_column = data_df["綾里"]
# Zスコアを計算
z_scores = stats.zscore(data_column)
# Zスコアが特定の閾値を超えるデータポイントを特定
threshold = 3 # 通常は2〜3の範囲が使用されます
outliers = (z_scores > threshold) | (z_scores < -threshold)
# 外れ値を表示
outlier_values = data_column[outliers]
print("外れ値:", outlier_values)
# 南鳥島
data_column = data_df["南鳥島"]
z_scores = stats.zscore(data_column)
threshold = 3
outliers = (z_scores > threshold) | (z_scores < -threshold)
outlier_values = data_column[outliers]
print("外れ値:", outlier_values)
# 与那国
data_column = data_df["与那国"]
z_scores = stats.zscore(data_column)
threshold = 3
outliers = (z_scores > threshold) | (z_scores < -threshold)
outlier_values = data_column[outliers]
print("外れ値:", outlier_values)
データの連続性を鑑み、それぞれの外れ値は前後の観測年月の平均値を入れることとする。
for index, row in data_df.iterrows():
if row['綾里'] == 99999.99:
# 前のデータと後のデータを取得
prev_row = data_df.loc[index - 1]
next_row = data_df.loc[index + 1]
# 平均値を計算
avg_value = (prev_row['綾里'] + next_row['綾里']) / 2
# 外れ値を平均値で置き換え
data_df.at[index, '綾里'] = avg_value
#南鳥島 外れ値が連続するため、以下の方法で行う #forループはデータフレーム内の各行を反復処理します 。
for index, row in data_df.iterrows():
#各行の綾里の値が 「999」である場合、前後のデータを確認するためのインデックスを設定します。
if row['南鳥島'] == 99999.99:
prev_index = index - 1
next_index = index + 1
#whileループを使用して 、前のデータと後のデータについて「999」でないデータを見つけるための処理を行います。
while prev_index >= 0 and data_df.at[prev_index, '南鳥島'] == 99999.99:
prev_index -= 1
while next_index < len(data_df) and data_df.at[next_index, '南鳥島'] == 99999.99:
next_index += 1 #前後のデータが見つかった場合 、それらのデータを使用して平均値を計算し、外れ値を平均値で置き換えます。
if prev_index >= 0 and next_index < len(data_df):
avg_value = (data_df.at[prev_index, '南鳥島'] + data_df.at[next_index, '南鳥島']) / 2
data_df.at[index, '南鳥島'] = avg_value
#与那国は綾里と同様
for index, row in data_df.iterrows():
if row['与那国'] == 99999.99:
# 前のデータと後のデータを取得
prev_row = data_df.loc[index - 1]
next_row = data_df.loc[index + 1]
# 平均値を計算
avg_value = (prev_row['与那国'] + next_row['与那国']) / 2
# 外れ値を平均値で置き換え
data_df.at[index, '与那国'] = avg_value
data_df.to_csv('cleaned_data.csv', index=False) # 置き換え後のデータを保存データの前処理は以上です。
以降は、モデルの作成と予測を行います。
データの可視化
データをグラフ化して、トレンド、季節性、周期性、ノイズなどの特性を観察します。これにより、どの種類のモデルを使用するかを決定するのに役立ちます。
plt.figure(figsize=(20, 8)) # グラフのサイズを調整
# 各列のデータをプロット
plt.plot(data_df["Date"], data_df["綾里"], label="Ryori")
plt.plot(data_df["Date"], data_df["南鳥島"], label="Minamitorishima")
plt.plot(data_df["Date"], data_df["与那国"], label="Yonaguni")
# X軸のラベルを１２か月ごとに表示
x_labels = data_df["Date"][::12] # １２か月ごとのラベルを取得
x_ticks = np.arange(len(data_df["Date"]))[::12] # １２か月ごとの位置を取得
#X軸のラベルを45度傾ける
plt.xticks(x_ticks, x_labels, rotation=45)
# グラフのタイトルと軸ラベルを設定
plt.title("Changes in Atmospheric Carbon Dioxide Concentration JAPAN")
plt.xlabel("DATE")
plt.ylabel("ppm")
# 凡例を表示
plt.legend()
# グリッドを表示
plt.grid(True)
# グラフを表示
plt.show()
今回使用するデータは数値が徐々に増加する傾向がみられる。
このデータを、時系列に並んだ1つのデータを長期的な傾向（トレンド）、周期的な変動（季節変動）、不規則な変動（残渣）に分解したうえで、それらのデータを統計的な手法を用いて将来を予測する。
今回使用するデータを処理前のデータ（原系列）、トレンド、季節変動、残渣に分けてみます。
# 日付列を作成
data_df['Date'] = data_df['観測年'].astype(str) + '-' + data_df['観測月'].astype(str).str.zfill(2)
# 日付列を日付型に変換し、インデックスに設定
data_df['Date'] = pd.to_datetime(data_df['Date'])
data_df.set_index('Date', inplace=True)
# 綾里の季節調整を行う
decomposition = sm.tsa.seasonal_decompose(data_df['綾里'], model='additive')
# トレンド、季節変動、残差を取得
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
# トレンド、季節変動、残差のグラフを描画
plt.figure(figsize=(8, 6))
plt.subplot(411)
plt.plot(data_df['綾里'], label='原系列')
plt.legend(loc='upper left')
plt.subplot(412)
plt.plot(trend, label='トレンド')
plt.legend(loc='upper left')
plt.subplot(413)
plt.plot(seasonal, label='季節変動')
plt.legend(loc='upper left')
plt.subplot(414)
plt.plot(residual, label='残差')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
南鳥島、与那国も同様にプロット
４種類に分解してみた結果、トレンドは”正のトレンド”、季節変動があり、一定周期であることがわかりました。
時系列データに特有な統計量である「自己相関係数」について見ていく。
同じ時系列データでの別々の時点同士での共分散で、これを様々な値間で比べられるようにしたものを「自己相関係数」という。簡単にいうと自己相関係数とは今のデータx年が過去のデータy年と比べてどれだけ似ているか（影響を受けているか）を数値で表した値となる。
自己相関を可視化した。
# 自己相関関数（ACF）を計算
acf_result = sm.tsa.acf(data_df['綾里'], nlags=40) # nlagsはラグの数を指定
# ACFのグラフを描画
plt.figure(figsize=(6, 3))
plt.stem(range(len(acf_result)), acf_result)
plt.title('ACF(ryouri)')
plt.xlabel('ラグ')
plt.ylabel('自己相関')
plt.show()
その結果、自己相関関係の図が波を描きながら右下がりになっており、時間的なパターンやトレンドが存在するデータセットに関して考えられる。自己相関が一定の周期的な変動を示し、その周期が時間とともに減少していることを示唆している。時間的な依存関係が存在し、その依存関係が時間が経つにつれて減少していることを示す。つまり、過去のデータが将来のデータに対して強い影響を持っていたが、時間が経つにつれてその影響が弱まっている可能性がある。
南鳥島、与那国でも同様の結果となった。
# 自己相関関数（ACF）を計算
acf_result = sm.tsa.acf(data_df['南鳥島'], nlags=40) # nlagsはラグの数を指定
# ACFのグラフを描画
plt.figure(figsize=(6, 3))
plt.stem(range(len(acf_result)), acf_result)
plt.title('ACF(minamitorishima)')
plt.xlabel('ラグ')
plt.ylabel('自己相関')
plt.show()
# 自己相関関数（ACF）を計算
acf_result = sm.tsa.acf(data_df['与那国'], nlags=40) # nlagsはラグの数を指定
# ACFのグラフを描画
plt.figure(figsize=(6, 3))
plt.stem(range(len(acf_result)), acf_result)
plt.title('ACF(yonaguni)')
plt.xlabel('ラグ')
plt.ylabel('自己相関')
plt.show()
モデルの作成と予測
今回は、前述の結果から、季節ARIMA（SARIMA:Seasonal Autoregressive Integrated Moving Average）モデルによって予測を行う。
〈SARIMAモデルとは〉
SARIMAモデルを説明する前に、ARIMAモデルについて整理する。
ARIMAモデルとは、「autoregressive integrated moving average」の略で、自己回帰モデル(ARモデル)、移動平均モデル(MAモデル)、和分モデル(Iモデル)の3モデルを組み合わせたモデルである。
SARIMAモデルとは、「Seasonal AutoRegressive Integrated Moving Average」の略で、ARIMAモデルに「季節的な周期パターン」を加えたモデルである。
つまり、季節変動があるデータに対してARIMAモデルを拡張した手法が、SARIMAモデルであり、具体的には、時系列方向にARIMAモデルを使い、かつ、周期方向にもARIMAモデルを使っているモデルである。
# 列名を変更
data_df.rename(columns={'綾里': 'ryori', '南鳥島': 'minamitorishima', '与那国': 'yonaguni'}, inplace=True)
# データをトレーニングデータとテストデータに分割
train_data = data_df[data_df['観測年'] < 2021] # 2021年より前をトレーニングデータとする
test_data = data_df[data_df['観測年'] == 2021] # 2021年のデータをテストデータとする
# SARIMAモデルをトレーニングして予測する関数を定義
def train_sarima_predict(train, test, location):
# SARIMAモデルを設定
order = (1, 1, 1) # (p, d, q)オーダーを設定
seasonal_order = (1, 1, 1, 12) # 季節性オーダーを設定
# SARIMAモデルをトレーニング
model = SARIMAX(train[location], order=order, seasonal_order=seasonal_order)
results = model.fit()
# テストデータに対する予測を行う
predictions = results.predict(start=len(train), end=len(train) + len(test) - 1, dynamic=False)
return predictions
# 一つのグラフに綾里、南鳥島、与那国の予測結果とテストデータをまとめて表示
plt.figure(figsize=(12, 6))
for location in ['ryori', 'minamitorishima', 'yonaguni']:
predictions = train_sarima_predict(train_data, test_data, location)
plt.plot(test_data.index, predictions, label=f'{location} SARIMA Predictions')
# テストデータをプロット
for location in ['ryori', 'minamitorishima', 'yonaguni']:
plt.plot(test_data.index, test_data[location], label=f'Test Data - {location}')
plt.legend()
plt.title('SARIMA Predictions and Test Data for Multiple Locations')
plt.show()
モデルの評価
予測の精度を確認。平均絶対誤差（MAE）、平均二乗誤差（MSE）により評価をしていく。
〈平均絶対誤差（Mean Absolute Error, MAE）〉
予測モデルの精度を評価するための一般的な指標の一つ。予測値と実際の値との差の絶対値を取り、それらの絶対値の平均を計算。
（参考）
精度が低い
MAEが0.5以上:
一般的には予測モデルの精度が低いと見なす。予測値と実際の値の誤差が大きく、改善が必要な場合がある。
精度が高い
MAEが0.1未満:
非常に高い予測精度を持つモデルとされる。予測値と実際の値の誤差が非常に小さいことを示す。
MAEが0.1から0.2:
高い予測精度を持つモデルとされる。誤差が小さく、実用的に非常に有用なモデルとされることが多い。
MAEが0.2から0.3:
予測精度はまずまずとされるが、改善の余地がある場合もある。
MAEが0.3から0.5:
一般的には中程度の予測精度とされ、モデルは利用可能だが、改善の検討が必要な場合もある。
〈平均二乗誤差（Mean Squared Error, MSE）〉
予測モデルの精度を評価するための一般的な指標の一つ。予測値と実際の値との差の二乗を取り、それらの二乗誤差の平均を計算。
（参考）
小さいMSEの基準（予測モデルの精度が高いと見なされる）
MSEが0未満: 実際には発生しない場合であり、理論的に不可能
MSEが0に近い: MSEが非常に小さい場合、予測モデルは非常に精度が高いと見なすことができる。予測値と実際の値がほぼ一致していることを示す。MSEが0.1未満: 非常に高い予測精度を持つモデル。誤差が非常に小さく、実用的に非常に有用なモデルとされることが多い。
大きいMSEの基準（予測モデルの精度が低いと見なされる）:
MSEが1以上: 一般的には予測モデルの精度が低いと見なす。予測値と実際の値の誤差が大きく、改善が必要。
MSEが10以上: 予測モデルの精度が非常に低いと見なされ、ほとんどの場合、モデルの改善が急務とされる。
#各場所 （綾里、南鳥島、与那国）に対して、予測値と実際の値の誤差を計算
for location in ['ryori', 'minamitorishima', 'yonaguni']:
predictions = train_sarima_predict(train_data, test_data, location)
actual = test_data[location]
mae = mean_absolute_error(actual, predictions)
mse = mean_squared_error(actual, predictions)
print(f'{location} - MAE: {mae}, MSE: {mse}')
#計算した誤差を表にまとめるために 、PandasのDataFrameを使用
error_data = pd.DataFrame(columns=['Location', 'MAE', 'MSE'])
for location in ['ryori', 'minamitorishima', 'yonaguni']:
predictions = train_sarima_predict(train_data, test_data, location)
actual = test_data[location]
mae = mean_absolute_error(actual, predictions)
mse = mean_squared_error(actual, predictions)
error_data = error_data.append({'Location': location, 'MAE': mae, 'MSE': mse}, ignore_index=True)
print(error_data)Location MAE MSE
0 ryori 0.950809 1.254300
1 minamitorishima 0.388628 0.241253
2 yonaguni 0.718216 0.636499
各場所（ryori、minamitorishima、yonaguni）のモデルを評価する。
ryori:
MAE: 0.950809
MSE: 1.254300
この場所のモデルは、MAEが0.95と比較的高い値を示しており、予測精度は中程度であると言える。また、MSEが1.25とやや高い値を示しており、モデルが改善の余地があるといえる。
minamitorishima:
MAE: 0.388628
MSE: 0.241253
この場所のモデルは、MAEが0.39と比較的低い値を示しており、予測精度は高いと言える。また、MSEが0.24と非常に小さい値を示していおり、モデルの予測精度は非常に高いと言える。
yonaguni:
MAE: 0.718216
MSE: 0.636499
この場所のモデルは、MAEが0.72と中程度の値を示しており、予測精度は一般的に受け入れられる範囲である。また、MSEが0.64と中程度の値を示しており、モデルの予測精度は良好と言えますが、改善の余地もあるかもしれない。
総括すると、minamitorishimaのモデルは他の場所に比べて非常に高い予測精度を持っており、ryoriのモデルは一般的な予測精度を持っている。yonaguniのモデルも予測精度は中程度であり、改善の余地はあるかもしれないが、一般的に受け入れられる範囲であると考える。
考察
・３つの地点で同じモデルを使用して予測を行えば、モデル評価は同様の結果が得られると考えいたが、実際は精度に差が出るものであるということを理解できた。測定値の地理的な環境や気温などが関係しているのかも知れない。今回の予測を行ったことで、予測のモデルについての理解を深めることができた。
・大気中の二酸化炭素濃度が、右肩あたりで毎年増加しているという事実を再認識することができた。
この記事が気に入ったらサポートをしてみませんか？