コロナが来なかった場合の国内映画興行収入予測

0. 目次

1. はじめに
2. データ内容
3. 分析の狙いと推測
4. 実行環境
5. データ分析の流れ
6. データ分析の結果
7. 終わりに

1. はじめに

昨今、様々な動画配信サービスの登場により、契約してしまえば映画館にもレンタルビデオ店にも行かず、自宅でドラマや映画が見放題という数年前までは考えられないような世界になった。
私自身も「NetFilx」「Disney+」「Prime Video」など複数のサービスを利用し、楽しんでいる。
これらにサービスを利用してる人なら1度は私と同じ経験をしたことがあるかもしれない。

最近まで上映してたはずなのに、もう配信が始まってる！！

すぐ配信するなら映画館に行かずに配信待てばよかった．．．

大迫力のスクリーン・高品質のサウンド・様々な映画予告を見ながら上映を待つワクワク感は、やはり映画館に足を運ばなければ体験出来ないが、自宅で気軽に映画が見られる時代となってからは映画館に行く機会は減り、最後に映画館に行ったのはいつか思い出せないほどになってしまった。

前置きが長くなってしまったが、「一般社団法人 日本映画製作者連盟」 の過去データ一覧表を利用させていただき、67年間分の映画館の軌跡をデータで確認し、データ分析を行ってみたいと思う。

2. データ内容

過去データ一覧表は、1955年～2021年までの67年間分の年度別データ7項目が下記の内容で公開されている。

• 西暦
• 映画館数 (スクリーン数)
• [邦画][洋画][合計]の公開本数
• 入場者数
• 興行収入
• 平均料金
• [邦画][洋画][合計]の配給収入
• [邦画][洋画]のシェア

配給収入は1999年で公開が終了しており、以降は興行収入のみの公開となっているため、配給収入は分析に使用せず、興行収入の合計を使用する。
2000年からは映画館の項目が映画館スクリーン数に変わっており、1999年以前は映画館数＝スクリーン数とカウントしていたが、2000年以降は1つの映画館が複数のスクリーンを所持しているシネマコンプレックスが主流になってきた状況を反映したものだと推測する。

3. 分析の狙いと推測

各パラメータの傾向を折れ線グラフで視覚化し、映画館数や公開本数、入場者数などから映画館の盛衰を確認し、予測に役立てたい。
データ分析を行う前の簡単な推測としては、【動画配信サービスの台頭】【新型コロナウィルスの流行】【SNSやYoutubeなどの娯楽の増加】【映画館料金増加】などの理由により、映画館の数や興行収入、入場者数は年々減っている可能性が高いと思う。

4. 実行環境

(1) PCスペック

• HP ENVY x360 Convertible 15-ee0-xxx
• AMD Ryzen 7 4700U with Radeon Graphics
• RAM 16GB
• Windows 11 Home Edition

(2) プログラミング言語

• Python3 (Ver 3.11.1)
• JupyterNotebook

(3) ライブラリ

• Numpay
• Pandas
• Matplotlib
• Scipy
• Seaborn
• Scikit-learn
• StatsModels

5.データ分析の流れ

[1]「一般社団法人日本映画製作者連盟 過去データ一覧表」をExcelにコピーし、余分なデータを削除、項目の整理を行い、CSV UTF-8(コンマ切)で保存。

CSVデータ内の表示形式に注意。
一部データのExcel表示形式が「通貨」になっており、文字列データだったため、修正するまで上手くデータが使用できなかった。

[2] Pandasでデータを読み込み、各データ毎の折れ線グラフを作成し、グラフの可視化によりデータにどのような傾向があるか確認。

qiita.py

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

pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', 10)
df = pd.read_csv('cinema_jp_data.csv')
display(df)

qiita.py

#映画館数のグラフ化--------------------------------------------------------------------
plt.subplots(figsize=(15,7))
plt.title('映画館数(スクリーン)', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(1500, 8000, step=500))
plt.grid(True)

