AIプログラミング初心者が日本の漁獲高を予測するモデルを作ってみた

漁業家系で育ったエンジニアがAIを使って、未来の日本の漁獲高を予測するモデルを作ってみた。

#目次
1.実装環境
2.利用データ
3.予測対象について
4.データの可視化
5.日本の漁獲高予測モデルを作成
6.まとめ
7.感想

#実装環境

• MacBook Pro
• Python3
• Visual Studio Code

#利用データ

• 日本の漁獲高
• 日本の漁業人口
• 隣国の人口
• 世界気温の変動
• 日本海の海面温度

#予測対象について

今回は、農林水産省の統計情報から取得した1999年〜2020年までの年度毎のデータ（図１）をから、各魚種の未来の漁獲高を予想するモデルを作成することにした。
※目的変数を各魚種の漁獲高とした
図１：日本の漁獲高(1999年〜2020年・44魚種）

データの可視化

まずは、全データの傾向を把握するために、データを可視化してみる。

日本の漁獲高

４4魚種をグラフ化して漁獲高の推移を並べると次の様になった。（X軸のラベルは、1999年〜2020年）

最終的には、これらのグラフと近似する予測値を出せるモデルを構築する。

日本の漁業人口

続いて、日本の漁業人口をグラフ化。

漁獲高と等しく減少傾向。目的変数（漁獲高）との相関が強そう。
画像からは読みとり難いが、データが不足している。2020年までのデータが欲しかったが取得できなかった。

時系列解析（SARIMAモデル）を利用し不足データを補強

時系列解析を行い、不足データを補強した。（以下画像はその結果）

プログラムは以下。

import warnings
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import os
import statsmodels.api as sm
import itertools

root_path = os.getcwd() +'/aidemy_blog'
data_path = root_path +'/maindata'

# データの読み込み
fisherman = pd.read_csv(data_path+'/fisherman.csv')
index = pd.date_range('1961-12-31','2017-12-31',freq='Y')
fisherman.index = index
del fisherman['year']
fisherman.index.name = 'year'

# 関数の定義
def selectparameter(DATA,s):
    p = d = q = range(0, 2)
    pdq = list(itertools.product(p, d, q))
    seasonal_pdq = [(x[0], x[1], x[2], s) for x in list(itertools.product(p, d, q))]
    parameters = []
    BICs = np.array([])
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(DATA,
                                            order=param,
                                            seasonal_order=param_seasonal)
                results = mod.fit()
                parameters.append([param, param_seasonal, results.bic])
                BICs = np.append(BICs,results.bic)
            except:
                continue
    return parameters[np.argmin(BICs)]

# 周期を埋めて最適なパラメーターを求めてください
params = selectparameter(fisherman,10)
print(params[0])
print(params[1])
# モデルの当てはめ
SARIMA_fisherman = sm.tsa.statespace.SARIMAX(fisherman,order=params[0],seasonal_order=params[1]).fit()

#構築したSARIMAモデルのBICを出力します
print(SARIMA_fisherman.bic)

#predに予測データを代入する
pred = SARIMA_fisherman.predict('2015-12-31', '2040-12-31')

#predデータの可視化
plt.plot(fisherman)
plt.savefig(data_path+'/nomal_fisherman.png')
plt.plot(pred,color='r')
plt.savefig(data_path+'/sarima_fisherman.png')

fisherman = pd.concat([fisherman,pred],axis=0)
fisherman.to_csv(data_path+'/sarima_fisherman.csv')

不足していたデータの補強に成功。

隣国の人口

隣国（中国、韓国、北朝鮮）の人口の推移をグラフ化。

突き抜けているのが中国。中国以外のデータが予測に影響を与えることは少なそう。

世界気温の変動

世界の基本の変動をグラフ化。

じわじわ平均気温が上がっているのが読み取れる。

日本海の海面温度

気象庁からデータを取得したかった。が、、、入手できたデータは次の画像データのみ。

このような画像データを、スクレイピングで1982年２月〜２０２２年1月分まで取得。（合計468枚）

さて、困った。。。これを数値化しないと分析に使えない。

教師あり学習（分類）で画像データを数値化

上記の海水温のサーモグラフィ画像を100エリアに分断し、各エリアの平均温度を求める事にした。幸いなことに判例が用意されているため、RGBの値から海水温の予測が可能である。

エリア分けのイメージ

まずは、次のコードで必要な部分（サーモグラフィ）のみにトリミング。

import numpy as np
import os
from PIL import Image

'''
気象庁から取得した海水温データをトリミングし、
450:450の画像データを生成する処理。
'''
def sea_img_trm(imput_path,save_path):
  img = np.array(Image.open(imput_path))
  trmimg = img[37:487,57:507]
  Image.fromarray(trmimg).save(save_path)

# 画像のトリミング
img_path = os.getcwd()+'/aidemy_blog/sea_img'
files = os.listdir(img_path)
for filename in files:
  save_path = img_path + "_trm/" + filename
  sea_img_trm(img_path+'/'+filename,save_path)

続いてトリミングした画像を45×45サイズに分断し100エリアに分ける。

import numpy as np
import os
from PIL import Image

'''
海水温画像データ(450:450)を(45:45)のエリアで分断し、
100個の画像に分ける処理。
'''
def mk(drpath):
  if not os.path.exists(drpath):
    # ディレクトリが存在しない場合、ディレクトリを作成する
    os.makedirs(drpath)

def sea_img_more_trm(imput_path,filename):
  img = np.array(Image.open(imput_path))
  count =0
  for x in range(1,11):
    x1 = (x-1)*45
    x2 = (x)*45
    for y in range(1,11):
      count +=1
      y1 =(y-1)*45
      y2 =(y)*45
      trmimg = img[x1:x2,y1:y2]
      if count <10:
        mk(os.getcwd()+'/aidemy_blog/sea_img_area/00'+str(count))
        save_path = os.getcwd()+'/aidemy_blog/sea_img_area/00'+str(count) + '/'+filename
      elif count < 100:
        mk(os.getcwd()+'/aidemy_blog/sea_img_area/0'+str(count))
        save_path = os.getcwd()+'/aidemy_blog/sea_img_area/0'+str(count) + '/'+filename
      else:
        mk(os.getcwd()+'/aidemy_blog/sea_img_area/'+str(count))
        save_path = os.getcwd()+'/aidemy_blog/sea_img_area/'+str(count) +'/'+filename
      Image.fromarray(trmimg).save(save_path)

# 画像のトリミング
img_path = os.getcwd()+'/aidemy_blog/sea_img_trm'
files = os.listdir(img_path)
for filename in files:
  sea_img_more_trm(img_path+'/'+filename,filename)

続いて、RGBの値から海水温を推測するために、凡例を数値化する。

import imp
import numpy as np
import os
from PIL import Image
import pandas as pd


work_path = os.getcwd() + '/aidemy_blog'
data_path = work_path + '/data'
model_path = work_path +'/model'
'''
海水温データの判例を数値化する処理。
温度(-2〜31):RGB のcsvデータを作成する
'''
img = np.array(Image.open(os.getcwd()+'/aidemy_blog/sea_img/0.png'))
print(img[0])

seawater_temperature=(33/img.shape[0])
j =0
rgb_temp_list=[]
for i in img:
  rgb_temp_list.append([round(31-(seawater_temperature*j),3),i[0,0],i[0,1],i[0,2]])
  j += 1

df = pd.DataFrame(np.array(rgb_temp_list),columns=['temperature','R','G','B'])
df.to_csv(data_path+'/rgb_temperature.csv')

これによって作成された「rgb_temperature.csv」は以下。

海面100エリアの画像と、温度とRGBの対応表ができた。
ここから、教師あり学習（分類）を使って海面100エリアの年毎の平均海温度（数値）を予測する。

k-nn法と、ランダムフォレスト、決定木でモデリングし、スコアを比較。

import pandas as pd
import os
import sklearn
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier

work_path = os.getcwd() + '/aidemy_blog'
data_path = work_path + '/data'
model_path = work_path +'/model'

#海水温とRGBの対応表を読み込む
img = pd.read_csv(data_path+'/rgb_temperature.csv')
#説明変数：RGB
X = img[['R','G','B']].astype(int)
#目的変数：温度
y = img['temperature'].astype(int)

# 訓練データと検証データに分ける
train_X, test_X, train_y, test_y = sklearn.model_selection.train_test_split(X,y,random_state=42)

print()
print('K-NN法')
# モデルの構築
model = KNeighborsClassifier()
# モデルの学習
model.fit(train_X, train_y)
# 正解率の算出
print(model.score(test_X, test_y))
print()
print('ランダムフォレスト')
# モデルの構築
model = RandomForestClassifier()
# モデルの学習
model.fit(train_X, train_y)
# 正解率の算出
print(model.score(test_X, test_y))
print()
print('決定木')
# モデルの構築
model = DecisionTreeClassifier()
# モデルの学習
model.fit(train_X, train_y)
# 正解率の算出
print(model.score(test_X, test_y))

結果はこちら。

どのモデルも遜色ない。
それぞれ予測データを作成して中身を比較し使うデータを決定する。
以下、K-NN法、ランダムフォレスト、決定木を用いた海面100エリアの海面温度予測。

import imp
import numpy as np
from PIL import Image
import statistics
import pandas as pd
import os
import re
from glob import glob
import pickle
import warnings

warnings.simplefilter('ignore')

work_path = os.getcwd() + '/aidemy_blog'
data_path = work_path + '/data'
model_path = work_path +'/model_sea'
sea_area_path = work_path + '/sea_img_area/'

labels =[]
indexs =[]
for dr in glob(sea_area_path+'/*/'):
  labels.append(re.sub('.*/','',dr[:-1]))
for f in  os.listdir(sea_area_path+'/001'):
  indexs.append(f.replace('.png',''))

labels = ['ym']+sorted(labels)
indexs = sorted(indexs)

clf_decisiontree = object
clf_knn = object
clf_randomforest = object
data_decisiontree =[]
data_knn = []
data_randomforest = []

#保存していた学習済みモデルを呼び出す
with open(model_path+'/model_sea_decisiontreee.pickle', mode='rb') as f:
  clf_decisiontree = pickle.load(f)
with open(model_path+'/model_sea_knn.pickle', mode='rb') as f:
  clf_knn = pickle.load(f)
with open(model_path+'/model_sea_randomforest.pickle', mode='rb') as f:
  clf_randomforest = pickle.load(f)

#海面100エリアのRGBデータをまとめる
for index in indexs:
  row_decisiontree=[index]
  row_knn=[index]
  row_randomforest=[index]

  for label in labels:
    if label != 'ym':
      img = np.array(Image.open(sea_area_path +label +'/'+ index+ '.png'))
      r =[]
      g =[]
      b =[]

      #45×45の画像のr,g,bの平均値を取得する
      for rgblist in img:
        for rgb in rgblist:
          r.append(rgb[0])
          g.append(rgb[1])
          b.append(rgb[2])
      rgbmean = [[statistics.mean(r),statistics.mean(g),statistics.mean(b)]]

      #RGBの平均値が等しい＝画像のグレーの部分は「99」にする
      if rgbmean[0][0] == rgbmean[0][1] and rgbmean[0][0] == rgbmean[0][2]:
        row_decisiontree.append(99)
        row_knn.append(99)
        row_randomforest.append(99)
      
      #グレーの部分以外は各モデルで温度を予測する
      else:
        row_decisiontree.append(clf_decisiontree.predict(rgbmean)[0])
        row_knn.append(clf_knn.predict(rgbmean)[0])
        row_randomforest.append(clf_randomforest.predict(rgbmean)[0])

  data_decisiontree.append(row_decisiontree)
  data_knn.append(row_knn)
  data_randomforest.append(row_randomforest)

#予測結果をcsvとして出力
data_decisiontree = pd.DataFrame(data_decisiontree,columns=labels)
data_decisiontree.to_csv(data_path+'/model_sea_decisiontreee.csv'.replace('.pickle','.csv'))

data_knn = pd.DataFrame(data_knn,columns=labels)
data_knn.to_csv(data_path+'/model_sea_knn.csv'.replace('.pickle','.csv'))

data_randomforest = pd.DataFrame(data_randomforest,columns=labels)
data_randomforest.to_csv(data_path+'/model_sea_randomforest.csv'.replace('.pickle','.csv'))

出力されたCSVは以下。

◆K-NN法

◆ランダムフォレスト

◆決定木

K-NN法の結果にのみ「０℃」以下の値が含まれているため、K-NNが最も実態に近い数値が含まれていると判断。今後の分析で利用する海水温データを、K-NN法のデータに決定。

#日本の漁獲高予測モデルを作成

さて、やっと本題に入る（この時点でお腹いっぱいだが、もう少し頑張ってみる事にした）。 分析する前に目的変数と、説明変数を整理。

目的変数の項目数

• 漁獲高：１個（魚種毎に１つ）

説明変数

• 日本海の海水温：100個（100エリア）
• 地球の気温推移：３個（全世界、北半球、南半球）
• 隣国の人口：3個（中国、韓国、北朝鮮）
• 日本の漁業人口：１個

合計：107個

この状態でモデリングすると、目的変数が多いため過学習となる。

説明変数を加工し数を減らす事にした。

加工後の説明変数

• 日本海の海水温：1個（100エリアの平均値を取る）
• 地球の気温推移：1個（全世界分のみ）
• 隣国の人口：1個（中国、韓国、北朝鮮の合計値）
• 日本の漁業人口：１個

合計：4個

目的変数（漁獲高）に対して、４つの説明変数でモデリングを行った。

モデリング

利用するモデルは、「LinearRegression」「Lasso」「Ridge」「ElasticNet」の４つ。各魚種の最もスコアの高いモデルを採用することにする。スコアが0以下のモデルは対象外とした。

import datetime
from email import header
import os
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
from sklearn.linear_model import ElasticNet
import matplotlib.pyplot as plt

#日本語表示用のフォントファミリーに変更
plt.rcParams['font.family'] = 'Hiragino Sans'
fig = plt.figure(figsize=(14,7))

def add_plot_fish(df,i,title):
  ax = fig.add_subplot(5,7,i)
  ax.set_title(title)
  ax.plot(df)

root_path = os.getcwd() +'/aidemy_blog_v2'
data_path = root_path +'/use_data'

<h2>#日本の漁獲量データを取得</h2>
fish_japan_all = pd.read_csv(data_path+'/jpn_fish.csv',index_col='year').fillna(0)
target_array =fish_japan_all.columns

#インデックス（年）の初期値と最終値を取得
first_year = fish_japan_all.index[0][0:4]
last_year = fish_japan_all.index[-1][0:4]
sd = datetime.date(int(first_year), 12, 31)
ed = datetime.date(int(last_year), 12, 31)
index = pd.date_range(sd,ed,freq='Y')

fish_japan_all.index = index
fish_japan_all.index.name = 'year'

#世界の平均気温上昇率を取得(first_year〜last_yearまで)
avg_temp_world = pd.read_csv(data_path+'/temp_world.csv',index_col='year').fillna(0)
avg_temp_world.index = index
avg_temp_world.index.name = 'year'

#日本の平均海水温(first_year〜last_yearまで)
avg_temp_jpn = pd.read_csv(data_path+'/temp_jpn_sea.csv',index_col='year').fillna(0)
avg_temp_jpn.index = index
avg_temp_jpn.index.name = 'year'

#隣国の人口推移を取得（first_year〜last_year：ロシアのぞく）
population = pd.read_csv(data_path+'/sum_pop_naver.csv',index_col='year').fillna(0)
population.index = index
population.index.name = 'year'

#日本の漁業人口推移を取得(first_year〜last_year)
fisherman = pd.read_csv(data_path+'/sarima_fisherman.csv',index_col='year').fillna(0)
fisherman.index = index
fisherman.index.name = 'year'

#■全データを説明変数にした場合ー説明変数のDataFrameをマージ
df_X_values = pd.merge(avg_temp_jpn,fisherman, on='year',how='outer')
df_X_values = pd.merge(df_X_values,avg_temp_world, on='year',how='outer')

#１魚種ずつモデルを作成する
values=[]
models={}

fig = plt.figure(figsize=(14,7))
plt.subplots_adjust(wspace=0.4, hspace=0.6)

for fish_name in target_array:
  val =[fish_name]
  #１魚種に絞ったDataFrameを取得（目的変数の部分)
  fish_japan_solo = fish_japan_all[fish_name]
  #各説明変数のDataFrameと結合
  data = pd.merge(fish_japan_solo,df_X_values, on='year',how='outer')

  #学習データと検証データに分ける
  X = data.values[:,1:]
  y = data.values[:,0]
  X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0, shuffle=True)

  best_value=0
  #LinearRegression,Lasso,Ride,ElasticNetで魚種毎にもっとスコアの高いモデルを選定する
  for model in [LinearRegression(),Lasso(),Ridge(),ElasticNet(l1_ratio=0.3)]:
      # fit関数を使って、モデルを構築してください。
      model.fit(X_train, y_train)
      ms = model.score(X_test, y_test)
      #スコアの高いモデルをmodelsに保持しておく
      if(best_value < ms):
        models[fish_name]=model
        best_value = ms

  values.append("《" + fish_name + "》-----------------------------------")
  values.append("--Method:" + models[fish_name].__class__.__name__ + "--")
  values.append("★Cross validatin scores:{}".format(best_value))