plt.plot(df['西暦'], df['映画館数 (スクリーン)'], linestyle='-', color='k', marker='o', markerfacecolor="k", markersize=4)
plt.show()


#公開本数のグラフ化--------------------------------------------------------------------
plt.subplots(figsize=(15,7))

plt.title('邦画・洋画別の公開本数', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(100, 750, step=50))
plt.grid(True)

plt.plot(df['西暦'], df['邦画 (本)'], linestyle='-', color='r', marker='o', markerfacecolor="r", markersize=4, label='邦画')
plt.plot(df['西暦'], df['洋画 (本)'], linestyle='-', color='b', marker='o', markerfacecolor="b", markersize=4, label='洋画')
plt.legend(loc='upper left', prop={"family":"MS Gothic"})
plt.show()


plt.subplots(figsize=(15,7))
plt.title('邦画・洋画の合計公開本数', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(450, 1350, step=50))
plt.grid(True)

plt.plot(df['西暦'], df['合計 (本)'], linestyle='-', color='g', marker='o', markerfacecolor="g",markersize=4)
plt.show()


#入場者数のグラフ化--------------------------------------------------------------------
plt.subplots(figsize=(15,7))
plt.title('入場者数(単位:万人)', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(5000, 140000, step=5000))
plt.grid(True)

plt.plot(df['西暦'], df['入場者数 (万人)'], linestyle='-', color='k', marker='o', markerfacecolor="k", markersize=4)
plt.show()


#興行収入のグラフ化--------------------------------------------------------------------
plt.subplots(figsize=(15,7))
plt.title('興行収入(単位:百万円)', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(45000, 300000, step=10000))
plt.grid(True)

plt.plot(df['西暦'], df['興行収入 (百万円)'], linestyle='-', color='k', marker='o', markerfacecolor="k", markersize=4)
plt.show()


#平均料金のグラフ化--------------------------------------------------------------------
plt.subplots(figsize=(15,7))
plt.title('平均料金(円)', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(0, 1500, step=100))
plt.grid(True)

plt.plot(df['西暦'], df['平均料金 (円)'], linestyle='-', color='k', marker='o', markerfacecolor="k", markersize=4)
plt.show()


#邦画・洋画のシェアグラフ化---------------------------------------------------------------
plt.subplots(figsize=(15,7))
plt.title('邦画・洋画別のシェア(%)', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(0, 100, step=5))
plt.grid(True)

plt.plot(df['西暦'], df['邦画率 (%)'], linestyle='-', color='r', marker='o', markerfacecolor="r", markersize=4, label='邦画率')
plt.plot(df['西暦'], df['洋画率 (%)'], linestyle='-', color='b', marker='o', markerfacecolor="b", markersize=4, label='洋画率')
plt.legend(loc='upper left', prop={"family":"MS Gothic"})
plt.show()

グラフを確認して気になったこと。

映画館数と入場者数のグラフ傾向と興行収入の傾向が異なっていることに違和感を覚えた。
1955年と2021年の円の価値には差があり、興行収入や平均料金が上昇傾向なのは当たり前なのではないか。と考えたため、「日本銀行HP」記載の消費者物価指数を用いて、2021年の消費者物価指数をベースに円価値を変換して再度グラフ化してみたいと思う。

計算例
消費者物価指数(持家の帰属家賃を除く総合)　99.7(2021年) ÷ 16.6(1955年) = 6.007倍

qiita.py

df_money = pd.read_csv('cinema_money_data.csv')
display(df_money)

qiita.py

#消費者物価指数変換前と後の平均料金のグラフ化-------------------------------------------------
plt.subplots(figsize=(15,7))
plt.title('消費者物価指数変換前と後の平均料金(円)', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(0, 1500, step=100))
plt.grid(True)

plt.plot(df_money['西暦'], df_money['原：平均料金 (円)'], linestyle='-', color='k', marker='o', markerfacecolor="k", markersize=4, label='原：平均料金')
plt.plot(df_money['西暦'], df_money['変：平均料金 (円)'], linestyle='-', color='y', marker='o', markerfacecolor="y", markersize=4, label='変：平均料金')
plt.legend(loc='upper left', prop={"family":"MS Gothic"})
plt.show()


#消費者物価指数変換前と後の平均料金のグラフ化-------------------------------------------------
plt.subplots(figsize=(15,7))
plt.title('消費者物価指数変換前と後の興行収入(単位:百万円)', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(45000, 600000, step=10000))
plt.grid(True)

plt.plot(df_money['西暦'], df_money['原：興行収入 (百万円)'], linestyle='-', color='k', marker='o', markerfacecolor="k", markersize=4, label='原：興行収入')
plt.plot(df_money['西暦'], df_money['変：興行収入 (百万円)'], linestyle='-', color='y', marker='o', markerfacecolor="y", markersize=4, label='変：興行収入')
plt.legend(loc='upper right', prop={"family":"MS Gothic"})
plt.show()

消費者物価指数を用いた変換を行ったことで、変換前と後で興行収入のグラフの傾向に大きな違いが出来た。
映画館数や入場者数のグラフとも形状が類似しているため、参考程度の変換ではあるが納得できる結果を得られたと思う。

[3] 今回のデータは時系列のため、SARIMAモデルを使用して時系列解析を行ってみる。

SARIMAモデルを用いた時系列データの分析手順は以下の通り
　1. データの読み込み
　2. データの整理　　←ここから始める。
　3. データの可視化
　4. データの周期の把握 (パラメーターsの決定)
　5. s以外のパラメーターの決定
　6. モデルの構築
　7. データとの予測とその可視化

SARIMAモデルとは？
SARIMAモデルとはARIMAモデルをさらに季節周期を持つ時系列データにも拡張できるようにしたモデル。
ARIMAモデルがARMAモデルへ原系列を階差系列に変換し適応したもので、ARMAモデルは定常過程にしか適応できませんがARIMAモデルは非定常過程（データの平均や分散が時間に依存している過程）にも適応可能。

SARIMAモデルの説明がされているURLを添付。
定常時系列の解析に使われるARMAモデル・SARIMAモデルとは？

2. データ整理
もしコロナが来なかった場合の映画館興行収入について予測してみたいと思うので、使用する時系列データは1955~2019年までに整理する。

qiita.py

drop_df = df.drop(['映画館数 (スクリーン)', '邦画 (本)', '洋画 (本)', '合計 (本)',
                   '平均料金 (円)', '邦画率 (%)', '洋画率 (%)', '入場者数 (万人)'], axis=1)
drop_df['年度'] = drop_df['西暦'].astype(str)
drop_df['年度'] = pd.DatetimeIndex(drop_df['年度']).year
drop_df = drop_df.drop(['西暦'], axis=1)
grouped = drop_df.groupby('年度')
df_pi = grouped.sum()
df_pi.add_suffix('').reset_index()
df_pi = df_pi.iloc[:65]
display(df_pi)

3. データの可視化は上記で行っているためコードは割愛。

4. データの周期の把握 (パラメーターsの決定)
データの周期の把握(パラメーターsの決定)については、偏自己相関係数を求め可視化して決定する。

qiita.py

fig = plt.figure(figsize=(15,5))

#自己相関の可視化
pi_acf = sm.tsa.stattools.acf(df_pi, nlags=30, fft=False)
ax1 = fig.add_subplot(1,2,1)
fig = sm.graphics.tsa.plot_acf(df_pi, lags=30, ax=ax1)

#偏自己相関の可視化
pi_pacf = sm.tsa.stattools.pacf(df_pi, nlags=30)
ax2 = fig.add_subplot(1,2,2)
fig = sm.graphics.tsa.plot_pacf(df_pi, lags=30, ax=ax2)
plt.show()

上図よりパラメータ s=6 で検証を行ってみたい。
また、sm.tsa.seasonal_decompose()を使用して[原系]、[トレンド]、 [季節変動]、 [残差]のデータパターンを確認する。

qiita.py