#作成したモデルで各魚種の予測をしてみる
c = 1
for fn,mdl in models.items():
  test_data = fish_japan_all[fn]
  pred_data = pd.DataFrame(mdl.predict(df_X_values))
  pred_data.index = index
  pred_data.columns =['pred-'+fn]
  df_plt = pd.merge(test_data,pred_data,on='year',how='outer')
  add_plot_fish(df_plt,c,fn)  
  c += 1

#結果の出力
df = pd.DataFrame(values)
df.to_csv(data_path+'/move_values_01.csv')
plt.savefig(data_path+'/pred_4.png')
print(len(models))

各魚種の最も高いスコアが出たモデルと、スコアは以下となった。

0,《合計》-----------------------------------
1,--Method:ElasticNet--
2,★Cross validatin scores:0.8636874598610523
3,《くろまぐろ》-----------------------------------
4,--Method:ElasticNet--
5,★Cross validatin scores:0.5080198088433892
6,《みなみまぐろ》-----------------------------------
7,--Method:LinearRegression--
8,★Cross validatin scores:0.3332223605886764
9,《めばち》-----------------------------------
10,--Method:ElasticNet--
11,★Cross validatin scores:0.837618717810202
12,《きはだ》-----------------------------------
13,--Method:Lasso--
14,★Cross validatin scores:0.4941346469628283
15,《その他のまぐろ類》-----------------------------------
16,--Method:Ridge--
17,★Cross validatin scores:0.7662042603490956
18,《まかじき》-----------------------------------
19,--Method:ElasticNet--
20,★Cross validatin scores:0.5867491850217903
21,《めかじき》-----------------------------------
22,--Method:LinearRegression--
23,★Cross validatin scores:0.5648103175077084
24,《くろかじき類》-----------------------------------
25,--Method:LinearRegression--
26,★Cross validatin scores:0.8803790847860316
27,《その他のかじき類》-----------------------------------
28,--Method:LinearRegression--
29,★Cross validatin scores:0.38190203917955146
30,《かつお》-----------------------------------
31,--Method:ElasticNet--
32,★Cross validatin scores:0.6821408337047192
33,《そうだがつお類》-----------------------------------
34,--Method:ElasticNet--
35,★Cross validatin scores:0.3739849715955603
36,《さめ類》-----------------------------------
37,--Method:LinearRegression--
38,★Cross validatin scores:0.355257860468284
39,《さけ類》-----------------------------------
40,--Method:Ridge--
41,★Cross validatin scores:0.6380437924773816
42,《ます類》-----------------------------------
43,--Method:Ridge--
44,★Cross validatin scores:0.7068615643579998
45,《このしろ》-----------------------------------
46,--Method:Ridge--
47,★Cross validatin scores:0.3329546063792177
48,《にしん》-----------------------------------
49,--Method:LinearRegression--
50,★Cross validatin scores:0.6651980672584288
51,《まいわし》-----------------------------------
52,--Method:Ridge--
53,★Cross validatin scores:0.7141038828178389
54,《かたくちいわし》-----------------------------------
55,--Method:ElasticNet--
56,★Cross validatin scores:0.5611845547187753
57,《まあじ》-----------------------------------
58,--Method:Lasso--
59,★Cross validatin scores:0.848520276853194
60,《むろあじ類》-----------------------------------
61,--Method:ElasticNet--
62,★Cross validatin scores:0.8735304202447156
63,《さんま》-----------------------------------
64,--Method:LinearRegression--
65,★Cross validatin scores:0.7834346104671887
66,《かれい類》-----------------------------------
67,--Method:LinearRegression--
68,★Cross validatin scores:0.9369443682329953
69,《すけとうだら》-----------------------------------
70,--Method:LinearRegression--
71,★Cross validatin scores:0.5058699736924356
72,《ほっけ》-----------------------------------
73,--Method:ElasticNet--
74,★Cross validatin scores:0.6005292012577094
75,《きちじ》-----------------------------------
76,--Method:Lasso--
77,★Cross validatin scores:0.5867116094444853
78,《にぎす類》-----------------------------------
79,--Method:LinearRegression--
80,★Cross validatin scores:0.9386541319488325
81,《あなご類》-----------------------------------
82,--Method:ElasticNet--
83,★Cross validatin scores:0.9314813533372198
84,《たちうお》-----------------------------------
85,--Method:ElasticNet--
86,★Cross validatin scores:0.8882533705562337
87,《いさき》-----------------------------------
88,--Method:Lasso--
89,★Cross validatin scores:0.8091891771258445
90,《すずき類》-----------------------------------
91,--Method:Ridge--
92,★Cross validatin scores:0.7177184115110604
93,《あまだい類》-----------------------------------
94,--Method:LinearRegression--
95,★Cross validatin scores:0.6496018877535136
96,《ふぐ類》-----------------------------------
97,--Method:LinearRegression--
98,★Cross validatin scores:0.6528189548975889
99,《その他の魚類》-----------------------------------
100,--Method:LinearRegression--
101,★Cross validatin scores:0.8996113881308885

モデルを利用した予測結果と、元データのグラフの比較は以下。

※ブルーのライン：元データ

※オレンジのライン：モデルの予測結果

#まとめ

魚種によっては精度の高い予測値が出ている事が分かる。
スコア0.8を超えるものもあるため、モデリングは上手くいった様子。

#感想

１変量の時系列解析、教師あり学習（回帰）、教師あり学習（分類）の３つをこの課題を通して実践できた。教材を何度も読み返し、数式の理解を深め、分析手法の使い所などを体験する良い機会となった。今後も学習を継続していきたい。