plt.rcParams['figure.figsize'] = [10, 10]

fig = sm.tsa.seasonal_decompose(df_pi, period=6).plot()
plt.show()

5. s以外のパラメーターの決定
PythonにはSARIMAモデルのパラメーター(p, d, q)(sp, sd, sq, s)を自動で最適化してくれる機能がないため、情報量規準(今回の場合は BIC(ベイズ情報量基準))によってどの値が最も適切なのか調べるプログラムを書き、最も良いパラメーターとそのBICを出力します。

qiita.py

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 print(parameters[np.argmin(BICs)])

selectparameter(df_pi, s=6)
>[(0, 1, 0), (0, 1, 1, 6), 1272.532880952859]

上記の[(p,d,q),(sp,sd,sq,s)]の数値を採用します。

6. モデルの構築
モデルの構築にはsm.tsa.statespace.SARIMAX(DATA,order=(p, d, q),seasonal_order=(sp, sd, sq, s)).fit()を用います。

qiita.py

SARIMA_performanceincome = sm.tsa.statespace.SARIMAX(df_pi,order=(0,1,0),seasonal_order=(0,1,1,6)).fit() 
print(SARIMA_performanceincome.bic)
>1272.532880952859

print(SARIMA_performanceincome.summary())

BICが同じことを確認し、結果は以下の通り。

qiita.py

                                     SARIMAX Results                                      
===========================================================================================
Dep. Variable:                          興行収入 (百万円)   No. Observations:                   65
Model:             SARIMAX(0, 1, 0)x(0, 1, [1], 6)   Log Likelihood                -632.206
Date:                             Wed, 11 Jan 2023   AIC                           1268.412
Time:                                     15:42:22   BIC                           1272.533
Sample:                                          0   HQIC                          1270.017
                                              - 65                                         
Covariance Type:                               opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ma.S.L6       -0.7483      0.174     -4.289      0.000      -1.090      -0.406
sigma2      2.319e+08   4.28e-11   5.42e+18      0.000    2.32e+08    2.32e+08
===================================================================================
Ljung-Box (L1) (Q):                   1.44   Jarque-Bera (JB):                 0.74
Prob(Q):                              0.23   Prob(JB):                         0.69
Heteroskedasticity (H):               3.06   Skew:                             0.13
Prob(H) (two-sided):                  0.02   Kurtosis:                         3.48
===================================================================================

7. データとの予測とその可視化
1975年から2021年までの予測をしてみる。

qiita.py

pred = SARIMA_performanceincome.predict(1975, 2021)
plt.plot(pred)
plt.show()

元の時系列データと予測データを同時出力して予測が正しいか確認してみる。(赤線が予測データ）

qiita.py

plt.plot(df_pi)
plt.plot(pred, color="r")
plt.show()

予測データが右肩上がりになっていたため、正しく学習出来ていないとは思ったが、元データと同時出力したことで明確になった。

原因予想としては、
「データ数が65個と少なすぎたため」
「年度別データだったため」
「季節周期がなくSARIMAモデルでの解析に向ていなかった」
「上記3つの理由ではなく、データ整理から予測までの流れの中で間違っている部分がある」
などが考えられる。

アプローチ方法をSARIMAから回帰分析に変更して再度予測してみたいと思う。

[4] 回帰分析を使って予測を行う。

回帰分析の手順は以下の通り
　1. データの読み込み
　2. データの整理
　3. 目的変数の決定
　4. 目的変数に影響を与えていそうな説明変数の決定
　5. 学習モデルの構築
　6. データとの予測とその可視化

回帰分析とは？
求めたい要素(目的変数)の値に対し、他要素(説明変数)の影響度を分析する手法。
説明変数の個数により、2種類に分けられる。
説明変数1つ　→　単回帰分析
説明変数複数　→　重回帰分析

1. データの読み込み

qiita.py

pip install japanize-matplotlib

qiita.py

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
from sklearn.linear_model import ElasticNet
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_squared_error
from scipy.stats import pearsonr
from scipy import stats
pd.set_option("display.max_columns", None)
pd.set_option("display.max_rows", 10)

qiita.py

df = pd.read_csv('cinema_jp_data.csv')
money_df = pd.read_csv('cinema_money_data.csv')

2. データの整理
興行収入と平均料金は消費者物価指数を元に変換した値に変更する。
2020年と2021年は除いたデータとする。

qiita.py

df = df.drop(['西暦', '平均料金 (円)', '興行収入 (百万円)'], axis=1)
money_df= money_df.drop(columns= money_df.columns[range(0, 5)], axis=1)

df_add = pd.concat([df, money_df], axis=1)
cinema_df = df_add.rename(columns={'変：平均料金 (円)': '平均料金 (円)', '変：興行収入 (百万円)': '興行収入 (百万円)'}).drop(index= df_add.index[65:], axis=0)
display(cinema_df)

3. 目的変数の決定
今回は「コロナが来なかった場合の興行収入」を予測しようと思うので、興行収入を目的変数に選ぶ。

4. 目的変数に影響を与えていそうな説明変数の決定
まず各パラメータの相関を確認する。

qiita.py

df_pearson = cinema_df.corr(method='pearson')

fig, ax = plt.subplots(figsize=(10, 8))
sns.heatmap(df_pearson, square=True, annot=True, fmt='.2f',annot_kws={'size':10},vmax=1, vmin=-1, center=0)
plt.title('Correlation coefficient (pearson)',fontsize=15)
plt.show()

「映画館数」「入場者数」「平均料金」 とは強い相関があることが確認できるが、興行収入に影響を与えていることは明白なので、その他で考える。
今回は、映画本数に着目して予測してみようと思うので、興行収入と映画本数について無相関検定を行ってみる。
有意水準は0.05を定める。

qiita.py

performance_income = cinema_df['興行収入 (百万円)'].values
total_films = cinema_df['合計 (本)'].values
japanese_films = cinema_df['邦画 (本)'].values
foreign_films = cinema_df['洋画 (本)'].values

result_1 = pearsonr(performance_income, total_films)
pi_value = result_1[0]
tf_value = result_1[1]
print('相関係数：', pi_value)
print('p値：', tf_value)
> 相関係数：-0.0030337689438222487
> p値：0.9808649816603856

result_2 = pearsonr(performance_income, japanese_films)
pi_value = result_2[0]
jf_value = result_2[1]
print('相関係数：', pi_value)
print('p値：', jf_value)
> 相関係数：0.46667079635116704
> p値：8.911616322683429e-05

result_3 = pearsonr(performance_income, foreign_films)
pi_value = result_3[0]
ft_value = result_3[1]
print('相関係数：', pi_value)
print('p値：', ft_value)
> 相関係数：-0.46415681994620006
> p値：9.837836149056026e-05

result_1より興行収入と合計本数は、
p値 > 有意水準（0.05）となり、2変数間に相関がない (相関係数が0である)という仮説は棄却されない。

result_2より興行収入と邦画本数は、
p値 < 有意水準（0.05）であるため帰無仮説は棄却され、二つの変数は無相関ではないだろうと考えられる。

result_3より興行収入と洋画本数は、
p値 < 有意水準（0.05）であるため帰無仮説は棄却され、二つの変数は無相関ではないだろうと考えられる。

5. 学習モデルの構築
複数の学習モデル(線形回帰、ラッソ回帰、リッジ回帰、ElasticNet回帰、ランダムフォレスト回帰)で比較し、
決定係数が最も高いモデルを用いる。

qiita.py

df_x = cinema_df.drop(cinema_df.columns[[0,3,4,5,6,7,8]].values, axis=1)
df_y = cinema_df['興行収入 (百万円)'].values

train_X, test_X, train_y, test_y = train_test_split(df_x, df_y, test_size = 0.2, random_state=42)
print(train_X.shape, test_X.shape, train_y.shape, test_y.shape)
>(52, 2) (13, 2) (52,) (13,)

qiita.py

max_score = 0 
best_model = ""
RMSE = 0

# 線形回帰
model_1 = LinearRegression()
model_1.fit(train_X, train_y)
pred_y_lr = model_1.predict(test_X)
score = model_1.score(test_X, test_y)
rmse_1 = mean_squared_error(test_y, pred_y_lr, squared=False)
if max_score < score: 
    max_score = score 
    best_model = "LinearRegression" 
    RMSE = rmse_1
    
print('LinearRegression Result')
print('R-squared:', model_1.score(test_X, test_y))
print('RMSE:', rmse_1)
#print(y_pred_lr)
print()

# ラッソ回帰
model_2 = Lasso()
model_2.fit(train_X, train_y)
pred_y_lasso = model_2.predict(test_X)
score = model_2.score(test_X, test_y)
rmse_2 = mean_squared_error(test_y, pred_y_lasso, squared=False)
if max_score < score:
    max_score = score
    best_model = "LassoRegression"
    RMSE = rmse_2
    
print('LassoRegression Result')
print('R-squared:', model_2.score(test_X, test_y))
print('RMSE:', rmse_2)
#print(y_pred_lasso)
print()

# リッジ回帰
model_3 = Ridge()
model_3.fit(train_X, train_y)
pred_y_ridge = model_3.predict(test_X)
score = model_3.score(test_X, test_y)
rmse_3 = mean_squared_error(test_y, pred_y_ridge, squared=False)
if max_score < score:
    max_score = score
    best_model = "RidgeRegression"
    RMSE = rmse_3

print('RidgeRegression Result')
print('R-squared:', model_3.score(test_X, test_y))
print('RMSE:', rmse_3)
#print(y_pred_ridge)
print()
    
# ElasticNet回帰
model_4 = ElasticNet(l1_ratio=0.5)
model_4.fit(train_X, train_y)
pred_y_en = model_4.predict(test_X)
score = model_4.score(test_X, test_y)
rmse_4 = mean_squared_error(test_y, pred_y_en, squared=False)
if max_score < score:
    max_score = score
    best_model = "ElasticNetRegression"
    RMSE = rmse_4

print('ElasticNetRegression Result')
print('R-squared:', model_4.score(test_X, test_y))
print('RMSE:', rmse_4)
#print(y_pred_en)
print()

# ランダムフォレスト回帰
parameters = {  
    'n_estimators': [10, 20, 30, 50, 100, 300],
    'max_features': ('sqrt', 'log2', None),
    'max_depth':    (10, 20, 30, 40, 50, None),
    'random_state': [42]
}
model = RandomForestRegressor()
gridsearch = GridSearchCV(estimator = model, param_grid = parameters)
gridsearch.fit(train_X, train_y)

model_5 = RandomForestRegressor(n_estimators=gridsearch.best_params_['n_estimators'], 
                              max_features=gridsearch.best_params_['max_features'], 
                              max_depth=gridsearch.best_params_['max_depth'], 
                              random_state=gridsearch.best_params_['random_state'])
model_5.fit(train_X, train_y)
pred_y_rfr = model_5.predict(test_X)
score = model_5.score(test_X, test_y)
rmse_5 = mean_squared_error(test_y, pred_y_rfr, squared=False)
if max_score < score:
    max_score = score
    best_model = "RandomForestRegressor"
    RMSE = rmse_5
    
print('RandomForestRegression Result')
print('R-squared:', model_5.score(test_X, test_y))
print('RMSE:', rmse_5)
#print(pred_y_rfr)

qiita.py

LinearRegression Result
R-squared: 0.6473896210941886
RMSE: 41843.00044997249

LassoRegression Result
R-squared: 0.647389583506369
RMSE: 41843.00268017785

RidgeRegression Result
R-squared: 0.6473893478051902
RMSE: 41843.01666507867

ElasticNetRegression Result
R-squared: 0.6473824956712984
RMSE: 41843.423221939556

RandomForestRegression Result
R-squared: 0.7458478088424674
RMSE: 35524.02790435272

qiita.py

print("BestModel:{}".format(best_model))
print("R-squared:{}".format(max_score))
print("RMSE:{}".format(RMSE))

> BestModel:RandomForestRegressor
> R-squared:0.7458478088424674
> RMSE:35524.02790435272

評価指標が決定係数のみだとモデルの評価が0～1の範囲であるため、直感的にわかりやすい一方で、
具体的なモデルの精度がわからないため、予測範囲を言いづらい側面がある。
そのため、二乗平均平方根誤差(RMSE)も評価指標とすることで、「±いくら」で予測できるため、体感的にわかりやすくなった。

6.データ分析の結果

5つの学習モデルで検証した結果、ランダムフォレスト回帰の精度が最も良いという答えが出た。
決定係数： 0.75
二乗平均平方根誤差: ± 355億2400万円
予測金額の誤差の幅がありすぎて実用的ではないと判断するが、2020年と2021年の映画本数で実値と比較したいと思う。

2020年　　邦画：506 本 　洋画：511 本　　興行収入： 1428億5600万円

qiita.py

model_5.predict(np.array([[506, 511]]))
> array([221785.4])

結果：2217億8540万円
誤差：789億2940万円

2021年　　邦画：490 本 　洋画：469 本　　興行収入： 1618億9300万円

qiita.py

model_5.predict(np.array([[490, 469]]))
> array([222915.7])

結果：2229億1570万円
誤差：612億2270万円

予測結果と実際の値をグラフで視覚化する。
cinema_df2 が予測値を追加したもの
cinema_df3 が実値を追加したもの

qiita.py

cinema_df2 = cinema_df['興行収入 (百万円)']
cinema_df2 = cinema_df2.append(pd.Series([221785, 222915]))
cinema_df2.index = range(1955, 2022)

cinema_df3 = cinema_df['興行収入 (百万円)']
cinema_df3 = cinema_df3.append(pd.Series([142856, 161893]))
cinema_df3.index = range(1955, 2022)

qiita.py

fig = plt.subplots(figsize=(15,7))

plt.title('興行収入の予測値と実値', fontname="MS Gothic", fontsize=15)
plt.xlabel('西暦', fontname="MS Gothic")
plt.xticks(np.arange(1955, 2021, step=5))
plt.yticks(np.arange(45000, 600000, step=10000))
plt.grid(True)

plt.plot(cinema_df2, linestyle='-', color='k', marker='o', markerfacecolor="k", markersize=4, label='予測値')
plt.plot(cinema_df3, linestyle='-', color='y', marker='o', markerfacecolor="y", markersize=4, label='実値')
plt.legend(loc='upper right', prop={"family":"MS Gothic"})

plt.show()

結果の推察
二乗平均平方根誤差が ±355億2400万円 に対して、実値の誤差が大きく異なった。
これは学習モデルがコロナ影響前の2019年までのデータでの学習結果であるため、実際はコロナ感染拡大の影響で予想以上に興行収入が落ち込んでしまったことが原因ではないかと思う。
学習モデルの精度問題はあるが、コロナの影響を大きく受けていることをデータで確認し、改めて映画館の今後が心配になってしまった。
近年は邦画アニメが好調の様で、映画1本で興行収入100億円を超えるような作品もあるため、今後興行収入がコロナ前に戻っていくことを願いたい。

7. 終わりに

自分の中でまだまだ理解出来ていないことが多く、間違っているところもあると思うが、自分の興味があるデータを取得し、一連のデータ分析の流れをアウトプット出来たことは良い勉強になった。
SARIMAモデルが使用できず、回帰分析に手法を変え結果出力まで至ったが、データをよく確認し理解することが重要だということを思い知った。

3ヶ月前まではPythonでコードの書き方すらわからなかった状態でも、学習を続けることでなんとか形にすることが出来たので、今後も学習を続けることで地道にでも成長していこうと思う。

あと、上映が終わる前に 『Dr.コトー診療所』 を観に行って、興行収入に貢献したい。